Annotation of rpl/lapack/lapack/dtrevc3.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DTREVC3
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTREVC3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
! 22: * VR, LDVR, MM, M, WORK, LWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER HOWMNY, SIDE
! 26: * INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * LOGICAL SELECT( * )
! 30: * DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
! 31: * $ WORK( * )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> DTREVC3 computes some or all of the right and/or left eigenvectors of
! 41: *> a real upper quasi-triangular matrix T.
! 42: *> Matrices of this type are produced by the Schur factorization of
! 43: *> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
! 44: *>
! 45: *> The right eigenvector x and the left eigenvector y of T corresponding
! 46: *> to an eigenvalue w are defined by:
! 47: *>
! 48: *> T*x = w*x, (y**T)*T = w*(y**T)
! 49: *>
! 50: *> where y**T denotes the transpose of the vector y.
! 51: *> The eigenvalues are not input to this routine, but are read directly
! 52: *> from the diagonal blocks of T.
! 53: *>
! 54: *> This routine returns the matrices X and/or Y of right and left
! 55: *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
! 56: *> input matrix. If Q is the orthogonal factor that reduces a matrix
! 57: *> A to Schur form T, then Q*X and Q*Y are the matrices of right and
! 58: *> left eigenvectors of A.
! 59: *>
! 60: *> This uses a Level 3 BLAS version of the back transformation.
! 61: *> \endverbatim
! 62: *
! 63: * Arguments:
! 64: * ==========
! 65: *
! 66: *> \param[in] SIDE
! 67: *> \verbatim
! 68: *> SIDE is CHARACTER*1
! 69: *> = 'R': compute right eigenvectors only;
! 70: *> = 'L': compute left eigenvectors only;
! 71: *> = 'B': compute both right and left eigenvectors.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] HOWMNY
! 75: *> \verbatim
! 76: *> HOWMNY is CHARACTER*1
! 77: *> = 'A': compute all right and/or left eigenvectors;
! 78: *> = 'B': compute all right and/or left eigenvectors,
! 79: *> backtransformed by the matrices in VR and/or VL;
! 80: *> = 'S': compute selected right and/or left eigenvectors,
! 81: *> as indicated by the logical array SELECT.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in,out] SELECT
! 85: *> \verbatim
! 86: *> SELECT is LOGICAL array, dimension (N)
! 87: *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
! 88: *> computed.
! 89: *> If w(j) is a real eigenvalue, the corresponding real
! 90: *> eigenvector is computed if SELECT(j) is .TRUE..
! 91: *> If w(j) and w(j+1) are the real and imaginary parts of a
! 92: *> complex eigenvalue, the corresponding complex eigenvector is
! 93: *> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
! 94: *> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
! 95: *> .FALSE..
! 96: *> Not referenced if HOWMNY = 'A' or 'B'.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] N
! 100: *> \verbatim
! 101: *> N is INTEGER
! 102: *> The order of the matrix T. N >= 0.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] T
! 106: *> \verbatim
! 107: *> T is DOUBLE PRECISION array, dimension (LDT,N)
! 108: *> The upper quasi-triangular matrix T in Schur canonical form.
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in] LDT
! 112: *> \verbatim
! 113: *> LDT is INTEGER
! 114: *> The leading dimension of the array T. LDT >= max(1,N).
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in,out] VL
! 118: *> \verbatim
! 119: *> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
! 120: *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
! 121: *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
! 122: *> of Schur vectors returned by DHSEQR).
! 123: *> On exit, if SIDE = 'L' or 'B', VL contains:
! 124: *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
! 125: *> if HOWMNY = 'B', the matrix Q*Y;
! 126: *> if HOWMNY = 'S', the left eigenvectors of T specified by
! 127: *> SELECT, stored consecutively in the columns
! 128: *> of VL, in the same order as their
! 129: *> eigenvalues.
! 130: *> A complex eigenvector corresponding to a complex eigenvalue
! 131: *> is stored in two consecutive columns, the first holding the
! 132: *> real part, and the second the imaginary part.
! 133: *> Not referenced if SIDE = 'R'.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in] LDVL
! 137: *> \verbatim
! 138: *> LDVL is INTEGER
! 139: *> The leading dimension of the array VL.
! 140: *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
! 141: *> \endverbatim
! 142: *>
! 143: *> \param[in,out] VR
! 144: *> \verbatim
! 145: *> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
! 146: *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
! 147: *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
! 148: *> of Schur vectors returned by DHSEQR).
! 149: *> On exit, if SIDE = 'R' or 'B', VR contains:
! 150: *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
! 151: *> if HOWMNY = 'B', the matrix Q*X;
! 152: *> if HOWMNY = 'S', the right eigenvectors of T specified by
! 153: *> SELECT, stored consecutively in the columns
! 154: *> of VR, in the same order as their
! 155: *> eigenvalues.
! 156: *> A complex eigenvector corresponding to a complex eigenvalue
! 157: *> is stored in two consecutive columns, the first holding the
! 158: *> real part and the second the imaginary part.
! 159: *> Not referenced if SIDE = 'L'.
! 160: *> \endverbatim
! 161: *>
! 162: *> \param[in] LDVR
! 163: *> \verbatim
! 164: *> LDVR is INTEGER
! 165: *> The leading dimension of the array VR.
! 166: *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
! 167: *> \endverbatim
! 168: *>
! 169: *> \param[in] MM
! 170: *> \verbatim
! 171: *> MM is INTEGER
! 172: *> The number of columns in the arrays VL and/or VR. MM >= M.
! 173: *> \endverbatim
! 174: *>
! 175: *> \param[out] M
! 176: *> \verbatim
! 177: *> M is INTEGER
! 178: *> The number of columns in the arrays VL and/or VR actually
! 179: *> used to store the eigenvectors.
! 180: *> If HOWMNY = 'A' or 'B', M is set to N.
! 181: *> Each selected real eigenvector occupies one column and each
! 182: *> selected complex eigenvector occupies two columns.
! 183: *> \endverbatim
! 184: *>
! 185: *> \param[out] WORK
! 186: *> \verbatim
! 187: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 188: *> \endverbatim
! 189: *>
! 190: *> \param[in] LWORK
! 191: *> \verbatim
! 192: *> LWORK is INTEGER
! 193: *> The dimension of array WORK. LWORK >= max(1,3*N).
! 194: *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
! 195: *> the optimal blocksize.
! 196: *>
! 197: *> If LWORK = -1, then a workspace query is assumed; the routine
! 198: *> only calculates the optimal size of the WORK array, returns
! 199: *> this value as the first entry of the WORK array, and no error
! 200: *> message related to LWORK is issued by XERBLA.
! 201: *> \endverbatim
! 202: *>
! 203: *> \param[out] INFO
! 204: *> \verbatim
! 205: *> INFO is INTEGER
! 206: *> = 0: successful exit
! 207: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 208: *> \endverbatim
! 209: *
! 210: * Authors:
! 211: * ========
! 212: *
! 213: *> \author Univ. of Tennessee
! 214: *> \author Univ. of California Berkeley
! 215: *> \author Univ. of Colorado Denver
! 216: *> \author NAG Ltd.
! 217: *
! 218: *> \date November 2011
! 219: *
! 220: * @precisions fortran d -> s
! 221: *
! 222: *> \ingroup doubleOTHERcomputational
! 223: *
! 224: *> \par Further Details:
! 225: * =====================
! 226: *>
! 227: *> \verbatim
! 228: *>
! 229: *> The algorithm used in this program is basically backward (forward)
! 230: *> substitution, with scaling to make the the code robust against
! 231: *> possible overflow.
! 232: *>
! 233: *> Each eigenvector is normalized so that the element of largest
! 234: *> magnitude has magnitude 1; here the magnitude of a complex number
! 235: *> (x,y) is taken to be |x| + |y|.
! 236: *> \endverbatim
! 237: *>
! 238: * =====================================================================
! 239: SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
! 240: $ VR, LDVR, MM, M, WORK, LWORK, INFO )
! 241: IMPLICIT NONE
! 242: *
! 243: * -- LAPACK computational routine (version 3.4.0) --
! 244: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 245: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 246: * November 2011
! 247: *
! 248: * .. Scalar Arguments ..
! 249: CHARACTER HOWMNY, SIDE
! 250: INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
! 251: * ..
! 252: * .. Array Arguments ..
! 253: LOGICAL SELECT( * )
! 254: DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
! 255: $ WORK( * )
! 256: * ..
! 257: *
! 258: * =====================================================================
! 259: *
! 260: * .. Parameters ..
! 261: DOUBLE PRECISION ZERO, ONE
! 262: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 263: INTEGER NBMIN, NBMAX
! 264: PARAMETER ( NBMIN = 8, NBMAX = 128 )
! 265: * ..
! 266: * .. Local Scalars ..
! 267: LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
! 268: $ RIGHTV, SOMEV
! 269: INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
! 270: $ IV, MAXWRK, NB, KI2
! 271: DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
! 272: $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
! 273: $ XNORM
! 274: * ..
! 275: * .. External Functions ..
! 276: LOGICAL LSAME
! 277: INTEGER IDAMAX, ILAENV
! 278: DOUBLE PRECISION DDOT, DLAMCH
! 279: EXTERNAL LSAME, IDAMAX, ILAENV, DDOT, DLAMCH
! 280: * ..
! 281: * .. External Subroutines ..
! 282: EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
! 283: * ..
! 284: * .. Intrinsic Functions ..
! 285: INTRINSIC ABS, MAX, SQRT
! 286: * ..
! 287: * .. Local Arrays ..
! 288: DOUBLE PRECISION X( 2, 2 )
! 289: INTEGER ISCOMPLEX( NBMAX )
! 290: * ..
! 291: * .. Executable Statements ..
! 292: *
! 293: * Decode and test the input parameters
! 294: *
! 295: BOTHV = LSAME( SIDE, 'B' )
! 296: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
! 297: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
! 298: *
! 299: ALLV = LSAME( HOWMNY, 'A' )
! 300: OVER = LSAME( HOWMNY, 'B' )
! 301: SOMEV = LSAME( HOWMNY, 'S' )
! 302: *
! 303: INFO = 0
! 304: NB = ILAENV( 1, 'DTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
! 305: MAXWRK = N + 2*N*NB
! 306: WORK(1) = MAXWRK
! 307: LQUERY = ( LWORK.EQ.-1 )
! 308: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
! 309: INFO = -1
! 310: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
! 311: INFO = -2
! 312: ELSE IF( N.LT.0 ) THEN
! 313: INFO = -4
! 314: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
! 315: INFO = -6
! 316: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
! 317: INFO = -8
! 318: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
! 319: INFO = -10
! 320: ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
! 321: INFO = -14
! 322: ELSE
! 323: *
! 324: * Set M to the number of columns required to store the selected
! 325: * eigenvectors, standardize the array SELECT if necessary, and
! 326: * test MM.
! 327: *
! 328: IF( SOMEV ) THEN
! 329: M = 0
! 330: PAIR = .FALSE.
! 331: DO 10 J = 1, N
! 332: IF( PAIR ) THEN
! 333: PAIR = .FALSE.
! 334: SELECT( J ) = .FALSE.
! 335: ELSE
! 336: IF( J.LT.N ) THEN
! 337: IF( T( J+1, J ).EQ.ZERO ) THEN
! 338: IF( SELECT( J ) )
! 339: $ M = M + 1
! 340: ELSE
! 341: PAIR = .TRUE.
! 342: IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
! 343: SELECT( J ) = .TRUE.
! 344: M = M + 2
! 345: END IF
! 346: END IF
! 347: ELSE
! 348: IF( SELECT( N ) )
! 349: $ M = M + 1
! 350: END IF
! 351: END IF
! 352: 10 CONTINUE
! 353: ELSE
! 354: M = N
! 355: END IF
! 356: *
! 357: IF( MM.LT.M ) THEN
! 358: INFO = -11
! 359: END IF
! 360: END IF
! 361: IF( INFO.NE.0 ) THEN
! 362: CALL XERBLA( 'DTREVC3', -INFO )
! 363: RETURN
! 364: ELSE IF( LQUERY ) THEN
! 365: RETURN
! 366: END IF
! 367: *
! 368: * Quick return if possible.
! 369: *
! 370: IF( N.EQ.0 )
! 371: $ RETURN
! 372: *
! 373: * Use blocked version of back-transformation if sufficient workspace.
! 374: * Zero-out the workspace to avoid potential NaN propagation.
! 375: *
! 376: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
! 377: NB = (LWORK - N) / (2*N)
! 378: NB = MIN( NB, NBMAX )
! 379: CALL DLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
! 380: ELSE
! 381: NB = 1
! 382: END IF
! 383: *
! 384: * Set the constants to control overflow.
! 385: *
! 386: UNFL = DLAMCH( 'Safe minimum' )
! 387: OVFL = ONE / UNFL
! 388: CALL DLABAD( UNFL, OVFL )
! 389: ULP = DLAMCH( 'Precision' )
! 390: SMLNUM = UNFL*( N / ULP )
! 391: BIGNUM = ( ONE-ULP ) / SMLNUM
! 392: *
! 393: * Compute 1-norm of each column of strictly upper triangular
! 394: * part of T to control overflow in triangular solver.
! 395: *
! 396: WORK( 1 ) = ZERO
! 397: DO 30 J = 2, N
! 398: WORK( J ) = ZERO
! 399: DO 20 I = 1, J - 1
! 400: WORK( J ) = WORK( J ) + ABS( T( I, J ) )
! 401: 20 CONTINUE
! 402: 30 CONTINUE
! 403: *
! 404: * Index IP is used to specify the real or complex eigenvalue:
! 405: * IP = 0, real eigenvalue,
! 406: * 1, first of conjugate complex pair: (wr,wi)
! 407: * -1, second of conjugate complex pair: (wr,wi)
! 408: * ISCOMPLEX array stores IP for each column in current block.
! 409: *
! 410: IF( RIGHTV ) THEN
! 411: *
! 412: * ============================================================
! 413: * Compute right eigenvectors.
! 414: *
! 415: * IV is index of column in current block.
! 416: * For complex right vector, uses IV-1 for real part and IV for complex part.
! 417: * Non-blocked version always uses IV=2;
! 418: * blocked version starts with IV=NB, goes down to 1 or 2.
! 419: * (Note the "0-th" column is used for 1-norms computed above.)
! 420: IV = 2
! 421: IF( NB.GT.2 ) THEN
! 422: IV = NB
! 423: END IF
! 424:
! 425: IP = 0
! 426: IS = M
! 427: DO 140 KI = N, 1, -1
! 428: IF( IP.EQ.-1 ) THEN
! 429: * previous iteration (ki+1) was second of conjugate pair,
! 430: * so this ki is first of conjugate pair; skip to end of loop
! 431: IP = 1
! 432: GO TO 140
! 433: ELSE IF( KI.EQ.1 ) THEN
! 434: * last column, so this ki must be real eigenvalue
! 435: IP = 0
! 436: ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
! 437: * zero on sub-diagonal, so this ki is real eigenvalue
! 438: IP = 0
! 439: ELSE
! 440: * non-zero on sub-diagonal, so this ki is second of conjugate pair
! 441: IP = -1
! 442: END IF
! 443:
! 444: IF( SOMEV ) THEN
! 445: IF( IP.EQ.0 ) THEN
! 446: IF( .NOT.SELECT( KI ) )
! 447: $ GO TO 140
! 448: ELSE
! 449: IF( .NOT.SELECT( KI-1 ) )
! 450: $ GO TO 140
! 451: END IF
! 452: END IF
! 453: *
! 454: * Compute the KI-th eigenvalue (WR,WI).
! 455: *
! 456: WR = T( KI, KI )
! 457: WI = ZERO
! 458: IF( IP.NE.0 )
! 459: $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
! 460: $ SQRT( ABS( T( KI-1, KI ) ) )
! 461: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
! 462: *
! 463: IF( IP.EQ.0 ) THEN
! 464: *
! 465: * --------------------------------------------------------
! 466: * Real right eigenvector
! 467: *
! 468: WORK( KI + IV*N ) = ONE
! 469: *
! 470: * Form right-hand side.
! 471: *
! 472: DO 50 K = 1, KI - 1
! 473: WORK( K + IV*N ) = -T( K, KI )
! 474: 50 CONTINUE
! 475: *
! 476: * Solve upper quasi-triangular system:
! 477: * [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
! 478: *
! 479: JNXT = KI - 1
! 480: DO 60 J = KI - 1, 1, -1
! 481: IF( J.GT.JNXT )
! 482: $ GO TO 60
! 483: J1 = J
! 484: J2 = J
! 485: JNXT = J - 1
! 486: IF( J.GT.1 ) THEN
! 487: IF( T( J, J-1 ).NE.ZERO ) THEN
! 488: J1 = J - 1
! 489: JNXT = J - 2
! 490: END IF
! 491: END IF
! 492: *
! 493: IF( J1.EQ.J2 ) THEN
! 494: *
! 495: * 1-by-1 diagonal block
! 496: *
! 497: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
! 498: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
! 499: $ ZERO, X, 2, SCALE, XNORM, IERR )
! 500: *
! 501: * Scale X(1,1) to avoid overflow when updating
! 502: * the right-hand side.
! 503: *
! 504: IF( XNORM.GT.ONE ) THEN
! 505: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
! 506: X( 1, 1 ) = X( 1, 1 ) / XNORM
! 507: SCALE = SCALE / XNORM
! 508: END IF
! 509: END IF
! 510: *
! 511: * Scale if necessary
! 512: *
! 513: IF( SCALE.NE.ONE )
! 514: $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
! 515: WORK( J+IV*N ) = X( 1, 1 )
! 516: *
! 517: * Update right-hand side
! 518: *
! 519: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
! 520: $ WORK( 1+IV*N ), 1 )
! 521: *
! 522: ELSE
! 523: *
! 524: * 2-by-2 diagonal block
! 525: *
! 526: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
! 527: $ T( J-1, J-1 ), LDT, ONE, ONE,
! 528: $ WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
! 529: $ SCALE, XNORM, IERR )
! 530: *
! 531: * Scale X(1,1) and X(2,1) to avoid overflow when
! 532: * updating the right-hand side.
! 533: *
! 534: IF( XNORM.GT.ONE ) THEN
! 535: BETA = MAX( WORK( J-1 ), WORK( J ) )
! 536: IF( BETA.GT.BIGNUM / XNORM ) THEN
! 537: X( 1, 1 ) = X( 1, 1 ) / XNORM
! 538: X( 2, 1 ) = X( 2, 1 ) / XNORM
! 539: SCALE = SCALE / XNORM
! 540: END IF
! 541: END IF
! 542: *
! 543: * Scale if necessary
! 544: *
! 545: IF( SCALE.NE.ONE )
! 546: $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
! 547: WORK( J-1+IV*N ) = X( 1, 1 )
! 548: WORK( J +IV*N ) = X( 2, 1 )
! 549: *
! 550: * Update right-hand side
! 551: *
! 552: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
! 553: $ WORK( 1+IV*N ), 1 )
! 554: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
! 555: $ WORK( 1+IV*N ), 1 )
! 556: END IF
! 557: 60 CONTINUE
! 558: *
! 559: * Copy the vector x or Q*x to VR and normalize.
! 560: *
! 561: IF( .NOT.OVER ) THEN
! 562: * ------------------------------
! 563: * no back-transform: copy x to VR and normalize.
! 564: CALL DCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
! 565: *
! 566: II = IDAMAX( KI, VR( 1, IS ), 1 )
! 567: REMAX = ONE / ABS( VR( II, IS ) )
! 568: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
! 569: *
! 570: DO 70 K = KI + 1, N
! 571: VR( K, IS ) = ZERO
! 572: 70 CONTINUE
! 573: *
! 574: ELSE IF( NB.EQ.1 ) THEN
! 575: * ------------------------------
! 576: * version 1: back-transform each vector with GEMV, Q*x.
! 577: IF( KI.GT.1 )
! 578: $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
! 579: $ WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
! 580: $ VR( 1, KI ), 1 )
! 581: *
! 582: II = IDAMAX( N, VR( 1, KI ), 1 )
! 583: REMAX = ONE / ABS( VR( II, KI ) )
! 584: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
! 585: *
! 586: ELSE
! 587: * ------------------------------
! 588: * version 2: back-transform block of vectors with GEMM
! 589: * zero out below vector
! 590: DO K = KI + 1, N
! 591: WORK( K + IV*N ) = ZERO
! 592: END DO
! 593: ISCOMPLEX( IV ) = IP
! 594: * back-transform and normalization is done below
! 595: END IF
! 596: ELSE
! 597: *
! 598: * --------------------------------------------------------
! 599: * Complex right eigenvector.
! 600: *
! 601: * Initial solve
! 602: * [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
! 603: * [ ( T(KI, KI-1) T(KI, KI) ) ]
! 604: *
! 605: IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
! 606: WORK( KI-1 + (IV-1)*N ) = ONE
! 607: WORK( KI + (IV )*N ) = WI / T( KI-1, KI )
! 608: ELSE
! 609: WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
! 610: WORK( KI + (IV )*N ) = ONE
! 611: END IF
! 612: WORK( KI + (IV-1)*N ) = ZERO
! 613: WORK( KI-1 + (IV )*N ) = ZERO
! 614: *
! 615: * Form right-hand side.
! 616: *
! 617: DO 80 K = 1, KI - 2
! 618: WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
! 619: WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(K,KI )
! 620: 80 CONTINUE
! 621: *
! 622: * Solve upper quasi-triangular system:
! 623: * [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
! 624: *
! 625: JNXT = KI - 2
! 626: DO 90 J = KI - 2, 1, -1
! 627: IF( J.GT.JNXT )
! 628: $ GO TO 90
! 629: J1 = J
! 630: J2 = J
! 631: JNXT = J - 1
! 632: IF( J.GT.1 ) THEN
! 633: IF( T( J, J-1 ).NE.ZERO ) THEN
! 634: J1 = J - 1
! 635: JNXT = J - 2
! 636: END IF
! 637: END IF
! 638: *
! 639: IF( J1.EQ.J2 ) THEN
! 640: *
! 641: * 1-by-1 diagonal block
! 642: *
! 643: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
! 644: $ LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
! 645: $ WR, WI, X, 2, SCALE, XNORM, IERR )
! 646: *
! 647: * Scale X(1,1) and X(1,2) to avoid overflow when
! 648: * updating the right-hand side.
! 649: *
! 650: IF( XNORM.GT.ONE ) THEN
! 651: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
! 652: X( 1, 1 ) = X( 1, 1 ) / XNORM
! 653: X( 1, 2 ) = X( 1, 2 ) / XNORM
! 654: SCALE = SCALE / XNORM
! 655: END IF
! 656: END IF
! 657: *
! 658: * Scale if necessary
! 659: *
! 660: IF( SCALE.NE.ONE ) THEN
! 661: CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
! 662: CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
! 663: END IF
! 664: WORK( J+(IV-1)*N ) = X( 1, 1 )
! 665: WORK( J+(IV )*N ) = X( 1, 2 )
! 666: *
! 667: * Update the right-hand side
! 668: *
! 669: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
! 670: $ WORK( 1+(IV-1)*N ), 1 )
! 671: CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
! 672: $ WORK( 1+(IV )*N ), 1 )
! 673: *
! 674: ELSE
! 675: *
! 676: * 2-by-2 diagonal block
! 677: *
! 678: CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
! 679: $ T( J-1, J-1 ), LDT, ONE, ONE,
! 680: $ WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
! 681: $ SCALE, XNORM, IERR )
! 682: *
! 683: * Scale X to avoid overflow when updating
! 684: * the right-hand side.
! 685: *
! 686: IF( XNORM.GT.ONE ) THEN
! 687: BETA = MAX( WORK( J-1 ), WORK( J ) )
! 688: IF( BETA.GT.BIGNUM / XNORM ) THEN
! 689: REC = ONE / XNORM
! 690: X( 1, 1 ) = X( 1, 1 )*REC
! 691: X( 1, 2 ) = X( 1, 2 )*REC
! 692: X( 2, 1 ) = X( 2, 1 )*REC
! 693: X( 2, 2 ) = X( 2, 2 )*REC
! 694: SCALE = SCALE*REC
! 695: END IF
! 696: END IF
! 697: *
! 698: * Scale if necessary
! 699: *
! 700: IF( SCALE.NE.ONE ) THEN
! 701: CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
! 702: CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
! 703: END IF
! 704: WORK( J-1+(IV-1)*N ) = X( 1, 1 )
! 705: WORK( J +(IV-1)*N ) = X( 2, 1 )
! 706: WORK( J-1+(IV )*N ) = X( 1, 2 )
! 707: WORK( J +(IV )*N ) = X( 2, 2 )
! 708: *
! 709: * Update the right-hand side
! 710: *
! 711: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
! 712: $ WORK( 1+(IV-1)*N ), 1 )
! 713: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
! 714: $ WORK( 1+(IV-1)*N ), 1 )
! 715: CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
! 716: $ WORK( 1+(IV )*N ), 1 )
! 717: CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
! 718: $ WORK( 1+(IV )*N ), 1 )
! 719: END IF
! 720: 90 CONTINUE
! 721: *
! 722: * Copy the vector x or Q*x to VR and normalize.
! 723: *
! 724: IF( .NOT.OVER ) THEN
! 725: * ------------------------------
! 726: * no back-transform: copy x to VR and normalize.
! 727: CALL DCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
! 728: CALL DCOPY( KI, WORK( 1+(IV )*N ), 1, VR(1,IS ), 1 )
! 729: *
! 730: EMAX = ZERO
! 731: DO 100 K = 1, KI
! 732: EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
! 733: $ ABS( VR( K, IS ) ) )
! 734: 100 CONTINUE
! 735: REMAX = ONE / EMAX
! 736: CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
! 737: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
! 738: *
! 739: DO 110 K = KI + 1, N
! 740: VR( K, IS-1 ) = ZERO
! 741: VR( K, IS ) = ZERO
! 742: 110 CONTINUE
! 743: *
! 744: ELSE IF( NB.EQ.1 ) THEN
! 745: * ------------------------------
! 746: * version 1: back-transform each vector with GEMV, Q*x.
! 747: IF( KI.GT.2 ) THEN
! 748: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
! 749: $ WORK( 1 + (IV-1)*N ), 1,
! 750: $ WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
! 751: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
! 752: $ WORK( 1 + (IV)*N ), 1,
! 753: $ WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
! 754: ELSE
! 755: CALL DSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
! 756: CALL DSCAL( N, WORK(KI +(IV )*N), VR(1,KI ), 1)
! 757: END IF
! 758: *
! 759: EMAX = ZERO
! 760: DO 120 K = 1, N
! 761: EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
! 762: $ ABS( VR( K, KI ) ) )
! 763: 120 CONTINUE
! 764: REMAX = ONE / EMAX
! 765: CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
! 766: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
! 767: *
! 768: ELSE
! 769: * ------------------------------
! 770: * version 2: back-transform block of vectors with GEMM
! 771: * zero out below vector
! 772: DO K = KI + 1, N
! 773: WORK( K + (IV-1)*N ) = ZERO
! 774: WORK( K + (IV )*N ) = ZERO
! 775: END DO
! 776: ISCOMPLEX( IV-1 ) = -IP
! 777: ISCOMPLEX( IV ) = IP
! 778: IV = IV - 1
! 779: * back-transform and normalization is done below
! 780: END IF
! 781: END IF
! 782:
! 783: IF( NB.GT.1 ) THEN
! 784: * --------------------------------------------------------
! 785: * Blocked version of back-transform
! 786: * For complex case, KI2 includes both vectors (KI-1 and KI)
! 787: IF( IP.EQ.0 ) THEN
! 788: KI2 = KI
! 789: ELSE
! 790: KI2 = KI - 1
! 791: END IF
! 792:
! 793: * Columns IV:NB of work are valid vectors.
! 794: * When the number of vectors stored reaches NB-1 or NB,
! 795: * or if this was last vector, do the GEMM
! 796: IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
! 797: CALL DGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
! 798: $ VR, LDVR,
! 799: $ WORK( 1 + (IV)*N ), N,
! 800: $ ZERO,
! 801: $ WORK( 1 + (NB+IV)*N ), N )
! 802: * normalize vectors
! 803: DO K = IV, NB
! 804: IF( ISCOMPLEX(K).EQ.0 ) THEN
! 805: * real eigenvector
! 806: II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
! 807: REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
! 808: ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
! 809: * first eigenvector of conjugate pair
! 810: EMAX = ZERO
! 811: DO II = 1, N
! 812: EMAX = MAX( EMAX,
! 813: $ ABS( WORK( II + (NB+K )*N ) )+
! 814: $ ABS( WORK( II + (NB+K+1)*N ) ) )
! 815: END DO
! 816: REMAX = ONE / EMAX
! 817: * else if ISCOMPLEX(K).EQ.-1
! 818: * second eigenvector of conjugate pair
! 819: * reuse same REMAX as previous K
! 820: END IF
! 821: CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
! 822: END DO
! 823: CALL DLACPY( 'F', N, NB-IV+1,
! 824: $ WORK( 1 + (NB+IV)*N ), N,
! 825: $ VR( 1, KI2 ), LDVR )
! 826: IV = NB
! 827: ELSE
! 828: IV = IV - 1
! 829: END IF
! 830: END IF ! blocked back-transform
! 831: *
! 832: IS = IS - 1
! 833: IF( IP.NE.0 )
! 834: $ IS = IS - 1
! 835: 140 CONTINUE
! 836: END IF
! 837:
! 838: IF( LEFTV ) THEN
! 839: *
! 840: * ============================================================
! 841: * Compute left eigenvectors.
! 842: *
! 843: * IV is index of column in current block.
! 844: * For complex left vector, uses IV for real part and IV+1 for complex part.
! 845: * Non-blocked version always uses IV=1;
! 846: * blocked version starts with IV=1, goes up to NB-1 or NB.
! 847: * (Note the "0-th" column is used for 1-norms computed above.)
! 848: IV = 1
! 849: IP = 0
! 850: IS = 1
! 851: DO 260 KI = 1, N
! 852: IF( IP.EQ.1 ) THEN
! 853: * previous iteration (ki-1) was first of conjugate pair,
! 854: * so this ki is second of conjugate pair; skip to end of loop
! 855: IP = -1
! 856: GO TO 260
! 857: ELSE IF( KI.EQ.N ) THEN
! 858: * last column, so this ki must be real eigenvalue
! 859: IP = 0
! 860: ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
! 861: * zero on sub-diagonal, so this ki is real eigenvalue
! 862: IP = 0
! 863: ELSE
! 864: * non-zero on sub-diagonal, so this ki is first of conjugate pair
! 865: IP = 1
! 866: END IF
! 867: *
! 868: IF( SOMEV ) THEN
! 869: IF( .NOT.SELECT( KI ) )
! 870: $ GO TO 260
! 871: END IF
! 872: *
! 873: * Compute the KI-th eigenvalue (WR,WI).
! 874: *
! 875: WR = T( KI, KI )
! 876: WI = ZERO
! 877: IF( IP.NE.0 )
! 878: $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
! 879: $ SQRT( ABS( T( KI+1, KI ) ) )
! 880: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
! 881: *
! 882: IF( IP.EQ.0 ) THEN
! 883: *
! 884: * --------------------------------------------------------
! 885: * Real left eigenvector
! 886: *
! 887: WORK( KI + IV*N ) = ONE
! 888: *
! 889: * Form right-hand side.
! 890: *
! 891: DO 160 K = KI + 1, N
! 892: WORK( K + IV*N ) = -T( KI, K )
! 893: 160 CONTINUE
! 894: *
! 895: * Solve transposed quasi-triangular system:
! 896: * [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
! 897: *
! 898: VMAX = ONE
! 899: VCRIT = BIGNUM
! 900: *
! 901: JNXT = KI + 1
! 902: DO 170 J = KI + 1, N
! 903: IF( J.LT.JNXT )
! 904: $ GO TO 170
! 905: J1 = J
! 906: J2 = J
! 907: JNXT = J + 1
! 908: IF( J.LT.N ) THEN
! 909: IF( T( J+1, J ).NE.ZERO ) THEN
! 910: J2 = J + 1
! 911: JNXT = J + 2
! 912: END IF
! 913: END IF
! 914: *
! 915: IF( J1.EQ.J2 ) THEN
! 916: *
! 917: * 1-by-1 diagonal block
! 918: *
! 919: * Scale if necessary to avoid overflow when forming
! 920: * the right-hand side.
! 921: *
! 922: IF( WORK( J ).GT.VCRIT ) THEN
! 923: REC = ONE / VMAX
! 924: CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
! 925: VMAX = ONE
! 926: VCRIT = BIGNUM
! 927: END IF
! 928: *
! 929: WORK( J+IV*N ) = WORK( J+IV*N ) -
! 930: $ DDOT( J-KI-1, T( KI+1, J ), 1,
! 931: $ WORK( KI+1+IV*N ), 1 )
! 932: *
! 933: * Solve [ T(J,J) - WR ]**T * X = WORK
! 934: *
! 935: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
! 936: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
! 937: $ ZERO, X, 2, SCALE, XNORM, IERR )
! 938: *
! 939: * Scale if necessary
! 940: *
! 941: IF( SCALE.NE.ONE )
! 942: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
! 943: WORK( J+IV*N ) = X( 1, 1 )
! 944: VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
! 945: VCRIT = BIGNUM / VMAX
! 946: *
! 947: ELSE
! 948: *
! 949: * 2-by-2 diagonal block
! 950: *
! 951: * Scale if necessary to avoid overflow when forming
! 952: * the right-hand side.
! 953: *
! 954: BETA = MAX( WORK( J ), WORK( J+1 ) )
! 955: IF( BETA.GT.VCRIT ) THEN
! 956: REC = ONE / VMAX
! 957: CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
! 958: VMAX = ONE
! 959: VCRIT = BIGNUM
! 960: END IF
! 961: *
! 962: WORK( J+IV*N ) = WORK( J+IV*N ) -
! 963: $ DDOT( J-KI-1, T( KI+1, J ), 1,
! 964: $ WORK( KI+1+IV*N ), 1 )
! 965: *
! 966: WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
! 967: $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
! 968: $ WORK( KI+1+IV*N ), 1 )
! 969: *
! 970: * Solve
! 971: * [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
! 972: * [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
! 973: *
! 974: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
! 975: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
! 976: $ ZERO, X, 2, SCALE, XNORM, IERR )
! 977: *
! 978: * Scale if necessary
! 979: *
! 980: IF( SCALE.NE.ONE )
! 981: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
! 982: WORK( J +IV*N ) = X( 1, 1 )
! 983: WORK( J+1+IV*N ) = X( 2, 1 )
! 984: *
! 985: VMAX = MAX( ABS( WORK( J +IV*N ) ),
! 986: $ ABS( WORK( J+1+IV*N ) ), VMAX )
! 987: VCRIT = BIGNUM / VMAX
! 988: *
! 989: END IF
! 990: 170 CONTINUE
! 991: *
! 992: * Copy the vector x or Q*x to VL and normalize.
! 993: *
! 994: IF( .NOT.OVER ) THEN
! 995: * ------------------------------
! 996: * no back-transform: copy x to VL and normalize.
! 997: CALL DCOPY( N-KI+1, WORK( KI + IV*N ), 1,
! 998: $ VL( KI, IS ), 1 )
! 999: *
! 1000: II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
! 1001: REMAX = ONE / ABS( VL( II, IS ) )
! 1002: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
! 1003: *
! 1004: DO 180 K = 1, KI - 1
! 1005: VL( K, IS ) = ZERO
! 1006: 180 CONTINUE
! 1007: *
! 1008: ELSE IF( NB.EQ.1 ) THEN
! 1009: * ------------------------------
! 1010: * version 1: back-transform each vector with GEMV, Q*x.
! 1011: IF( KI.LT.N )
! 1012: $ CALL DGEMV( 'N', N, N-KI, ONE,
! 1013: $ VL( 1, KI+1 ), LDVL,
! 1014: $ WORK( KI+1 + IV*N ), 1,
! 1015: $ WORK( KI + IV*N ), VL( 1, KI ), 1 )
! 1016: *
! 1017: II = IDAMAX( N, VL( 1, KI ), 1 )
! 1018: REMAX = ONE / ABS( VL( II, KI ) )
! 1019: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
! 1020: *
! 1021: ELSE
! 1022: * ------------------------------
! 1023: * version 2: back-transform block of vectors with GEMM
! 1024: * zero out above vector
! 1025: * could go from KI-NV+1 to KI-1
! 1026: DO K = 1, KI - 1
! 1027: WORK( K + IV*N ) = ZERO
! 1028: END DO
! 1029: ISCOMPLEX( IV ) = IP
! 1030: * back-transform and normalization is done below
! 1031: END IF
! 1032: ELSE
! 1033: *
! 1034: * --------------------------------------------------------
! 1035: * Complex left eigenvector.
! 1036: *
! 1037: * Initial solve:
! 1038: * [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
! 1039: * [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
! 1040: *
! 1041: IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
! 1042: WORK( KI + (IV )*N ) = WI / T( KI, KI+1 )
! 1043: WORK( KI+1 + (IV+1)*N ) = ONE
! 1044: ELSE
! 1045: WORK( KI + (IV )*N ) = ONE
! 1046: WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
! 1047: END IF
! 1048: WORK( KI+1 + (IV )*N ) = ZERO
! 1049: WORK( KI + (IV+1)*N ) = ZERO
! 1050: *
! 1051: * Form right-hand side.
! 1052: *
! 1053: DO 190 K = KI + 2, N
! 1054: WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(KI, K)
! 1055: WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
! 1056: 190 CONTINUE
! 1057: *
! 1058: * Solve transposed quasi-triangular system:
! 1059: * [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
! 1060: *
! 1061: VMAX = ONE
! 1062: VCRIT = BIGNUM
! 1063: *
! 1064: JNXT = KI + 2
! 1065: DO 200 J = KI + 2, N
! 1066: IF( J.LT.JNXT )
! 1067: $ GO TO 200
! 1068: J1 = J
! 1069: J2 = J
! 1070: JNXT = J + 1
! 1071: IF( J.LT.N ) THEN
! 1072: IF( T( J+1, J ).NE.ZERO ) THEN
! 1073: J2 = J + 1
! 1074: JNXT = J + 2
! 1075: END IF
! 1076: END IF
! 1077: *
! 1078: IF( J1.EQ.J2 ) THEN
! 1079: *
! 1080: * 1-by-1 diagonal block
! 1081: *
! 1082: * Scale if necessary to avoid overflow when
! 1083: * forming the right-hand side elements.
! 1084: *
! 1085: IF( WORK( J ).GT.VCRIT ) THEN
! 1086: REC = ONE / VMAX
! 1087: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
! 1088: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
! 1089: VMAX = ONE
! 1090: VCRIT = BIGNUM
! 1091: END IF
! 1092: *
! 1093: WORK( J+(IV )*N ) = WORK( J+(IV)*N ) -
! 1094: $ DDOT( J-KI-2, T( KI+2, J ), 1,
! 1095: $ WORK( KI+2+(IV)*N ), 1 )
! 1096: WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
! 1097: $ DDOT( J-KI-2, T( KI+2, J ), 1,
! 1098: $ WORK( KI+2+(IV+1)*N ), 1 )
! 1099: *
! 1100: * Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
! 1101: *
! 1102: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
! 1103: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
! 1104: $ -WI, X, 2, SCALE, XNORM, IERR )
! 1105: *
! 1106: * Scale if necessary
! 1107: *
! 1108: IF( SCALE.NE.ONE ) THEN
! 1109: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
! 1110: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
! 1111: END IF
! 1112: WORK( J+(IV )*N ) = X( 1, 1 )
! 1113: WORK( J+(IV+1)*N ) = X( 1, 2 )
! 1114: VMAX = MAX( ABS( WORK( J+(IV )*N ) ),
! 1115: $ ABS( WORK( J+(IV+1)*N ) ), VMAX )
! 1116: VCRIT = BIGNUM / VMAX
! 1117: *
! 1118: ELSE
! 1119: *
! 1120: * 2-by-2 diagonal block
! 1121: *
! 1122: * Scale if necessary to avoid overflow when forming
! 1123: * the right-hand side elements.
! 1124: *
! 1125: BETA = MAX( WORK( J ), WORK( J+1 ) )
! 1126: IF( BETA.GT.VCRIT ) THEN
! 1127: REC = ONE / VMAX
! 1128: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
! 1129: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
! 1130: VMAX = ONE
! 1131: VCRIT = BIGNUM
! 1132: END IF
! 1133: *
! 1134: WORK( J +(IV )*N ) = WORK( J+(IV)*N ) -
! 1135: $ DDOT( J-KI-2, T( KI+2, J ), 1,
! 1136: $ WORK( KI+2+(IV)*N ), 1 )
! 1137: *
! 1138: WORK( J +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
! 1139: $ DDOT( J-KI-2, T( KI+2, J ), 1,
! 1140: $ WORK( KI+2+(IV+1)*N ), 1 )
! 1141: *
! 1142: WORK( J+1+(IV )*N ) = WORK( J+1+(IV)*N ) -
! 1143: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
! 1144: $ WORK( KI+2+(IV)*N ), 1 )
! 1145: *
! 1146: WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
! 1147: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
! 1148: $ WORK( KI+2+(IV+1)*N ), 1 )
! 1149: *
! 1150: * Solve 2-by-2 complex linear equation
! 1151: * [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
! 1152: * [ (T(j+1,j) T(j+1,j+1)) ]
! 1153: *
! 1154: CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
! 1155: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
! 1156: $ -WI, X, 2, SCALE, XNORM, IERR )
! 1157: *
! 1158: * Scale if necessary
! 1159: *
! 1160: IF( SCALE.NE.ONE ) THEN
! 1161: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
! 1162: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
! 1163: END IF
! 1164: WORK( J +(IV )*N ) = X( 1, 1 )
! 1165: WORK( J +(IV+1)*N ) = X( 1, 2 )
! 1166: WORK( J+1+(IV )*N ) = X( 2, 1 )
! 1167: WORK( J+1+(IV+1)*N ) = X( 2, 2 )
! 1168: VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
! 1169: $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
! 1170: $ VMAX )
! 1171: VCRIT = BIGNUM / VMAX
! 1172: *
! 1173: END IF
! 1174: 200 CONTINUE
! 1175: *
! 1176: * Copy the vector x or Q*x to VL and normalize.
! 1177: *
! 1178: IF( .NOT.OVER ) THEN
! 1179: * ------------------------------
! 1180: * no back-transform: copy x to VL and normalize.
! 1181: CALL DCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
! 1182: $ VL( KI, IS ), 1 )
! 1183: CALL DCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
! 1184: $ VL( KI, IS+1 ), 1 )
! 1185: *
! 1186: EMAX = ZERO
! 1187: DO 220 K = KI, N
! 1188: EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
! 1189: $ ABS( VL( K, IS+1 ) ) )
! 1190: 220 CONTINUE
! 1191: REMAX = ONE / EMAX
! 1192: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
! 1193: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
! 1194: *
! 1195: DO 230 K = 1, KI - 1
! 1196: VL( K, IS ) = ZERO
! 1197: VL( K, IS+1 ) = ZERO
! 1198: 230 CONTINUE
! 1199: *
! 1200: ELSE IF( NB.EQ.1 ) THEN
! 1201: * ------------------------------
! 1202: * version 1: back-transform each vector with GEMV, Q*x.
! 1203: IF( KI.LT.N-1 ) THEN
! 1204: CALL DGEMV( 'N', N, N-KI-1, ONE,
! 1205: $ VL( 1, KI+2 ), LDVL,
! 1206: $ WORK( KI+2 + (IV)*N ), 1,
! 1207: $ WORK( KI + (IV)*N ),
! 1208: $ VL( 1, KI ), 1 )
! 1209: CALL DGEMV( 'N', N, N-KI-1, ONE,
! 1210: $ VL( 1, KI+2 ), LDVL,
! 1211: $ WORK( KI+2 + (IV+1)*N ), 1,
! 1212: $ WORK( KI+1 + (IV+1)*N ),
! 1213: $ VL( 1, KI+1 ), 1 )
! 1214: ELSE
! 1215: CALL DSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
! 1216: CALL DSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
! 1217: END IF
! 1218: *
! 1219: EMAX = ZERO
! 1220: DO 240 K = 1, N
! 1221: EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
! 1222: $ ABS( VL( K, KI+1 ) ) )
! 1223: 240 CONTINUE
! 1224: REMAX = ONE / EMAX
! 1225: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
! 1226: CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
! 1227: *
! 1228: ELSE
! 1229: * ------------------------------
! 1230: * version 2: back-transform block of vectors with GEMM
! 1231: * zero out above vector
! 1232: * could go from KI-NV+1 to KI-1
! 1233: DO K = 1, KI - 1
! 1234: WORK( K + (IV )*N ) = ZERO
! 1235: WORK( K + (IV+1)*N ) = ZERO
! 1236: END DO
! 1237: ISCOMPLEX( IV ) = IP
! 1238: ISCOMPLEX( IV+1 ) = -IP
! 1239: IV = IV + 1
! 1240: * back-transform and normalization is done below
! 1241: END IF
! 1242: END IF
! 1243:
! 1244: IF( NB.GT.1 ) THEN
! 1245: * --------------------------------------------------------
! 1246: * Blocked version of back-transform
! 1247: * For complex case, KI2 includes both vectors (KI and KI+1)
! 1248: IF( IP.EQ.0 ) THEN
! 1249: KI2 = KI
! 1250: ELSE
! 1251: KI2 = KI + 1
! 1252: END IF
! 1253:
! 1254: * Columns 1:IV of work are valid vectors.
! 1255: * When the number of vectors stored reaches NB-1 or NB,
! 1256: * or if this was last vector, do the GEMM
! 1257: IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
! 1258: CALL DGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
! 1259: $ VL( 1, KI2-IV+1 ), LDVL,
! 1260: $ WORK( KI2-IV+1 + (1)*N ), N,
! 1261: $ ZERO,
! 1262: $ WORK( 1 + (NB+1)*N ), N )
! 1263: * normalize vectors
! 1264: DO K = 1, IV
! 1265: IF( ISCOMPLEX(K).EQ.0) THEN
! 1266: * real eigenvector
! 1267: II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
! 1268: REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
! 1269: ELSE IF( ISCOMPLEX(K).EQ.1) THEN
! 1270: * first eigenvector of conjugate pair
! 1271: EMAX = ZERO
! 1272: DO II = 1, N
! 1273: EMAX = MAX( EMAX,
! 1274: $ ABS( WORK( II + (NB+K )*N ) )+
! 1275: $ ABS( WORK( II + (NB+K+1)*N ) ) )
! 1276: END DO
! 1277: REMAX = ONE / EMAX
! 1278: * else if ISCOMPLEX(K).EQ.-1
! 1279: * second eigenvector of conjugate pair
! 1280: * reuse same REMAX as previous K
! 1281: END IF
! 1282: CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
! 1283: END DO
! 1284: CALL DLACPY( 'F', N, IV,
! 1285: $ WORK( 1 + (NB+1)*N ), N,
! 1286: $ VL( 1, KI2-IV+1 ), LDVL )
! 1287: IV = 1
! 1288: ELSE
! 1289: IV = IV + 1
! 1290: END IF
! 1291: END IF ! blocked back-transform
! 1292: *
! 1293: IS = IS + 1
! 1294: IF( IP.NE.0 )
! 1295: $ IS = IS + 1
! 1296: 260 CONTINUE
! 1297: END IF
! 1298: *
! 1299: RETURN
! 1300: *
! 1301: * End of DTREVC3
! 1302: *
! 1303: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrevc3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack