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