![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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