Annotation of rpl/lapack/lapack/dtgevc.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
! 2: $ LDVL, VR, 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, LDP, LDS, LDVL, LDVR, M, MM, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: LOGICAL SELECT( * )
! 15: DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
! 16: $ VR( LDVR, * ), WORK( * )
! 17: * ..
! 18: *
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * DTGEVC computes some or all of the right and/or left eigenvectors of
! 24: * a pair of real matrices (S,P), where S is a quasi-triangular matrix
! 25: * and P is upper triangular. Matrix pairs of this type are produced by
! 26: * the generalized Schur factorization of a matrix pair (A,B):
! 27: *
! 28: * A = Q*S*Z**T, B = Q*P*Z**T
! 29: *
! 30: * as computed by DGGHRD + DHGEQZ.
! 31: *
! 32: * The right eigenvector x and the left eigenvector y of (S,P)
! 33: * corresponding to an eigenvalue w are defined by:
! 34: *
! 35: * S*x = w*P*x, (y**H)*S = w*(y**H)*P,
! 36: *
! 37: * where y**H denotes the conjugate tranpose of y.
! 38: * The eigenvalues are not input to this routine, but are computed
! 39: * directly from the diagonal blocks of S and P.
! 40: *
! 41: * This routine returns the matrices X and/or Y of right and left
! 42: * eigenvectors of (S,P), or the products Z*X and/or Q*Y,
! 43: * where Z and Q are input matrices.
! 44: * If Q and Z are the orthogonal factors from the generalized Schur
! 45: * factorization of a matrix pair (A,B), then Z*X and Q*Y
! 46: * are the matrices of right and left eigenvectors of (A,B).
! 47: *
! 48: * Arguments
! 49: * =========
! 50: *
! 51: * SIDE (input) CHARACTER*1
! 52: * = 'R': compute right eigenvectors only;
! 53: * = 'L': compute left eigenvectors only;
! 54: * = 'B': compute both right and left eigenvectors.
! 55: *
! 56: * HOWMNY (input) CHARACTER*1
! 57: * = 'A': compute all right and/or left eigenvectors;
! 58: * = 'B': compute all right and/or left eigenvectors,
! 59: * backtransformed by the matrices in VR and/or VL;
! 60: * = 'S': compute selected right and/or left eigenvectors,
! 61: * specified by the logical array SELECT.
! 62: *
! 63: * SELECT (input) LOGICAL array, dimension (N)
! 64: * If HOWMNY='S', SELECT specifies the eigenvectors to be
! 65: * computed. If w(j) is a real eigenvalue, the corresponding
! 66: * real eigenvector is computed if SELECT(j) is .TRUE..
! 67: * If w(j) and w(j+1) are the real and imaginary parts of a
! 68: * complex eigenvalue, the corresponding complex eigenvector
! 69: * is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
! 70: * and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
! 71: * set to .FALSE..
! 72: * Not referenced if HOWMNY = 'A' or 'B'.
! 73: *
! 74: * N (input) INTEGER
! 75: * The order of the matrices S and P. N >= 0.
! 76: *
! 77: * S (input) DOUBLE PRECISION array, dimension (LDS,N)
! 78: * The upper quasi-triangular matrix S from a generalized Schur
! 79: * factorization, as computed by DHGEQZ.
! 80: *
! 81: * LDS (input) INTEGER
! 82: * The leading dimension of array S. LDS >= max(1,N).
! 83: *
! 84: * P (input) DOUBLE PRECISION array, dimension (LDP,N)
! 85: * The upper triangular matrix P from a generalized Schur
! 86: * factorization, as computed by DHGEQZ.
! 87: * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
! 88: * of S must be in positive diagonal form.
! 89: *
! 90: * LDP (input) INTEGER
! 91: * The leading dimension of array P. LDP >= max(1,N).
! 92: *
! 93: * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
! 94: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
! 95: * contain an N-by-N matrix Q (usually the orthogonal matrix Q
! 96: * of left Schur vectors returned by DHGEQZ).
! 97: * On exit, if SIDE = 'L' or 'B', VL contains:
! 98: * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
! 99: * if HOWMNY = 'B', the matrix Q*Y;
! 100: * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
! 101: * SELECT, stored consecutively in the columns of
! 102: * VL, in the same order as their eigenvalues.
! 103: *
! 104: * A complex eigenvector corresponding to a complex eigenvalue
! 105: * is stored in two consecutive columns, the first holding the
! 106: * real part, and the second the imaginary part.
! 107: *
! 108: * Not referenced if SIDE = 'R'.
! 109: *
! 110: * LDVL (input) INTEGER
! 111: * The leading dimension of array VL. LDVL >= 1, and if
! 112: * SIDE = 'L' or 'B', LDVL >= N.
! 113: *
! 114: * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
! 115: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
! 116: * contain an N-by-N matrix Z (usually the orthogonal matrix Z
! 117: * of right Schur vectors returned by DHGEQZ).
! 118: *
! 119: * On exit, if SIDE = 'R' or 'B', VR contains:
! 120: * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
! 121: * if HOWMNY = 'B' or 'b', the matrix Z*X;
! 122: * if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
! 123: * specified by SELECT, stored consecutively in the
! 124: * columns of VR, in the same order as their
! 125: * eigenvalues.
! 126: *
! 127: * A complex eigenvector corresponding to a complex eigenvalue
! 128: * is stored in two consecutive columns, the first holding the
! 129: * real part and the second the imaginary part.
! 130: *
! 131: * Not referenced if SIDE = 'L'.
! 132: *
! 133: * LDVR (input) INTEGER
! 134: * The leading dimension of the array VR. LDVR >= 1, and if
! 135: * SIDE = 'R' or 'B', LDVR >= N.
! 136: *
! 137: * MM (input) INTEGER
! 138: * The number of columns in the arrays VL and/or VR. MM >= M.
! 139: *
! 140: * M (output) INTEGER
! 141: * The number of columns in the arrays VL and/or VR actually
! 142: * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
! 143: * is set to N. Each selected real eigenvector occupies one
! 144: * column and each selected complex eigenvector occupies two
! 145: * columns.
! 146: *
! 147: * WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
! 148: *
! 149: * INFO (output) INTEGER
! 150: * = 0: successful exit.
! 151: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 152: * > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex
! 153: * eigenvalue.
! 154: *
! 155: * Further Details
! 156: * ===============
! 157: *
! 158: * Allocation of workspace:
! 159: * ---------- -- ---------
! 160: *
! 161: * WORK( j ) = 1-norm of j-th column of A, above the diagonal
! 162: * WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
! 163: * WORK( 2*N+1:3*N ) = real part of eigenvector
! 164: * WORK( 3*N+1:4*N ) = imaginary part of eigenvector
! 165: * WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
! 166: * WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
! 167: *
! 168: * Rowwise vs. columnwise solution methods:
! 169: * ------- -- ---------- -------- -------
! 170: *
! 171: * Finding a generalized eigenvector consists basically of solving the
! 172: * singular triangular system
! 173: *
! 174: * (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
! 175: *
! 176: * Consider finding the i-th right eigenvector (assume all eigenvalues
! 177: * are real). The equation to be solved is:
! 178: * n i
! 179: * 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
! 180: * k=j k=j
! 181: *
! 182: * where C = (A - w B) (The components v(i+1:n) are 0.)
! 183: *
! 184: * The "rowwise" method is:
! 185: *
! 186: * (1) v(i) := 1
! 187: * for j = i-1,. . .,1:
! 188: * i
! 189: * (2) compute s = - sum C(j,k) v(k) and
! 190: * k=j+1
! 191: *
! 192: * (3) v(j) := s / C(j,j)
! 193: *
! 194: * Step 2 is sometimes called the "dot product" step, since it is an
! 195: * inner product between the j-th row and the portion of the eigenvector
! 196: * that has been computed so far.
! 197: *
! 198: * The "columnwise" method consists basically in doing the sums
! 199: * for all the rows in parallel. As each v(j) is computed, the
! 200: * contribution of v(j) times the j-th column of C is added to the
! 201: * partial sums. Since FORTRAN arrays are stored columnwise, this has
! 202: * the advantage that at each step, the elements of C that are accessed
! 203: * are adjacent to one another, whereas with the rowwise method, the
! 204: * elements accessed at a step are spaced LDS (and LDP) words apart.
! 205: *
! 206: * When finding left eigenvectors, the matrix in question is the
! 207: * transpose of the one in storage, so the rowwise method then
! 208: * actually accesses columns of A and B at each step, and so is the
! 209: * preferred method.
! 210: *
! 211: * =====================================================================
! 212: *
! 213: * .. Parameters ..
! 214: DOUBLE PRECISION ZERO, ONE, SAFETY
! 215: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
! 216: $ SAFETY = 1.0D+2 )
! 217: * ..
! 218: * .. Local Scalars ..
! 219: LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
! 220: $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB
! 221: INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
! 222: $ J, JA, JC, JE, JR, JW, NA, NW
! 223: DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
! 224: $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
! 225: $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
! 226: $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
! 227: $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
! 228: $ XSCALE
! 229: * ..
! 230: * .. Local Arrays ..
! 231: DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
! 232: $ SUMP( 2, 2 )
! 233: * ..
! 234: * .. External Functions ..
! 235: LOGICAL LSAME
! 236: DOUBLE PRECISION DLAMCH
! 237: EXTERNAL LSAME, DLAMCH
! 238: * ..
! 239: * .. External Subroutines ..
! 240: EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA
! 241: * ..
! 242: * .. Intrinsic Functions ..
! 243: INTRINSIC ABS, MAX, MIN
! 244: * ..
! 245: * .. Executable Statements ..
! 246: *
! 247: * Decode and Test the input parameters
! 248: *
! 249: IF( LSAME( HOWMNY, 'A' ) ) THEN
! 250: IHWMNY = 1
! 251: ILALL = .TRUE.
! 252: ILBACK = .FALSE.
! 253: ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
! 254: IHWMNY = 2
! 255: ILALL = .FALSE.
! 256: ILBACK = .FALSE.
! 257: ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
! 258: IHWMNY = 3
! 259: ILALL = .TRUE.
! 260: ILBACK = .TRUE.
! 261: ELSE
! 262: IHWMNY = -1
! 263: ILALL = .TRUE.
! 264: END IF
! 265: *
! 266: IF( LSAME( SIDE, 'R' ) ) THEN
! 267: ISIDE = 1
! 268: COMPL = .FALSE.
! 269: COMPR = .TRUE.
! 270: ELSE IF( LSAME( SIDE, 'L' ) ) THEN
! 271: ISIDE = 2
! 272: COMPL = .TRUE.
! 273: COMPR = .FALSE.
! 274: ELSE IF( LSAME( SIDE, 'B' ) ) THEN
! 275: ISIDE = 3
! 276: COMPL = .TRUE.
! 277: COMPR = .TRUE.
! 278: ELSE
! 279: ISIDE = -1
! 280: END IF
! 281: *
! 282: INFO = 0
! 283: IF( ISIDE.LT.0 ) THEN
! 284: INFO = -1
! 285: ELSE IF( IHWMNY.LT.0 ) THEN
! 286: INFO = -2
! 287: ELSE IF( N.LT.0 ) THEN
! 288: INFO = -4
! 289: ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
! 290: INFO = -6
! 291: ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
! 292: INFO = -8
! 293: END IF
! 294: IF( INFO.NE.0 ) THEN
! 295: CALL XERBLA( 'DTGEVC', -INFO )
! 296: RETURN
! 297: END IF
! 298: *
! 299: * Count the number of eigenvectors to be computed
! 300: *
! 301: IF( .NOT.ILALL ) THEN
! 302: IM = 0
! 303: ILCPLX = .FALSE.
! 304: DO 10 J = 1, N
! 305: IF( ILCPLX ) THEN
! 306: ILCPLX = .FALSE.
! 307: GO TO 10
! 308: END IF
! 309: IF( J.LT.N ) THEN
! 310: IF( S( J+1, J ).NE.ZERO )
! 311: $ ILCPLX = .TRUE.
! 312: END IF
! 313: IF( ILCPLX ) THEN
! 314: IF( SELECT( J ) .OR. SELECT( J+1 ) )
! 315: $ IM = IM + 2
! 316: ELSE
! 317: IF( SELECT( J ) )
! 318: $ IM = IM + 1
! 319: END IF
! 320: 10 CONTINUE
! 321: ELSE
! 322: IM = N
! 323: END IF
! 324: *
! 325: * Check 2-by-2 diagonal blocks of A, B
! 326: *
! 327: ILABAD = .FALSE.
! 328: ILBBAD = .FALSE.
! 329: DO 20 J = 1, N - 1
! 330: IF( S( J+1, J ).NE.ZERO ) THEN
! 331: IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
! 332: $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
! 333: IF( J.LT.N-1 ) THEN
! 334: IF( S( J+2, J+1 ).NE.ZERO )
! 335: $ ILABAD = .TRUE.
! 336: END IF
! 337: END IF
! 338: 20 CONTINUE
! 339: *
! 340: IF( ILABAD ) THEN
! 341: INFO = -5
! 342: ELSE IF( ILBBAD ) THEN
! 343: INFO = -7
! 344: ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
! 345: INFO = -10
! 346: ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
! 347: INFO = -12
! 348: ELSE IF( MM.LT.IM ) THEN
! 349: INFO = -13
! 350: END IF
! 351: IF( INFO.NE.0 ) THEN
! 352: CALL XERBLA( 'DTGEVC', -INFO )
! 353: RETURN
! 354: END IF
! 355: *
! 356: * Quick return if possible
! 357: *
! 358: M = IM
! 359: IF( N.EQ.0 )
! 360: $ RETURN
! 361: *
! 362: * Machine Constants
! 363: *
! 364: SAFMIN = DLAMCH( 'Safe minimum' )
! 365: BIG = ONE / SAFMIN
! 366: CALL DLABAD( SAFMIN, BIG )
! 367: ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
! 368: SMALL = SAFMIN*N / ULP
! 369: BIG = ONE / SMALL
! 370: BIGNUM = ONE / ( SAFMIN*N )
! 371: *
! 372: * Compute the 1-norm of each column of the strictly upper triangular
! 373: * part (i.e., excluding all elements belonging to the diagonal
! 374: * blocks) of A and B to check for possible overflow in the
! 375: * triangular solver.
! 376: *
! 377: ANORM = ABS( S( 1, 1 ) )
! 378: IF( N.GT.1 )
! 379: $ ANORM = ANORM + ABS( S( 2, 1 ) )
! 380: BNORM = ABS( P( 1, 1 ) )
! 381: WORK( 1 ) = ZERO
! 382: WORK( N+1 ) = ZERO
! 383: *
! 384: DO 50 J = 2, N
! 385: TEMP = ZERO
! 386: TEMP2 = ZERO
! 387: IF( S( J, J-1 ).EQ.ZERO ) THEN
! 388: IEND = J - 1
! 389: ELSE
! 390: IEND = J - 2
! 391: END IF
! 392: DO 30 I = 1, IEND
! 393: TEMP = TEMP + ABS( S( I, J ) )
! 394: TEMP2 = TEMP2 + ABS( P( I, J ) )
! 395: 30 CONTINUE
! 396: WORK( J ) = TEMP
! 397: WORK( N+J ) = TEMP2
! 398: DO 40 I = IEND + 1, MIN( J+1, N )
! 399: TEMP = TEMP + ABS( S( I, J ) )
! 400: TEMP2 = TEMP2 + ABS( P( I, J ) )
! 401: 40 CONTINUE
! 402: ANORM = MAX( ANORM, TEMP )
! 403: BNORM = MAX( BNORM, TEMP2 )
! 404: 50 CONTINUE
! 405: *
! 406: ASCALE = ONE / MAX( ANORM, SAFMIN )
! 407: BSCALE = ONE / MAX( BNORM, SAFMIN )
! 408: *
! 409: * Left eigenvectors
! 410: *
! 411: IF( COMPL ) THEN
! 412: IEIG = 0
! 413: *
! 414: * Main loop over eigenvalues
! 415: *
! 416: ILCPLX = .FALSE.
! 417: DO 220 JE = 1, N
! 418: *
! 419: * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
! 420: * (b) this would be the second of a complex pair.
! 421: * Check for complex eigenvalue, so as to be sure of which
! 422: * entry(-ies) of SELECT to look at.
! 423: *
! 424: IF( ILCPLX ) THEN
! 425: ILCPLX = .FALSE.
! 426: GO TO 220
! 427: END IF
! 428: NW = 1
! 429: IF( JE.LT.N ) THEN
! 430: IF( S( JE+1, JE ).NE.ZERO ) THEN
! 431: ILCPLX = .TRUE.
! 432: NW = 2
! 433: END IF
! 434: END IF
! 435: IF( ILALL ) THEN
! 436: ILCOMP = .TRUE.
! 437: ELSE IF( ILCPLX ) THEN
! 438: ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
! 439: ELSE
! 440: ILCOMP = SELECT( JE )
! 441: END IF
! 442: IF( .NOT.ILCOMP )
! 443: $ GO TO 220
! 444: *
! 445: * Decide if (a) singular pencil, (b) real eigenvalue, or
! 446: * (c) complex eigenvalue.
! 447: *
! 448: IF( .NOT.ILCPLX ) THEN
! 449: IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
! 450: $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
! 451: *
! 452: * Singular matrix pencil -- return unit eigenvector
! 453: *
! 454: IEIG = IEIG + 1
! 455: DO 60 JR = 1, N
! 456: VL( JR, IEIG ) = ZERO
! 457: 60 CONTINUE
! 458: VL( IEIG, IEIG ) = ONE
! 459: GO TO 220
! 460: END IF
! 461: END IF
! 462: *
! 463: * Clear vector
! 464: *
! 465: DO 70 JR = 1, NW*N
! 466: WORK( 2*N+JR ) = ZERO
! 467: 70 CONTINUE
! 468: * T
! 469: * Compute coefficients in ( a A - b B ) y = 0
! 470: * a is ACOEF
! 471: * b is BCOEFR + i*BCOEFI
! 472: *
! 473: IF( .NOT.ILCPLX ) THEN
! 474: *
! 475: * Real eigenvalue
! 476: *
! 477: TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
! 478: $ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
! 479: SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
! 480: SBETA = ( TEMP*P( JE, JE ) )*BSCALE
! 481: ACOEF = SBETA*ASCALE
! 482: BCOEFR = SALFAR*BSCALE
! 483: BCOEFI = ZERO
! 484: *
! 485: * Scale to avoid underflow
! 486: *
! 487: SCALE = ONE
! 488: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
! 489: LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
! 490: $ SMALL
! 491: IF( LSA )
! 492: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
! 493: IF( LSB )
! 494: $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
! 495: $ MIN( BNORM, BIG ) )
! 496: IF( LSA .OR. LSB ) THEN
! 497: SCALE = MIN( SCALE, ONE /
! 498: $ ( SAFMIN*MAX( ONE, ABS( ACOEF ),
! 499: $ ABS( BCOEFR ) ) ) )
! 500: IF( LSA ) THEN
! 501: ACOEF = ASCALE*( SCALE*SBETA )
! 502: ELSE
! 503: ACOEF = SCALE*ACOEF
! 504: END IF
! 505: IF( LSB ) THEN
! 506: BCOEFR = BSCALE*( SCALE*SALFAR )
! 507: ELSE
! 508: BCOEFR = SCALE*BCOEFR
! 509: END IF
! 510: END IF
! 511: ACOEFA = ABS( ACOEF )
! 512: BCOEFA = ABS( BCOEFR )
! 513: *
! 514: * First component is 1
! 515: *
! 516: WORK( 2*N+JE ) = ONE
! 517: XMAX = ONE
! 518: ELSE
! 519: *
! 520: * Complex eigenvalue
! 521: *
! 522: CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
! 523: $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
! 524: $ BCOEFI )
! 525: BCOEFI = -BCOEFI
! 526: IF( BCOEFI.EQ.ZERO ) THEN
! 527: INFO = JE
! 528: RETURN
! 529: END IF
! 530: *
! 531: * Scale to avoid over/underflow
! 532: *
! 533: ACOEFA = ABS( ACOEF )
! 534: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
! 535: SCALE = ONE
! 536: IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
! 537: $ SCALE = ( SAFMIN / ULP ) / ACOEFA
! 538: IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
! 539: $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
! 540: IF( SAFMIN*ACOEFA.GT.ASCALE )
! 541: $ SCALE = ASCALE / ( SAFMIN*ACOEFA )
! 542: IF( SAFMIN*BCOEFA.GT.BSCALE )
! 543: $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
! 544: IF( SCALE.NE.ONE ) THEN
! 545: ACOEF = SCALE*ACOEF
! 546: ACOEFA = ABS( ACOEF )
! 547: BCOEFR = SCALE*BCOEFR
! 548: BCOEFI = SCALE*BCOEFI
! 549: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
! 550: END IF
! 551: *
! 552: * Compute first two components of eigenvector
! 553: *
! 554: TEMP = ACOEF*S( JE+1, JE )
! 555: TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
! 556: TEMP2I = -BCOEFI*P( JE, JE )
! 557: IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
! 558: WORK( 2*N+JE ) = ONE
! 559: WORK( 3*N+JE ) = ZERO
! 560: WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
! 561: WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
! 562: ELSE
! 563: WORK( 2*N+JE+1 ) = ONE
! 564: WORK( 3*N+JE+1 ) = ZERO
! 565: TEMP = ACOEF*S( JE, JE+1 )
! 566: WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
! 567: $ S( JE+1, JE+1 ) ) / TEMP
! 568: WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
! 569: END IF
! 570: XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
! 571: $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
! 572: END IF
! 573: *
! 574: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
! 575: *
! 576: * T
! 577: * Triangular solve of (a A - b B) y = 0
! 578: *
! 579: * T
! 580: * (rowwise in (a A - b B) , or columnwise in (a A - b B) )
! 581: *
! 582: IL2BY2 = .FALSE.
! 583: *
! 584: DO 160 J = JE + NW, N
! 585: IF( IL2BY2 ) THEN
! 586: IL2BY2 = .FALSE.
! 587: GO TO 160
! 588: END IF
! 589: *
! 590: NA = 1
! 591: BDIAG( 1 ) = P( J, J )
! 592: IF( J.LT.N ) THEN
! 593: IF( S( J+1, J ).NE.ZERO ) THEN
! 594: IL2BY2 = .TRUE.
! 595: BDIAG( 2 ) = P( J+1, J+1 )
! 596: NA = 2
! 597: END IF
! 598: END IF
! 599: *
! 600: * Check whether scaling is necessary for dot products
! 601: *
! 602: XSCALE = ONE / MAX( ONE, XMAX )
! 603: TEMP = MAX( WORK( J ), WORK( N+J ),
! 604: $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
! 605: IF( IL2BY2 )
! 606: $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
! 607: $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
! 608: IF( TEMP.GT.BIGNUM*XSCALE ) THEN
! 609: DO 90 JW = 0, NW - 1
! 610: DO 80 JR = JE, J - 1
! 611: WORK( ( JW+2 )*N+JR ) = XSCALE*
! 612: $ WORK( ( JW+2 )*N+JR )
! 613: 80 CONTINUE
! 614: 90 CONTINUE
! 615: XMAX = XMAX*XSCALE
! 616: END IF
! 617: *
! 618: * Compute dot products
! 619: *
! 620: * j-1
! 621: * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
! 622: * k=je
! 623: *
! 624: * To reduce the op count, this is done as
! 625: *
! 626: * _ j-1 _ j-1
! 627: * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
! 628: * k=je k=je
! 629: *
! 630: * which may cause underflow problems if A or B are close
! 631: * to underflow. (E.g., less than SMALL.)
! 632: *
! 633: *
! 634: * A series of compiler directives to defeat vectorization
! 635: * for the next loop
! 636: *
! 637: *$PL$ CMCHAR=' '
! 638: CDIR$ NEXTSCALAR
! 639: C$DIR SCALAR
! 640: CDIR$ NEXT SCALAR
! 641: CVD$L NOVECTOR
! 642: CDEC$ NOVECTOR
! 643: CVD$ NOVECTOR
! 644: *VDIR NOVECTOR
! 645: *VOCL LOOP,SCALAR
! 646: CIBM PREFER SCALAR
! 647: *$PL$ CMCHAR='*'
! 648: *
! 649: DO 120 JW = 1, NW
! 650: *
! 651: *$PL$ CMCHAR=' '
! 652: CDIR$ NEXTSCALAR
! 653: C$DIR SCALAR
! 654: CDIR$ NEXT SCALAR
! 655: CVD$L NOVECTOR
! 656: CDEC$ NOVECTOR
! 657: CVD$ NOVECTOR
! 658: *VDIR NOVECTOR
! 659: *VOCL LOOP,SCALAR
! 660: CIBM PREFER SCALAR
! 661: *$PL$ CMCHAR='*'
! 662: *
! 663: DO 110 JA = 1, NA
! 664: SUMS( JA, JW ) = ZERO
! 665: SUMP( JA, JW ) = ZERO
! 666: *
! 667: DO 100 JR = JE, J - 1
! 668: SUMS( JA, JW ) = SUMS( JA, JW ) +
! 669: $ S( JR, J+JA-1 )*
! 670: $ WORK( ( JW+1 )*N+JR )
! 671: SUMP( JA, JW ) = SUMP( JA, JW ) +
! 672: $ P( JR, J+JA-1 )*
! 673: $ WORK( ( JW+1 )*N+JR )
! 674: 100 CONTINUE
! 675: 110 CONTINUE
! 676: 120 CONTINUE
! 677: *
! 678: *$PL$ CMCHAR=' '
! 679: CDIR$ NEXTSCALAR
! 680: C$DIR SCALAR
! 681: CDIR$ NEXT SCALAR
! 682: CVD$L NOVECTOR
! 683: CDEC$ NOVECTOR
! 684: CVD$ NOVECTOR
! 685: *VDIR NOVECTOR
! 686: *VOCL LOOP,SCALAR
! 687: CIBM PREFER SCALAR
! 688: *$PL$ CMCHAR='*'
! 689: *
! 690: DO 130 JA = 1, NA
! 691: IF( ILCPLX ) THEN
! 692: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
! 693: $ BCOEFR*SUMP( JA, 1 ) -
! 694: $ BCOEFI*SUMP( JA, 2 )
! 695: SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
! 696: $ BCOEFR*SUMP( JA, 2 ) +
! 697: $ BCOEFI*SUMP( JA, 1 )
! 698: ELSE
! 699: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
! 700: $ BCOEFR*SUMP( JA, 1 )
! 701: END IF
! 702: 130 CONTINUE
! 703: *
! 704: * T
! 705: * Solve ( a A - b B ) y = SUM(,)
! 706: * with scaling and perturbation of the denominator
! 707: *
! 708: CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
! 709: $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
! 710: $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
! 711: $ IINFO )
! 712: IF( SCALE.LT.ONE ) THEN
! 713: DO 150 JW = 0, NW - 1
! 714: DO 140 JR = JE, J - 1
! 715: WORK( ( JW+2 )*N+JR ) = SCALE*
! 716: $ WORK( ( JW+2 )*N+JR )
! 717: 140 CONTINUE
! 718: 150 CONTINUE
! 719: XMAX = SCALE*XMAX
! 720: END IF
! 721: XMAX = MAX( XMAX, TEMP )
! 722: 160 CONTINUE
! 723: *
! 724: * Copy eigenvector to VL, back transforming if
! 725: * HOWMNY='B'.
! 726: *
! 727: IEIG = IEIG + 1
! 728: IF( ILBACK ) THEN
! 729: DO 170 JW = 0, NW - 1
! 730: CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
! 731: $ WORK( ( JW+2 )*N+JE ), 1, ZERO,
! 732: $ WORK( ( JW+4 )*N+1 ), 1 )
! 733: 170 CONTINUE
! 734: CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
! 735: $ LDVL )
! 736: IBEG = 1
! 737: ELSE
! 738: CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
! 739: $ LDVL )
! 740: IBEG = JE
! 741: END IF
! 742: *
! 743: * Scale eigenvector
! 744: *
! 745: XMAX = ZERO
! 746: IF( ILCPLX ) THEN
! 747: DO 180 J = IBEG, N
! 748: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
! 749: $ ABS( VL( J, IEIG+1 ) ) )
! 750: 180 CONTINUE
! 751: ELSE
! 752: DO 190 J = IBEG, N
! 753: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
! 754: 190 CONTINUE
! 755: END IF
! 756: *
! 757: IF( XMAX.GT.SAFMIN ) THEN
! 758: XSCALE = ONE / XMAX
! 759: *
! 760: DO 210 JW = 0, NW - 1
! 761: DO 200 JR = IBEG, N
! 762: VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
! 763: 200 CONTINUE
! 764: 210 CONTINUE
! 765: END IF
! 766: IEIG = IEIG + NW - 1
! 767: *
! 768: 220 CONTINUE
! 769: END IF
! 770: *
! 771: * Right eigenvectors
! 772: *
! 773: IF( COMPR ) THEN
! 774: IEIG = IM + 1
! 775: *
! 776: * Main loop over eigenvalues
! 777: *
! 778: ILCPLX = .FALSE.
! 779: DO 500 JE = N, 1, -1
! 780: *
! 781: * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
! 782: * (b) this would be the second of a complex pair.
! 783: * Check for complex eigenvalue, so as to be sure of which
! 784: * entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
! 785: * or SELECT(JE-1).
! 786: * If this is a complex pair, the 2-by-2 diagonal block
! 787: * corresponding to the eigenvalue is in rows/columns JE-1:JE
! 788: *
! 789: IF( ILCPLX ) THEN
! 790: ILCPLX = .FALSE.
! 791: GO TO 500
! 792: END IF
! 793: NW = 1
! 794: IF( JE.GT.1 ) THEN
! 795: IF( S( JE, JE-1 ).NE.ZERO ) THEN
! 796: ILCPLX = .TRUE.
! 797: NW = 2
! 798: END IF
! 799: END IF
! 800: IF( ILALL ) THEN
! 801: ILCOMP = .TRUE.
! 802: ELSE IF( ILCPLX ) THEN
! 803: ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
! 804: ELSE
! 805: ILCOMP = SELECT( JE )
! 806: END IF
! 807: IF( .NOT.ILCOMP )
! 808: $ GO TO 500
! 809: *
! 810: * Decide if (a) singular pencil, (b) real eigenvalue, or
! 811: * (c) complex eigenvalue.
! 812: *
! 813: IF( .NOT.ILCPLX ) THEN
! 814: IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
! 815: $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
! 816: *
! 817: * Singular matrix pencil -- unit eigenvector
! 818: *
! 819: IEIG = IEIG - 1
! 820: DO 230 JR = 1, N
! 821: VR( JR, IEIG ) = ZERO
! 822: 230 CONTINUE
! 823: VR( IEIG, IEIG ) = ONE
! 824: GO TO 500
! 825: END IF
! 826: END IF
! 827: *
! 828: * Clear vector
! 829: *
! 830: DO 250 JW = 0, NW - 1
! 831: DO 240 JR = 1, N
! 832: WORK( ( JW+2 )*N+JR ) = ZERO
! 833: 240 CONTINUE
! 834: 250 CONTINUE
! 835: *
! 836: * Compute coefficients in ( a A - b B ) x = 0
! 837: * a is ACOEF
! 838: * b is BCOEFR + i*BCOEFI
! 839: *
! 840: IF( .NOT.ILCPLX ) THEN
! 841: *
! 842: * Real eigenvalue
! 843: *
! 844: TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
! 845: $ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
! 846: SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
! 847: SBETA = ( TEMP*P( JE, JE ) )*BSCALE
! 848: ACOEF = SBETA*ASCALE
! 849: BCOEFR = SALFAR*BSCALE
! 850: BCOEFI = ZERO
! 851: *
! 852: * Scale to avoid underflow
! 853: *
! 854: SCALE = ONE
! 855: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
! 856: LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
! 857: $ SMALL
! 858: IF( LSA )
! 859: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
! 860: IF( LSB )
! 861: $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
! 862: $ MIN( BNORM, BIG ) )
! 863: IF( LSA .OR. LSB ) THEN
! 864: SCALE = MIN( SCALE, ONE /
! 865: $ ( SAFMIN*MAX( ONE, ABS( ACOEF ),
! 866: $ ABS( BCOEFR ) ) ) )
! 867: IF( LSA ) THEN
! 868: ACOEF = ASCALE*( SCALE*SBETA )
! 869: ELSE
! 870: ACOEF = SCALE*ACOEF
! 871: END IF
! 872: IF( LSB ) THEN
! 873: BCOEFR = BSCALE*( SCALE*SALFAR )
! 874: ELSE
! 875: BCOEFR = SCALE*BCOEFR
! 876: END IF
! 877: END IF
! 878: ACOEFA = ABS( ACOEF )
! 879: BCOEFA = ABS( BCOEFR )
! 880: *
! 881: * First component is 1
! 882: *
! 883: WORK( 2*N+JE ) = ONE
! 884: XMAX = ONE
! 885: *
! 886: * Compute contribution from column JE of A and B to sum
! 887: * (See "Further Details", above.)
! 888: *
! 889: DO 260 JR = 1, JE - 1
! 890: WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
! 891: $ ACOEF*S( JR, JE )
! 892: 260 CONTINUE
! 893: ELSE
! 894: *
! 895: * Complex eigenvalue
! 896: *
! 897: CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
! 898: $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
! 899: $ BCOEFI )
! 900: IF( BCOEFI.EQ.ZERO ) THEN
! 901: INFO = JE - 1
! 902: RETURN
! 903: END IF
! 904: *
! 905: * Scale to avoid over/underflow
! 906: *
! 907: ACOEFA = ABS( ACOEF )
! 908: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
! 909: SCALE = ONE
! 910: IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
! 911: $ SCALE = ( SAFMIN / ULP ) / ACOEFA
! 912: IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
! 913: $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
! 914: IF( SAFMIN*ACOEFA.GT.ASCALE )
! 915: $ SCALE = ASCALE / ( SAFMIN*ACOEFA )
! 916: IF( SAFMIN*BCOEFA.GT.BSCALE )
! 917: $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
! 918: IF( SCALE.NE.ONE ) THEN
! 919: ACOEF = SCALE*ACOEF
! 920: ACOEFA = ABS( ACOEF )
! 921: BCOEFR = SCALE*BCOEFR
! 922: BCOEFI = SCALE*BCOEFI
! 923: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
! 924: END IF
! 925: *
! 926: * Compute first two components of eigenvector
! 927: * and contribution to sums
! 928: *
! 929: TEMP = ACOEF*S( JE, JE-1 )
! 930: TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
! 931: TEMP2I = -BCOEFI*P( JE, JE )
! 932: IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
! 933: WORK( 2*N+JE ) = ONE
! 934: WORK( 3*N+JE ) = ZERO
! 935: WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
! 936: WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
! 937: ELSE
! 938: WORK( 2*N+JE-1 ) = ONE
! 939: WORK( 3*N+JE-1 ) = ZERO
! 940: TEMP = ACOEF*S( JE-1, JE )
! 941: WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
! 942: $ S( JE-1, JE-1 ) ) / TEMP
! 943: WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
! 944: END IF
! 945: *
! 946: XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
! 947: $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
! 948: *
! 949: * Compute contribution from columns JE and JE-1
! 950: * of A and B to the sums.
! 951: *
! 952: CREALA = ACOEF*WORK( 2*N+JE-1 )
! 953: CIMAGA = ACOEF*WORK( 3*N+JE-1 )
! 954: CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
! 955: $ BCOEFI*WORK( 3*N+JE-1 )
! 956: CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
! 957: $ BCOEFR*WORK( 3*N+JE-1 )
! 958: CRE2A = ACOEF*WORK( 2*N+JE )
! 959: CIM2A = ACOEF*WORK( 3*N+JE )
! 960: CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
! 961: CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
! 962: DO 270 JR = 1, JE - 2
! 963: WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
! 964: $ CREALB*P( JR, JE-1 ) -
! 965: $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
! 966: WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
! 967: $ CIMAGB*P( JR, JE-1 ) -
! 968: $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
! 969: 270 CONTINUE
! 970: END IF
! 971: *
! 972: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
! 973: *
! 974: * Columnwise triangular solve of (a A - b B) x = 0
! 975: *
! 976: IL2BY2 = .FALSE.
! 977: DO 370 J = JE - NW, 1, -1
! 978: *
! 979: * If a 2-by-2 block, is in position j-1:j, wait until
! 980: * next iteration to process it (when it will be j:j+1)
! 981: *
! 982: IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
! 983: IF( S( J, J-1 ).NE.ZERO ) THEN
! 984: IL2BY2 = .TRUE.
! 985: GO TO 370
! 986: END IF
! 987: END IF
! 988: BDIAG( 1 ) = P( J, J )
! 989: IF( IL2BY2 ) THEN
! 990: NA = 2
! 991: BDIAG( 2 ) = P( J+1, J+1 )
! 992: ELSE
! 993: NA = 1
! 994: END IF
! 995: *
! 996: * Compute x(j) (and x(j+1), if 2-by-2 block)
! 997: *
! 998: CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
! 999: $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
! 1000: $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
! 1001: $ IINFO )
! 1002: IF( SCALE.LT.ONE ) THEN
! 1003: *
! 1004: DO 290 JW = 0, NW - 1
! 1005: DO 280 JR = 1, JE
! 1006: WORK( ( JW+2 )*N+JR ) = SCALE*
! 1007: $ WORK( ( JW+2 )*N+JR )
! 1008: 280 CONTINUE
! 1009: 290 CONTINUE
! 1010: END IF
! 1011: XMAX = MAX( SCALE*XMAX, TEMP )
! 1012: *
! 1013: DO 310 JW = 1, NW
! 1014: DO 300 JA = 1, NA
! 1015: WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
! 1016: 300 CONTINUE
! 1017: 310 CONTINUE
! 1018: *
! 1019: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
! 1020: *
! 1021: IF( J.GT.1 ) THEN
! 1022: *
! 1023: * Check whether scaling is necessary for sum.
! 1024: *
! 1025: XSCALE = ONE / MAX( ONE, XMAX )
! 1026: TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
! 1027: IF( IL2BY2 )
! 1028: $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
! 1029: $ WORK( N+J+1 ) )
! 1030: TEMP = MAX( TEMP, ACOEFA, BCOEFA )
! 1031: IF( TEMP.GT.BIGNUM*XSCALE ) THEN
! 1032: *
! 1033: DO 330 JW = 0, NW - 1
! 1034: DO 320 JR = 1, JE
! 1035: WORK( ( JW+2 )*N+JR ) = XSCALE*
! 1036: $ WORK( ( JW+2 )*N+JR )
! 1037: 320 CONTINUE
! 1038: 330 CONTINUE
! 1039: XMAX = XMAX*XSCALE
! 1040: END IF
! 1041: *
! 1042: * Compute the contributions of the off-diagonals of
! 1043: * column j (and j+1, if 2-by-2 block) of A and B to the
! 1044: * sums.
! 1045: *
! 1046: *
! 1047: DO 360 JA = 1, NA
! 1048: IF( ILCPLX ) THEN
! 1049: CREALA = ACOEF*WORK( 2*N+J+JA-1 )
! 1050: CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
! 1051: CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
! 1052: $ BCOEFI*WORK( 3*N+J+JA-1 )
! 1053: CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
! 1054: $ BCOEFR*WORK( 3*N+J+JA-1 )
! 1055: DO 340 JR = 1, J - 1
! 1056: WORK( 2*N+JR ) = WORK( 2*N+JR ) -
! 1057: $ CREALA*S( JR, J+JA-1 ) +
! 1058: $ CREALB*P( JR, J+JA-1 )
! 1059: WORK( 3*N+JR ) = WORK( 3*N+JR ) -
! 1060: $ CIMAGA*S( JR, J+JA-1 ) +
! 1061: $ CIMAGB*P( JR, J+JA-1 )
! 1062: 340 CONTINUE
! 1063: ELSE
! 1064: CREALA = ACOEF*WORK( 2*N+J+JA-1 )
! 1065: CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
! 1066: DO 350 JR = 1, J - 1
! 1067: WORK( 2*N+JR ) = WORK( 2*N+JR ) -
! 1068: $ CREALA*S( JR, J+JA-1 ) +
! 1069: $ CREALB*P( JR, J+JA-1 )
! 1070: 350 CONTINUE
! 1071: END IF
! 1072: 360 CONTINUE
! 1073: END IF
! 1074: *
! 1075: IL2BY2 = .FALSE.
! 1076: 370 CONTINUE
! 1077: *
! 1078: * Copy eigenvector to VR, back transforming if
! 1079: * HOWMNY='B'.
! 1080: *
! 1081: IEIG = IEIG - NW
! 1082: IF( ILBACK ) THEN
! 1083: *
! 1084: DO 410 JW = 0, NW - 1
! 1085: DO 380 JR = 1, N
! 1086: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
! 1087: $ VR( JR, 1 )
! 1088: 380 CONTINUE
! 1089: *
! 1090: * A series of compiler directives to defeat
! 1091: * vectorization for the next loop
! 1092: *
! 1093: *
! 1094: DO 400 JC = 2, JE
! 1095: DO 390 JR = 1, N
! 1096: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
! 1097: $ WORK( ( JW+2 )*N+JC )*VR( JR, JC )
! 1098: 390 CONTINUE
! 1099: 400 CONTINUE
! 1100: 410 CONTINUE
! 1101: *
! 1102: DO 430 JW = 0, NW - 1
! 1103: DO 420 JR = 1, N
! 1104: VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
! 1105: 420 CONTINUE
! 1106: 430 CONTINUE
! 1107: *
! 1108: IEND = N
! 1109: ELSE
! 1110: DO 450 JW = 0, NW - 1
! 1111: DO 440 JR = 1, N
! 1112: VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
! 1113: 440 CONTINUE
! 1114: 450 CONTINUE
! 1115: *
! 1116: IEND = JE
! 1117: END IF
! 1118: *
! 1119: * Scale eigenvector
! 1120: *
! 1121: XMAX = ZERO
! 1122: IF( ILCPLX ) THEN
! 1123: DO 460 J = 1, IEND
! 1124: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
! 1125: $ ABS( VR( J, IEIG+1 ) ) )
! 1126: 460 CONTINUE
! 1127: ELSE
! 1128: DO 470 J = 1, IEND
! 1129: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
! 1130: 470 CONTINUE
! 1131: END IF
! 1132: *
! 1133: IF( XMAX.GT.SAFMIN ) THEN
! 1134: XSCALE = ONE / XMAX
! 1135: DO 490 JW = 0, NW - 1
! 1136: DO 480 JR = 1, IEND
! 1137: VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
! 1138: 480 CONTINUE
! 1139: 490 CONTINUE
! 1140: END IF
! 1141: 500 CONTINUE
! 1142: END IF
! 1143: *
! 1144: RETURN
! 1145: *
! 1146: * End of DTGEVC
! 1147: *
! 1148: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgevc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack