Annotation of rpl/lapack/lapack/dtrevc.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrevc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack