src/LAPACK/zsteqr.f

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