![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 2: $ LDVR, MM, M, WORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER HOWMNY, SIDE 11: INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 12: * .. 13: * .. Array Arguments .. 14: LOGICAL SELECT( * ) 15: DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 16: $ WORK( * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * DTREVC computes some or all of the right and/or left eigenvectors of 23: * a real upper quasi-triangular matrix T. 24: * Matrices of this type are produced by the Schur factorization of 25: * a real general matrix: A = Q*T*Q**T, as computed by DHSEQR. 26: * 27: * The right eigenvector x and the left eigenvector y of T corresponding 28: * to an eigenvalue w are defined by: 29: * 30: * T*x = w*x, (y**H)*T = w*(y**H) 31: * 32: * where y**H denotes the conjugate transpose of y. 33: * The eigenvalues are not input to this routine, but are read directly 34: * from the diagonal blocks of T. 35: * 36: * This routine returns the matrices X and/or Y of right and left 37: * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 38: * input matrix. If Q is the orthogonal factor that reduces a matrix 39: * A to Schur form T, then Q*X and Q*Y are the matrices of right and 40: * left eigenvectors of A. 41: * 42: * Arguments 43: * ========= 44: * 45: * SIDE (input) CHARACTER*1 46: * = 'R': compute right eigenvectors only; 47: * = 'L': compute left eigenvectors only; 48: * = 'B': compute both right and left eigenvectors. 49: * 50: * HOWMNY (input) CHARACTER*1 51: * = 'A': compute all right and/or left eigenvectors; 52: * = 'B': compute all right and/or left eigenvectors, 53: * backtransformed by the matrices in VR and/or VL; 54: * = 'S': compute selected right and/or left eigenvectors, 55: * as indicated by the logical array SELECT. 56: * 57: * SELECT (input/output) LOGICAL array, dimension (N) 58: * If HOWMNY = 'S', SELECT specifies the eigenvectors to be 59: * computed. 60: * If w(j) is a real eigenvalue, the corresponding real 61: * eigenvector is computed if SELECT(j) is .TRUE.. 62: * If w(j) and w(j+1) are the real and imaginary parts of a 63: * complex eigenvalue, the corresponding complex eigenvector is 64: * computed if either SELECT(j) or SELECT(j+1) is .TRUE., and 65: * on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to 66: * .FALSE.. 67: * Not referenced if HOWMNY = 'A' or 'B'. 68: * 69: * N (input) INTEGER 70: * The order of the matrix T. N >= 0. 71: * 72: * T (input) DOUBLE PRECISION array, dimension (LDT,N) 73: * The upper quasi-triangular matrix T in Schur canonical form. 74: * 75: * LDT (input) INTEGER 76: * The leading dimension of the array T. LDT >= max(1,N). 77: * 78: * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) 79: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 80: * contain an N-by-N matrix Q (usually the orthogonal matrix Q 81: * of Schur vectors returned by DHSEQR). 82: * On exit, if SIDE = 'L' or 'B', VL contains: 83: * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 84: * if HOWMNY = 'B', the matrix Q*Y; 85: * if HOWMNY = 'S', the left eigenvectors of T specified by 86: * SELECT, stored consecutively in the columns 87: * of VL, in the same order as their 88: * eigenvalues. 89: * A complex eigenvector corresponding to a complex eigenvalue 90: * is stored in two consecutive columns, the first holding the 91: * real part, and the second the imaginary part. 92: * Not referenced if SIDE = 'R'. 93: * 94: * LDVL (input) INTEGER 95: * The leading dimension of the array VL. LDVL >= 1, and if 96: * SIDE = 'L' or 'B', LDVL >= N. 97: * 98: * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) 99: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 100: * contain an N-by-N matrix Q (usually the orthogonal matrix Q 101: * of Schur vectors returned by DHSEQR). 102: * On exit, if SIDE = 'R' or 'B', VR contains: 103: * if HOWMNY = 'A', the matrix X of right eigenvectors of T; 104: * if HOWMNY = 'B', the matrix Q*X; 105: * if HOWMNY = 'S', the right eigenvectors of T specified by 106: * SELECT, stored consecutively in the columns 107: * of VR, in the same order as their 108: * eigenvalues. 109: * A complex eigenvector corresponding to a complex eigenvalue 110: * is stored in two consecutive columns, the first holding the 111: * real part and the second the imaginary part. 112: * Not referenced if SIDE = 'L'. 113: * 114: * LDVR (input) INTEGER 115: * The leading dimension of the array VR. LDVR >= 1, and if 116: * SIDE = 'R' or 'B', LDVR >= N. 117: * 118: * MM (input) INTEGER 119: * The number of columns in the arrays VL and/or VR. MM >= M. 120: * 121: * M (output) INTEGER 122: * The number of columns in the arrays VL and/or VR actually 123: * used to store the eigenvectors. 124: * If HOWMNY = 'A' or 'B', M is set to N. 125: * Each selected real eigenvector occupies one column and each 126: * selected complex eigenvector occupies two columns. 127: * 128: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 129: * 130: * INFO (output) INTEGER 131: * = 0: successful exit 132: * < 0: if INFO = -i, the i-th argument had an illegal value 133: * 134: * Further Details 135: * =============== 136: * 137: * The algorithm used in this program is basically backward (forward) 138: * substitution, with scaling to make the the code robust against 139: * possible overflow. 140: * 141: * Each eigenvector is normalized so that the element of largest 142: * magnitude has magnitude 1; here the magnitude of a complex number 143: * (x,y) is taken to be |x| + |y|. 144: * 145: * ===================================================================== 146: * 147: * .. Parameters .. 148: DOUBLE PRECISION ZERO, ONE 149: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 150: * .. 151: * .. Local Scalars .. 152: LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV 153: INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2 154: DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE, 155: $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR, 156: $ XNORM 157: * .. 158: * .. External Functions .. 159: LOGICAL LSAME 160: INTEGER IDAMAX 161: DOUBLE PRECISION DDOT, DLAMCH 162: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 163: * .. 164: * .. External Subroutines .. 165: EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA 166: * .. 167: * .. Intrinsic Functions .. 168: INTRINSIC ABS, MAX, SQRT 169: * .. 170: * .. Local Arrays .. 171: DOUBLE PRECISION X( 2, 2 ) 172: * .. 173: * .. Executable Statements .. 174: * 175: * Decode and test the input parameters 176: * 177: BOTHV = LSAME( SIDE, 'B' ) 178: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 179: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 180: * 181: ALLV = LSAME( HOWMNY, 'A' ) 182: OVER = LSAME( HOWMNY, 'B' ) 183: SOMEV = LSAME( HOWMNY, 'S' ) 184: * 185: INFO = 0 186: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 187: INFO = -1 188: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 189: INFO = -2 190: ELSE IF( N.LT.0 ) THEN 191: INFO = -4 192: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 193: INFO = -6 194: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 195: INFO = -8 196: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 197: INFO = -10 198: ELSE 199: * 200: * Set M to the number of columns required to store the selected 201: * eigenvectors, standardize the array SELECT if necessary, and 202: * test MM. 203: * 204: IF( SOMEV ) THEN 205: M = 0 206: PAIR = .FALSE. 207: DO 10 J = 1, N 208: IF( PAIR ) THEN 209: PAIR = .FALSE. 210: SELECT( J ) = .FALSE. 211: ELSE 212: IF( J.LT.N ) THEN 213: IF( T( J+1, J ).EQ.ZERO ) THEN 214: IF( SELECT( J ) ) 215: $ M = M + 1 216: ELSE 217: PAIR = .TRUE. 218: IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN 219: SELECT( J ) = .TRUE. 220: M = M + 2 221: END IF 222: END IF 223: ELSE 224: IF( SELECT( N ) ) 225: $ M = M + 1 226: END IF 227: END IF 228: 10 CONTINUE 229: ELSE 230: M = N 231: END IF 232: * 233: IF( MM.LT.M ) THEN 234: INFO = -11 235: END IF 236: END IF 237: IF( INFO.NE.0 ) THEN 238: CALL XERBLA( 'DTREVC', -INFO ) 239: RETURN 240: END IF 241: * 242: * Quick return if possible. 243: * 244: IF( N.EQ.0 ) 245: $ RETURN 246: * 247: * Set the constants to control overflow. 248: * 249: UNFL = DLAMCH( 'Safe minimum' ) 250: OVFL = ONE / UNFL 251: CALL DLABAD( UNFL, OVFL ) 252: ULP = DLAMCH( 'Precision' ) 253: SMLNUM = UNFL*( N / ULP ) 254: BIGNUM = ( ONE-ULP ) / SMLNUM 255: * 256: * Compute 1-norm of each column of strictly upper triangular 257: * part of T to control overflow in triangular solver. 258: * 259: WORK( 1 ) = ZERO 260: DO 30 J = 2, N 261: WORK( J ) = ZERO 262: DO 20 I = 1, J - 1 263: WORK( J ) = WORK( J ) + ABS( T( I, J ) ) 264: 20 CONTINUE 265: 30 CONTINUE 266: * 267: * Index IP is used to specify the real or complex eigenvalue: 268: * IP = 0, real eigenvalue, 269: * 1, first of conjugate complex pair: (wr,wi) 270: * -1, second of conjugate complex pair: (wr,wi) 271: * 272: N2 = 2*N 273: * 274: IF( RIGHTV ) THEN 275: * 276: * Compute right eigenvectors. 277: * 278: IP = 0 279: IS = M 280: DO 140 KI = N, 1, -1 281: * 282: IF( IP.EQ.1 ) 283: $ GO TO 130 284: IF( KI.EQ.1 ) 285: $ GO TO 40 286: IF( T( KI, KI-1 ).EQ.ZERO ) 287: $ GO TO 40 288: IP = -1 289: * 290: 40 CONTINUE 291: IF( SOMEV ) THEN 292: IF( IP.EQ.0 ) THEN 293: IF( .NOT.SELECT( KI ) ) 294: $ GO TO 130 295: ELSE 296: IF( .NOT.SELECT( KI-1 ) ) 297: $ GO TO 130 298: END IF 299: END IF 300: * 301: * Compute the KI-th eigenvalue (WR,WI). 302: * 303: WR = T( KI, KI ) 304: WI = ZERO 305: IF( IP.NE.0 ) 306: $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* 307: $ SQRT( ABS( T( KI-1, KI ) ) ) 308: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 309: * 310: IF( IP.EQ.0 ) THEN 311: * 312: * Real right eigenvector 313: * 314: WORK( KI+N ) = ONE 315: * 316: * Form right-hand side 317: * 318: DO 50 K = 1, KI - 1 319: WORK( K+N ) = -T( K, KI ) 320: 50 CONTINUE 321: * 322: * Solve the upper quasi-triangular system: 323: * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. 324: * 325: JNXT = KI - 1 326: DO 60 J = KI - 1, 1, -1 327: IF( J.GT.JNXT ) 328: $ GO TO 60 329: J1 = J 330: J2 = J 331: JNXT = J - 1 332: IF( J.GT.1 ) THEN 333: IF( T( J, J-1 ).NE.ZERO ) THEN 334: J1 = J - 1 335: JNXT = J - 2 336: END IF 337: END IF 338: * 339: IF( J1.EQ.J2 ) THEN 340: * 341: * 1-by-1 diagonal block 342: * 343: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 344: $ LDT, ONE, ONE, WORK( J+N ), N, WR, 345: $ ZERO, X, 2, SCALE, XNORM, IERR ) 346: * 347: * Scale X(1,1) to avoid overflow when updating 348: * the right-hand side. 349: * 350: IF( XNORM.GT.ONE ) THEN 351: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 352: X( 1, 1 ) = X( 1, 1 ) / XNORM 353: SCALE = SCALE / XNORM 354: END IF 355: END IF 356: * 357: * Scale if necessary 358: * 359: IF( SCALE.NE.ONE ) 360: $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 361: WORK( J+N ) = X( 1, 1 ) 362: * 363: * Update right-hand side 364: * 365: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 366: $ WORK( 1+N ), 1 ) 367: * 368: ELSE 369: * 370: * 2-by-2 diagonal block 371: * 372: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, 373: $ T( J-1, J-1 ), LDT, ONE, ONE, 374: $ WORK( J-1+N ), N, WR, ZERO, X, 2, 375: $ SCALE, XNORM, IERR ) 376: * 377: * Scale X(1,1) and X(2,1) to avoid overflow when 378: * updating the right-hand side. 379: * 380: IF( XNORM.GT.ONE ) THEN 381: BETA = MAX( WORK( J-1 ), WORK( J ) ) 382: IF( BETA.GT.BIGNUM / XNORM ) THEN 383: X( 1, 1 ) = X( 1, 1 ) / XNORM 384: X( 2, 1 ) = X( 2, 1 ) / XNORM 385: SCALE = SCALE / XNORM 386: END IF 387: END IF 388: * 389: * Scale if necessary 390: * 391: IF( SCALE.NE.ONE ) 392: $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 393: WORK( J-1+N ) = X( 1, 1 ) 394: WORK( J+N ) = X( 2, 1 ) 395: * 396: * Update right-hand side 397: * 398: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 399: $ WORK( 1+N ), 1 ) 400: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 401: $ WORK( 1+N ), 1 ) 402: END IF 403: 60 CONTINUE 404: * 405: * Copy the vector x or Q*x to VR and normalize. 406: * 407: IF( .NOT.OVER ) THEN 408: CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) 409: * 410: II = IDAMAX( KI, VR( 1, IS ), 1 ) 411: REMAX = ONE / ABS( VR( II, IS ) ) 412: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 413: * 414: DO 70 K = KI + 1, N 415: VR( K, IS ) = ZERO 416: 70 CONTINUE 417: ELSE 418: IF( KI.GT.1 ) 419: $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR, 420: $ WORK( 1+N ), 1, WORK( KI+N ), 421: $ VR( 1, KI ), 1 ) 422: * 423: II = IDAMAX( N, VR( 1, KI ), 1 ) 424: REMAX = ONE / ABS( VR( II, KI ) ) 425: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 426: END IF 427: * 428: ELSE 429: * 430: * Complex right eigenvector. 431: * 432: * Initial solve 433: * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. 434: * [ (T(KI,KI-1) T(KI,KI) ) ] 435: * 436: IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN 437: WORK( KI-1+N ) = ONE 438: WORK( KI+N2 ) = WI / T( KI-1, KI ) 439: ELSE 440: WORK( KI-1+N ) = -WI / T( KI, KI-1 ) 441: WORK( KI+N2 ) = ONE 442: END IF 443: WORK( KI+N ) = ZERO 444: WORK( KI-1+N2 ) = ZERO 445: * 446: * Form right-hand side 447: * 448: DO 80 K = 1, KI - 2 449: WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) 450: WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) 451: 80 CONTINUE 452: * 453: * Solve upper quasi-triangular system: 454: * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) 455: * 456: JNXT = KI - 2 457: DO 90 J = KI - 2, 1, -1 458: IF( J.GT.JNXT ) 459: $ GO TO 90 460: J1 = J 461: J2 = J 462: JNXT = J - 1 463: IF( J.GT.1 ) THEN 464: IF( T( J, J-1 ).NE.ZERO ) THEN 465: J1 = J - 1 466: JNXT = J - 2 467: END IF 468: END IF 469: * 470: IF( J1.EQ.J2 ) THEN 471: * 472: * 1-by-1 diagonal block 473: * 474: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 475: $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, 476: $ X, 2, SCALE, XNORM, IERR ) 477: * 478: * Scale X(1,1) and X(1,2) to avoid overflow when 479: * updating the right-hand side. 480: * 481: IF( XNORM.GT.ONE ) THEN 482: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 483: X( 1, 1 ) = X( 1, 1 ) / XNORM 484: X( 1, 2 ) = X( 1, 2 ) / XNORM 485: SCALE = SCALE / XNORM 486: END IF 487: END IF 488: * 489: * Scale if necessary 490: * 491: IF( SCALE.NE.ONE ) THEN 492: CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 493: CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 494: END IF 495: WORK( J+N ) = X( 1, 1 ) 496: WORK( J+N2 ) = X( 1, 2 ) 497: * 498: * Update the right-hand side 499: * 500: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 501: $ WORK( 1+N ), 1 ) 502: CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, 503: $ WORK( 1+N2 ), 1 ) 504: * 505: ELSE 506: * 507: * 2-by-2 diagonal block 508: * 509: CALL DLALN2( .FALSE., 2, 2, SMIN, ONE, 510: $ T( J-1, J-1 ), LDT, ONE, ONE, 511: $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, 512: $ XNORM, IERR ) 513: * 514: * Scale X to avoid overflow when updating 515: * the right-hand side. 516: * 517: IF( XNORM.GT.ONE ) THEN 518: BETA = MAX( WORK( J-1 ), WORK( J ) ) 519: IF( BETA.GT.BIGNUM / XNORM ) THEN 520: REC = ONE / XNORM 521: X( 1, 1 ) = X( 1, 1 )*REC 522: X( 1, 2 ) = X( 1, 2 )*REC 523: X( 2, 1 ) = X( 2, 1 )*REC 524: X( 2, 2 ) = X( 2, 2 )*REC 525: SCALE = SCALE*REC 526: END IF 527: END IF 528: * 529: * Scale if necessary 530: * 531: IF( SCALE.NE.ONE ) THEN 532: CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 533: CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 534: END IF 535: WORK( J-1+N ) = X( 1, 1 ) 536: WORK( J+N ) = X( 2, 1 ) 537: WORK( J-1+N2 ) = X( 1, 2 ) 538: WORK( J+N2 ) = X( 2, 2 ) 539: * 540: * Update the right-hand side 541: * 542: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 543: $ WORK( 1+N ), 1 ) 544: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 545: $ WORK( 1+N ), 1 ) 546: CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, 547: $ WORK( 1+N2 ), 1 ) 548: CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, 549: $ WORK( 1+N2 ), 1 ) 550: END IF 551: 90 CONTINUE 552: * 553: * Copy the vector x or Q*x to VR and normalize. 554: * 555: IF( .NOT.OVER ) THEN 556: CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) 557: CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) 558: * 559: EMAX = ZERO 560: DO 100 K = 1, KI 561: EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ 562: $ ABS( VR( K, IS ) ) ) 563: 100 CONTINUE 564: * 565: REMAX = ONE / EMAX 566: CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) 567: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 568: * 569: DO 110 K = KI + 1, N 570: VR( K, IS-1 ) = ZERO 571: VR( K, IS ) = ZERO 572: 110 CONTINUE 573: * 574: ELSE 575: * 576: IF( KI.GT.2 ) THEN 577: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 578: $ WORK( 1+N ), 1, WORK( KI-1+N ), 579: $ VR( 1, KI-1 ), 1 ) 580: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 581: $ WORK( 1+N2 ), 1, WORK( KI+N2 ), 582: $ VR( 1, KI ), 1 ) 583: ELSE 584: CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) 585: CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) 586: END IF 587: * 588: EMAX = ZERO 589: DO 120 K = 1, N 590: EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ 591: $ ABS( VR( K, KI ) ) ) 592: 120 CONTINUE 593: REMAX = ONE / EMAX 594: CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) 595: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 596: END IF 597: END IF 598: * 599: IS = IS - 1 600: IF( IP.NE.0 ) 601: $ IS = IS - 1 602: 130 CONTINUE 603: IF( IP.EQ.1 ) 604: $ IP = 0 605: IF( IP.EQ.-1 ) 606: $ IP = 1 607: 140 CONTINUE 608: END IF 609: * 610: IF( LEFTV ) THEN 611: * 612: * Compute left eigenvectors. 613: * 614: IP = 0 615: IS = 1 616: DO 260 KI = 1, N 617: * 618: IF( IP.EQ.-1 ) 619: $ GO TO 250 620: IF( KI.EQ.N ) 621: $ GO TO 150 622: IF( T( KI+1, KI ).EQ.ZERO ) 623: $ GO TO 150 624: IP = 1 625: * 626: 150 CONTINUE 627: IF( SOMEV ) THEN 628: IF( .NOT.SELECT( KI ) ) 629: $ GO TO 250 630: END IF 631: * 632: * Compute the KI-th eigenvalue (WR,WI). 633: * 634: WR = T( KI, KI ) 635: WI = ZERO 636: IF( IP.NE.0 ) 637: $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* 638: $ SQRT( ABS( T( KI+1, KI ) ) ) 639: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 640: * 641: IF( IP.EQ.0 ) THEN 642: * 643: * Real left eigenvector. 644: * 645: WORK( KI+N ) = ONE 646: * 647: * Form right-hand side 648: * 649: DO 160 K = KI + 1, N 650: WORK( K+N ) = -T( KI, K ) 651: 160 CONTINUE 652: * 653: * Solve the quasi-triangular system: 654: * (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK 655: * 656: VMAX = ONE 657: VCRIT = BIGNUM 658: * 659: JNXT = KI + 1 660: DO 170 J = KI + 1, N 661: IF( J.LT.JNXT ) 662: $ GO TO 170 663: J1 = J 664: J2 = J 665: JNXT = J + 1 666: IF( J.LT.N ) THEN 667: IF( T( J+1, J ).NE.ZERO ) THEN 668: J2 = J + 1 669: JNXT = J + 2 670: END IF 671: END IF 672: * 673: IF( J1.EQ.J2 ) THEN 674: * 675: * 1-by-1 diagonal block 676: * 677: * Scale if necessary to avoid overflow when forming 678: * the right-hand side. 679: * 680: IF( WORK( J ).GT.VCRIT ) THEN 681: REC = ONE / VMAX 682: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 683: VMAX = ONE 684: VCRIT = BIGNUM 685: END IF 686: * 687: WORK( J+N ) = WORK( J+N ) - 688: $ DDOT( J-KI-1, T( KI+1, J ), 1, 689: $ WORK( KI+1+N ), 1 ) 690: * 691: * Solve (T(J,J)-WR)'*X = WORK 692: * 693: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 694: $ LDT, ONE, ONE, WORK( J+N ), N, WR, 695: $ ZERO, X, 2, SCALE, XNORM, IERR ) 696: * 697: * Scale if necessary 698: * 699: IF( SCALE.NE.ONE ) 700: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 701: WORK( J+N ) = X( 1, 1 ) 702: VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) 703: VCRIT = BIGNUM / VMAX 704: * 705: ELSE 706: * 707: * 2-by-2 diagonal block 708: * 709: * Scale if necessary to avoid overflow when forming 710: * the right-hand side. 711: * 712: BETA = MAX( WORK( J ), WORK( J+1 ) ) 713: IF( BETA.GT.VCRIT ) THEN 714: REC = ONE / VMAX 715: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 716: VMAX = ONE 717: VCRIT = BIGNUM 718: END IF 719: * 720: WORK( J+N ) = WORK( J+N ) - 721: $ DDOT( J-KI-1, T( KI+1, J ), 1, 722: $ WORK( KI+1+N ), 1 ) 723: * 724: WORK( J+1+N ) = WORK( J+1+N ) - 725: $ DDOT( J-KI-1, T( KI+1, J+1 ), 1, 726: $ WORK( KI+1+N ), 1 ) 727: * 728: * Solve 729: * [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 ) 730: * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) 731: * 732: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), 733: $ LDT, ONE, ONE, WORK( J+N ), N, WR, 734: $ ZERO, X, 2, SCALE, XNORM, IERR ) 735: * 736: * Scale if necessary 737: * 738: IF( SCALE.NE.ONE ) 739: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 740: WORK( J+N ) = X( 1, 1 ) 741: WORK( J+1+N ) = X( 2, 1 ) 742: * 743: VMAX = MAX( ABS( WORK( J+N ) ), 744: $ ABS( WORK( J+1+N ) ), VMAX ) 745: VCRIT = BIGNUM / VMAX 746: * 747: END IF 748: 170 CONTINUE 749: * 750: * Copy the vector x or Q*x to VL and normalize. 751: * 752: IF( .NOT.OVER ) THEN 753: CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 754: * 755: II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 756: REMAX = ONE / ABS( VL( II, IS ) ) 757: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 758: * 759: DO 180 K = 1, KI - 1 760: VL( K, IS ) = ZERO 761: 180 CONTINUE 762: * 763: ELSE 764: * 765: IF( KI.LT.N ) 766: $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, 767: $ WORK( KI+1+N ), 1, WORK( KI+N ), 768: $ VL( 1, KI ), 1 ) 769: * 770: II = IDAMAX( N, VL( 1, KI ), 1 ) 771: REMAX = ONE / ABS( VL( II, KI ) ) 772: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 773: * 774: END IF 775: * 776: ELSE 777: * 778: * Complex left eigenvector. 779: * 780: * Initial solve: 781: * ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0. 782: * ((T(KI+1,KI) T(KI+1,KI+1)) ) 783: * 784: IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN 785: WORK( KI+N ) = WI / T( KI, KI+1 ) 786: WORK( KI+1+N2 ) = ONE 787: ELSE 788: WORK( KI+N ) = ONE 789: WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) 790: END IF 791: WORK( KI+1+N ) = ZERO 792: WORK( KI+N2 ) = ZERO 793: * 794: * Form right-hand side 795: * 796: DO 190 K = KI + 2, N 797: WORK( K+N ) = -WORK( KI+N )*T( KI, K ) 798: WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) 799: 190 CONTINUE 800: * 801: * Solve complex quasi-triangular system: 802: * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 803: * 804: VMAX = ONE 805: VCRIT = BIGNUM 806: * 807: JNXT = KI + 2 808: DO 200 J = KI + 2, N 809: IF( J.LT.JNXT ) 810: $ GO TO 200 811: J1 = J 812: J2 = J 813: JNXT = J + 1 814: IF( J.LT.N ) THEN 815: IF( T( J+1, J ).NE.ZERO ) THEN 816: J2 = J + 1 817: JNXT = J + 2 818: END IF 819: END IF 820: * 821: IF( J1.EQ.J2 ) THEN 822: * 823: * 1-by-1 diagonal block 824: * 825: * Scale if necessary to avoid overflow when 826: * forming the right-hand side elements. 827: * 828: IF( WORK( J ).GT.VCRIT ) THEN 829: REC = ONE / VMAX 830: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 831: CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 832: VMAX = ONE 833: VCRIT = BIGNUM 834: END IF 835: * 836: WORK( J+N ) = WORK( J+N ) - 837: $ DDOT( J-KI-2, T( KI+2, J ), 1, 838: $ WORK( KI+2+N ), 1 ) 839: WORK( J+N2 ) = WORK( J+N2 ) - 840: $ DDOT( J-KI-2, T( KI+2, J ), 1, 841: $ WORK( KI+2+N2 ), 1 ) 842: * 843: * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 844: * 845: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 846: $ LDT, ONE, ONE, WORK( J+N ), N, WR, 847: $ -WI, X, 2, SCALE, XNORM, IERR ) 848: * 849: * Scale if necessary 850: * 851: IF( SCALE.NE.ONE ) THEN 852: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 853: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 854: END IF 855: WORK( J+N ) = X( 1, 1 ) 856: WORK( J+N2 ) = X( 1, 2 ) 857: VMAX = MAX( ABS( WORK( J+N ) ), 858: $ ABS( WORK( J+N2 ) ), VMAX ) 859: VCRIT = BIGNUM / VMAX 860: * 861: ELSE 862: * 863: * 2-by-2 diagonal block 864: * 865: * Scale if necessary to avoid overflow when forming 866: * the right-hand side elements. 867: * 868: BETA = MAX( WORK( J ), WORK( J+1 ) ) 869: IF( BETA.GT.VCRIT ) THEN 870: REC = ONE / VMAX 871: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 872: CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 873: VMAX = ONE 874: VCRIT = BIGNUM 875: END IF 876: * 877: WORK( J+N ) = WORK( J+N ) - 878: $ DDOT( J-KI-2, T( KI+2, J ), 1, 879: $ WORK( KI+2+N ), 1 ) 880: * 881: WORK( J+N2 ) = WORK( J+N2 ) - 882: $ DDOT( J-KI-2, T( KI+2, J ), 1, 883: $ WORK( KI+2+N2 ), 1 ) 884: * 885: WORK( J+1+N ) = WORK( J+1+N ) - 886: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 887: $ WORK( KI+2+N ), 1 ) 888: * 889: WORK( J+1+N2 ) = WORK( J+1+N2 ) - 890: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 891: $ WORK( KI+2+N2 ), 1 ) 892: * 893: * Solve 2-by-2 complex linear equation 894: * ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B 895: * ([T(j+1,j) T(j+1,j+1)] ) 896: * 897: CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), 898: $ LDT, ONE, ONE, WORK( J+N ), N, WR, 899: $ -WI, X, 2, SCALE, XNORM, IERR ) 900: * 901: * Scale if necessary 902: * 903: IF( SCALE.NE.ONE ) THEN 904: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 905: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 906: END IF 907: WORK( J+N ) = X( 1, 1 ) 908: WORK( J+N2 ) = X( 1, 2 ) 909: WORK( J+1+N ) = X( 2, 1 ) 910: WORK( J+1+N2 ) = X( 2, 2 ) 911: VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), 912: $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) 913: VCRIT = BIGNUM / VMAX 914: * 915: END IF 916: 200 CONTINUE 917: * 918: * Copy the vector x or Q*x to VL and normalize. 919: * 920: IF( .NOT.OVER ) THEN 921: CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 922: CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), 923: $ 1 ) 924: * 925: EMAX = ZERO 926: DO 220 K = KI, N 927: EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ 928: $ ABS( VL( K, IS+1 ) ) ) 929: 220 CONTINUE 930: REMAX = ONE / EMAX 931: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 932: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) 933: * 934: DO 230 K = 1, KI - 1 935: VL( K, IS ) = ZERO 936: VL( K, IS+1 ) = ZERO 937: 230 CONTINUE 938: ELSE 939: IF( KI.LT.N-1 ) THEN 940: CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 941: $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), 942: $ VL( 1, KI ), 1 ) 943: CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 944: $ LDVL, WORK( KI+2+N2 ), 1, 945: $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 946: ELSE 947: CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) 948: CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 949: END IF 950: * 951: EMAX = ZERO 952: DO 240 K = 1, N 953: EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ 954: $ ABS( VL( K, KI+1 ) ) ) 955: 240 CONTINUE 956: REMAX = ONE / EMAX 957: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 958: CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) 959: * 960: END IF 961: * 962: END IF 963: * 964: IS = IS + 1 965: IF( IP.NE.0 ) 966: $ IS = IS + 1 967: 250 CONTINUE 968: IF( IP.EQ.-1 ) 969: $ IP = 0 970: IF( IP.EQ.1 ) 971: $ IP = -1 972: * 973: 260 CONTINUE 974: * 975: END IF 976: * 977: RETURN 978: * 979: * End of DTREVC 980: * 981: END