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