src/LAPACK/dsteqr.f

   1       SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          COMPZ
   9       INTEGER            INFO, LDZ, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  19 *  symmetric tridiagonal matrix using the implicit QL or QR method.
  20 *  The eigenvectors of a full or band symmetric matrix can also be found
  21 *  if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
  22 *  tridiagonal form.
  23 *
  24 *  Arguments
  25 *  =========
  26 *
  27 *  COMPZ   (input) CHARACTER*1
  28 *          = 'N':  Compute eigenvalues only.
  29 *          = 'V':  Compute eigenvalues and eigenvectors of the original
  30 *                  symmetric matrix.  On entry, Z must contain the
  31 *                  orthogonal matrix used to reduce the original matrix
  32 *                  to tridiagonal form.
  33 *          = 'I':  Compute eigenvalues and eigenvectors of the
  34 *                  tridiagonal matrix.  Z is initialized to the identity
  35 *                  matrix.
  36 *
  37 *  N       (input) INTEGER
  38 *          The order of the matrix.  N >= 0.
  39 *
  40 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
  41 *          On entry, the diagonal elements of the tridiagonal matrix.
  42 *          On exit, if INFO = 0, the eigenvalues in ascending order.
  43 *
  44 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
  45 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
  46 *          matrix.
  47 *          On exit, E has been destroyed.
  48 *
  49 *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
  50 *          On entry, if  COMPZ = 'V', then Z contains the orthogonal
  51 *          matrix used in the reduction to tridiagonal form.
  52 *          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
  53 *          orthonormal eigenvectors of the original symmetric matrix,
  54 *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  55 *          of the symmetric tridiagonal matrix.
  56 *          If COMPZ = 'N', then Z is not referenced.
  57 *
  58 *  LDZ     (input) INTEGER
  59 *          The leading dimension of the array Z.  LDZ >= 1, and if
  60 *          eigenvectors are desired, then  LDZ >= max(1,N).
  61 *
  62 *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
  63 *          If COMPZ = 'N', then WORK is not referenced.
  64 *
  65 *  INFO    (output) INTEGER
  66 *          = 0:  successful exit
  67 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  68 *          > 0:  the algorithm has failed to find all the eigenvalues in
  69 *                a total of 30*N iterations; if INFO = i, then i
  70 *                elements of E have not converged to zero; on exit, D
  71 *                and E contain the elements of a symmetric tridiagonal
  72 *                matrix which is orthogonally similar to the original
  73 *                matrix.
  74 *
  75 *  =====================================================================
  76 *
  77 *     .. Parameters ..
  78       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
  79       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  80      $                   THREE = 3.0D0 )
  81       INTEGER            MAXIT
  82       PARAMETER          ( MAXIT = 30 )
  83 *     ..
  84 *     .. Local Scalars ..
  85       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
  86      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
  87      $                   NM1, NMAXIT
  88       DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
  89      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
  90 *     ..
  91 *     .. External Functions ..
  92       LOGICAL            LSAME
  93       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
  94       EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
  95 *     ..
  96 *     .. External Subroutines ..
  97       EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR,
  98      $                   DLASRT, DSWAP, XERBLA
  99 *     ..
 100 *     .. Intrinsic Functions ..
 101       INTRINSIC          ABS, MAX, SIGN, SQRT
 102 *     ..
 103 *     .. Executable Statements ..
 104 *
 105 *     Test the input parameters.
 106 *
 107       INFO = 0
 108 *
 109       IF( LSAME( COMPZ, 'N' ) ) THEN
 110          ICOMPZ = 0
 111       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 112          ICOMPZ = 1
 113       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 114          ICOMPZ = 2
 115       ELSE
 116          ICOMPZ = -1
 117       END IF
 118       IF( ICOMPZ.LT.0 ) THEN
 119          INFO = -1
 120       ELSE IF( N.LT.0 ) THEN
 121          INFO = -2
 122       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
 123      $         N ) ) ) THEN
 124          INFO = -6
 125       END IF
 126       IF( INFO.NE.0 ) THEN
 127          CALL XERBLA( 'DSTEQR', -INFO )
 128          RETURN
 129       END IF
 130 *
 131 *     Quick return if possible
 132 *
 133       IF( N.EQ.0 )
 134      $   RETURN
 135 *
 136       IF( N.EQ.1 ) THEN
 137          IF( ICOMPZ.EQ.2 )
 138      $      Z( 1, 1 ) = ONE
 139          RETURN
 140       END IF
 141 *
 142 *     Determine the unit roundoff and over/underflow thresholds.
 143 *
 144       EPS = DLAMCH( 'E' )
 145       EPS2 = EPS**2
 146       SAFMIN = DLAMCH( 'S' )
 147       SAFMAX = ONE / SAFMIN
 148       SSFMAX = SQRT( SAFMAX ) / THREE
 149       SSFMIN = SQRT( SAFMIN ) / EPS2
 150 *
 151 *     Compute the eigenvalues and eigenvectors of the tridiagonal
 152 *     matrix.
 153 *
 154       IF( ICOMPZ.EQ.2 )
 155      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
 156 *
 157       NMAXIT = N*MAXIT
 158       JTOT = 0
 159 *
 160 *     Determine where the matrix splits and choose QL or QR iteration
 161 *     for each block, according to whether top or bottom diagonal
 162 *     element is smaller.
 163 *
 164       L1 = 1
 165       NM1 = N - 1
 166 *
 167    10 CONTINUE
 168       IF( L1.GT.N )
 169      $   GO TO 160
 170       IF( L1.GT.1 )
 171      $   E( L1-1 ) = ZERO
 172       IF( L1.LE.NM1 ) THEN
 173          DO 20 M = L1, NM1
 174             TST = ABS( E( M ) )
 175             IF( TST.EQ.ZERO )
 176      $         GO TO 30
 177             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
 178      $          1 ) ) ) )*EPS ) THEN
 179                E( M ) = ZERO
 180                GO TO 30
 181             END IF
 182    20    CONTINUE
 183       END IF
 184       M = N
 185 *
 186    30 CONTINUE
 187       L = L1
 188       LSV = L
 189       LEND = M
 190       LENDSV = LEND
 191       L1 = M + 1
 192       IF( LEND.EQ.L )
 193      $   GO TO 10
 194 *
 195 *     Scale submatrix in rows and columns L to LEND
 196 *
 197       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
 198       ISCALE = 0
 199       IF( ANORM.EQ.ZERO )
 200      $   GO TO 10
 201       IF( ANORM.GT.SSFMAX ) THEN
 202          ISCALE = 1
 203          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
 204      $                INFO )
 205          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
 206      $                INFO )
 207       ELSE IF( ANORM.LT.SSFMIN ) THEN
 208          ISCALE = 2
 209          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
 210      $                INFO )
 211          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
 212      $                INFO )
 213       END IF
 214 *
 215 *     Choose between QL and QR iteration
 216 *
 217       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
 218          LEND = LSV
 219          L = LENDSV
 220       END IF
 221 *
 222       IF( LEND.GT.L ) THEN
 223 *
 224 *        QL Iteration
 225 *
 226 *        Look for small subdiagonal element.
 227 *
 228    40    CONTINUE
 229          IF( L.NE.LEND ) THEN
 230             LENDM1 = LEND - 1
 231             DO 50 M = L, LENDM1
 232                TST = ABS( E( M ) )**2
 233                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
 234      $             SAFMIN )GO TO 60
 235    50       CONTINUE
 236          END IF
 237 *
 238          M = LEND
 239 *
 240    60    CONTINUE
 241          IF( M.LT.LEND )
 242      $      E( M ) = ZERO
 243          P = D( L )
 244          IF( M.EQ.L )
 245      $      GO TO 80
 246 *
 247 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 248 *        to compute its eigensystem.
 249 *
 250          IF( M.EQ.L+1 ) THEN
 251             IF( ICOMPZ.GT.0 ) THEN
 252                CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
 253                WORK( L ) = C
 254                WORK( N-1+L ) = S
 255                CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ),
 256      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
 257             ELSE
 258                CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
 259             END IF
 260             D( L ) = RT1
 261             D( L+1 ) = RT2
 262             E( L ) = ZERO
 263             L = L + 2
 264             IF( L.LE.LEND )
 265      $         GO TO 40
 266             GO TO 140
 267          END IF
 268 *
 269          IF( JTOT.EQ.NMAXIT )
 270      $      GO TO 140
 271          JTOT = JTOT + 1
 272 *
 273 *        Form shift.
 274 *
 275          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
 276          R = DLAPY2( G, ONE )
 277          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
 278 *
 279          S = ONE
 280          C = ONE
 281          P = ZERO
 282 *
 283 *        Inner loop
 284 *
 285          MM1 = M - 1
 286          DO 70 I = MM1, L, -1
 287             F = S*E( I )
 288             B = C*E( I )
 289             CALL DLARTG( G, F, C, S, R )
 290             IF( I.NE.M-1 )
 291      $         E( I+1 ) = R
 292             G = D( I+1 ) - P
 293             R = ( D( I )-G )*S + TWO*C*B
 294             P = S*R
 295             D( I+1 ) = G + P
 296             G = C*R - B
 297 *
 298 *           If eigenvectors are desired, then save rotations.
 299 *
 300             IF( ICOMPZ.GT.0 ) THEN
 301                WORK( I ) = C
 302                WORK( N-1+I ) = -S
 303             END IF
 304 *
 305    70    CONTINUE
 306 *
 307 *        If eigenvectors are desired, then apply saved rotations.
 308 *
 309          IF( ICOMPZ.GT.0 ) THEN
 310             MM = M - L + 1
 311             CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
 312      $                  Z( 1, L ), LDZ )
 313          END IF
 314 *
 315          D( L ) = D( L ) - P
 316          E( L ) = G
 317          GO TO 40
 318 *
 319 *        Eigenvalue found.
 320 *
 321    80    CONTINUE
 322          D( L ) = P
 323 *
 324          L = L + 1
 325          IF( L.LE.LEND )
 326      $      GO TO 40
 327          GO TO 140
 328 *
 329       ELSE
 330 *
 331 *        QR Iteration
 332 *
 333 *        Look for small superdiagonal element.
 334 *
 335    90    CONTINUE
 336          IF( L.NE.LEND ) THEN
 337             LENDP1 = LEND + 1
 338             DO 100 M = L, LENDP1, -1
 339                TST = ABS( E( M-1 ) )**2
 340                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
 341      $             SAFMIN )GO TO 110
 342   100       CONTINUE
 343          END IF
 344 *
 345          M = LEND
 346 *
 347   110    CONTINUE
 348          IF( M.GT.LEND )
 349      $      E( M-1 ) = ZERO
 350          P = D( L )
 351          IF( M.EQ.L )
 352      $      GO TO 130
 353 *
 354 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 355 *        to compute its eigensystem.
 356 *
 357          IF( M.EQ.L-1 ) THEN
 358             IF( ICOMPZ.GT.0 ) THEN
 359                CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
 360                WORK( M ) = C
 361                WORK( N-1+M ) = S
 362                CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ),
 363      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
 364             ELSE
 365                CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
 366             END IF
 367             D( L-1 ) = RT1
 368             D( L ) = RT2
 369             E( L-1 ) = ZERO
 370             L = L - 2
 371             IF( L.GE.LEND )
 372      $         GO TO 90
 373             GO TO 140
 374          END IF
 375 *
 376          IF( JTOT.EQ.NMAXIT )
 377      $      GO TO 140
 378          JTOT = JTOT + 1
 379 *
 380 *        Form shift.
 381 *
 382          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
 383          R = DLAPY2( G, ONE )
 384          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
 385 *
 386          S = ONE
 387          C = ONE
 388          P = ZERO
 389 *
 390 *        Inner loop
 391 *
 392          LM1 = L - 1
 393          DO 120 I = M, LM1
 394             F = S*E( I )
 395             B = C*E( I )
 396             CALL DLARTG( G, F, C, S, R )
 397             IF( I.NE.M )
 398      $         E( I-1 ) = R
 399             G = D( I ) - P
 400             R = ( D( I+1 )-G )*S + TWO*C*B
 401             P = S*R
 402             D( I ) = G + P
 403             G = C*R - B
 404 *
 405 *           If eigenvectors are desired, then save rotations.
 406 *
 407             IF( ICOMPZ.GT.0 ) THEN
 408                WORK( I ) = C
 409                WORK( N-1+I ) = S
 410             END IF
 411 *
 412   120    CONTINUE
 413 *
 414 *        If eigenvectors are desired, then apply saved rotations.
 415 *
 416          IF( ICOMPZ.GT.0 ) THEN
 417             MM = L - M + 1
 418             CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
 419      $                  Z( 1, M ), LDZ )
 420          END IF
 421 *
 422          D( L ) = D( L ) - P
 423          E( LM1 ) = G
 424          GO TO 90
 425 *
 426 *        Eigenvalue found.
 427 *
 428   130    CONTINUE
 429          D( L ) = P
 430 *
 431          L = L - 1
 432          IF( L.GE.LEND )
 433      $      GO TO 90
 434          GO TO 140
 435 *
 436       END IF
 437 *
 438 *     Undo scaling if necessary
 439 *
 440   140 CONTINUE
 441       IF( ISCALE.EQ.1 ) THEN
 442          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
 443      $                D( LSV ), N, INFO )
 444          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
 445      $                N, INFO )
 446       ELSE IF( ISCALE.EQ.2 ) THEN
 447          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
 448      $                D( LSV ), N, INFO )
 449          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
 450      $                N, INFO )
 451       END IF
 452 *
 453 *     Check for no convergence to an eigenvalue after a total
 454 *     of N*MAXIT iterations.
 455 *
 456       IF( JTOT.LT.NMAXIT )
 457      $   GO TO 10
 458       DO 150 I = 1, N - 1
 459          IF( E( I ).NE.ZERO )
 460      $      INFO = INFO + 1
 461   150 CONTINUE
 462       GO TO 190
 463 *
 464 *     Order eigenvalues and eigenvectors.
 465 *
 466   160 CONTINUE
 467       IF( ICOMPZ.EQ.0 ) THEN
 468 *
 469 *        Use Quick Sort
 470 *
 471          CALL DLASRT( 'I', N, D, INFO )
 472 *
 473       ELSE
 474 *
 475 *        Use Selection Sort to minimize swaps of eigenvectors
 476 *
 477          DO 180 II = 2, N
 478             I = II - 1
 479             K = I
 480             P = D( I )
 481             DO 170 J = II, N
 482                IF( D( J ).LT.P ) THEN
 483                   K = J
 484                   P = D( J )
 485                END IF
 486   170       CONTINUE
 487             IF( K.NE.I ) THEN
 488                D( K ) = D( I )
 489                D( I ) = P
 490                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
 491             END IF
 492   180    CONTINUE
 493       END IF
 494 *
 495   190 CONTINUE
 496       RETURN
 497 *
 498 *     End of DSTEQR
 499 *
 500       END