Annotation of rpl/lapack/lapack/ztgevc.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
! 2: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, 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 RWORK( * )
! 16: COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
! 17: $ VR( LDVR, * ), WORK( * )
! 18: * ..
! 19: *
! 20: *
! 21: * Purpose
! 22: * =======
! 23: *
! 24: * ZTGEVC computes some or all of the right and/or left eigenvectors of
! 25: * a pair of complex matrices (S,P), where S and P are upper triangular.
! 26: * Matrix pairs of this type are produced by the generalized Schur
! 27: * factorization of a complex matrix pair (A,B):
! 28: *
! 29: * A = Q*S*Z**H, B = Q*P*Z**H
! 30: *
! 31: * as computed by ZGGHRD + ZHGEQZ.
! 32: *
! 33: * The right eigenvector x and the left eigenvector y of (S,P)
! 34: * corresponding to an eigenvalue w are defined by:
! 35: *
! 36: * S*x = w*P*x, (y**H)*S = w*(y**H)*P,
! 37: *
! 38: * where y**H denotes the conjugate tranpose of y.
! 39: * The eigenvalues are not input to this routine, but are computed
! 40: * directly from the diagonal elements of S and P.
! 41: *
! 42: * This routine returns the matrices X and/or Y of right and left
! 43: * eigenvectors of (S,P), or the products Z*X and/or Q*Y,
! 44: * where Z and Q are input matrices.
! 45: * If Q and Z are the unitary factors from the generalized Schur
! 46: * factorization of a matrix pair (A,B), then Z*X and Q*Y
! 47: * are the matrices of right and left eigenvectors of (A,B).
! 48: *
! 49: * Arguments
! 50: * =========
! 51: *
! 52: * SIDE (input) CHARACTER*1
! 53: * = 'R': compute right eigenvectors only;
! 54: * = 'L': compute left eigenvectors only;
! 55: * = 'B': compute both right and left eigenvectors.
! 56: *
! 57: * HOWMNY (input) CHARACTER*1
! 58: * = 'A': compute all right and/or left eigenvectors;
! 59: * = 'B': compute all right and/or left eigenvectors,
! 60: * backtransformed by the matrices in VR and/or VL;
! 61: * = 'S': compute selected right and/or left eigenvectors,
! 62: * specified by the logical array SELECT.
! 63: *
! 64: * SELECT (input) LOGICAL array, dimension (N)
! 65: * If HOWMNY='S', SELECT specifies the eigenvectors to be
! 66: * computed. The eigenvector corresponding to the j-th
! 67: * eigenvalue is computed if SELECT(j) = .TRUE..
! 68: * Not referenced if HOWMNY = 'A' or 'B'.
! 69: *
! 70: * N (input) INTEGER
! 71: * The order of the matrices S and P. N >= 0.
! 72: *
! 73: * S (input) COMPLEX*16 array, dimension (LDS,N)
! 74: * The upper triangular matrix S from a generalized Schur
! 75: * factorization, as computed by ZHGEQZ.
! 76: *
! 77: * LDS (input) INTEGER
! 78: * The leading dimension of array S. LDS >= max(1,N).
! 79: *
! 80: * P (input) COMPLEX*16 array, dimension (LDP,N)
! 81: * The upper triangular matrix P from a generalized Schur
! 82: * factorization, as computed by ZHGEQZ. P must have real
! 83: * diagonal elements.
! 84: *
! 85: * LDP (input) INTEGER
! 86: * The leading dimension of array P. LDP >= max(1,N).
! 87: *
! 88: * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM)
! 89: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
! 90: * contain an N-by-N matrix Q (usually the unitary matrix Q
! 91: * of left Schur vectors returned by ZHGEQZ).
! 92: * On exit, if SIDE = 'L' or 'B', VL contains:
! 93: * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
! 94: * if HOWMNY = 'B', the matrix Q*Y;
! 95: * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
! 96: * SELECT, stored consecutively in the columns of
! 97: * VL, in the same order as their eigenvalues.
! 98: * Not referenced if SIDE = 'R'.
! 99: *
! 100: * LDVL (input) INTEGER
! 101: * The leading dimension of array VL. LDVL >= 1, and if
! 102: * SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
! 103: *
! 104: * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM)
! 105: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
! 106: * contain an N-by-N matrix Q (usually the unitary matrix Z
! 107: * of right Schur vectors returned by ZHGEQZ).
! 108: * On exit, if SIDE = 'R' or 'B', VR contains:
! 109: * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
! 110: * if HOWMNY = 'B', the matrix Z*X;
! 111: * if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
! 112: * SELECT, stored consecutively in the columns of
! 113: * VR, in the same order as their eigenvalues.
! 114: * Not referenced if SIDE = 'L'.
! 115: *
! 116: * LDVR (input) INTEGER
! 117: * The leading dimension of the array VR. LDVR >= 1, and if
! 118: * SIDE = 'R' or 'B', LDVR >= N.
! 119: *
! 120: * MM (input) INTEGER
! 121: * The number of columns in the arrays VL and/or VR. MM >= M.
! 122: *
! 123: * M (output) INTEGER
! 124: * The number of columns in the arrays VL and/or VR actually
! 125: * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
! 126: * is set to N. Each selected eigenvector occupies one column.
! 127: *
! 128: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
! 129: *
! 130: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 131: *
! 132: * INFO (output) INTEGER
! 133: * = 0: successful exit.
! 134: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 135: *
! 136: * =====================================================================
! 137: *
! 138: * .. Parameters ..
! 139: DOUBLE PRECISION ZERO, ONE
! 140: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 141: COMPLEX*16 CZERO, CONE
! 142: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 143: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 144: * ..
! 145: * .. Local Scalars ..
! 146: LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
! 147: $ LSA, LSB
! 148: INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
! 149: $ J, JE, JR
! 150: DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
! 151: $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
! 152: $ SCALE, SMALL, TEMP, ULP, XMAX
! 153: COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
! 154: * ..
! 155: * .. External Functions ..
! 156: LOGICAL LSAME
! 157: DOUBLE PRECISION DLAMCH
! 158: COMPLEX*16 ZLADIV
! 159: EXTERNAL LSAME, DLAMCH, ZLADIV
! 160: * ..
! 161: * .. External Subroutines ..
! 162: EXTERNAL DLABAD, XERBLA, ZGEMV
! 163: * ..
! 164: * .. Intrinsic Functions ..
! 165: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
! 166: * ..
! 167: * .. Statement Functions ..
! 168: DOUBLE PRECISION ABS1
! 169: * ..
! 170: * .. Statement Function definitions ..
! 171: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
! 172: * ..
! 173: * .. Executable Statements ..
! 174: *
! 175: * Decode and Test the input parameters
! 176: *
! 177: IF( LSAME( HOWMNY, 'A' ) ) THEN
! 178: IHWMNY = 1
! 179: ILALL = .TRUE.
! 180: ILBACK = .FALSE.
! 181: ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
! 182: IHWMNY = 2
! 183: ILALL = .FALSE.
! 184: ILBACK = .FALSE.
! 185: ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
! 186: IHWMNY = 3
! 187: ILALL = .TRUE.
! 188: ILBACK = .TRUE.
! 189: ELSE
! 190: IHWMNY = -1
! 191: END IF
! 192: *
! 193: IF( LSAME( SIDE, 'R' ) ) THEN
! 194: ISIDE = 1
! 195: COMPL = .FALSE.
! 196: COMPR = .TRUE.
! 197: ELSE IF( LSAME( SIDE, 'L' ) ) THEN
! 198: ISIDE = 2
! 199: COMPL = .TRUE.
! 200: COMPR = .FALSE.
! 201: ELSE IF( LSAME( SIDE, 'B' ) ) THEN
! 202: ISIDE = 3
! 203: COMPL = .TRUE.
! 204: COMPR = .TRUE.
! 205: ELSE
! 206: ISIDE = -1
! 207: END IF
! 208: *
! 209: INFO = 0
! 210: IF( ISIDE.LT.0 ) THEN
! 211: INFO = -1
! 212: ELSE IF( IHWMNY.LT.0 ) THEN
! 213: INFO = -2
! 214: ELSE IF( N.LT.0 ) THEN
! 215: INFO = -4
! 216: ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
! 217: INFO = -6
! 218: ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
! 219: INFO = -8
! 220: END IF
! 221: IF( INFO.NE.0 ) THEN
! 222: CALL XERBLA( 'ZTGEVC', -INFO )
! 223: RETURN
! 224: END IF
! 225: *
! 226: * Count the number of eigenvectors
! 227: *
! 228: IF( .NOT.ILALL ) THEN
! 229: IM = 0
! 230: DO 10 J = 1, N
! 231: IF( SELECT( J ) )
! 232: $ IM = IM + 1
! 233: 10 CONTINUE
! 234: ELSE
! 235: IM = N
! 236: END IF
! 237: *
! 238: * Check diagonal of B
! 239: *
! 240: ILBBAD = .FALSE.
! 241: DO 20 J = 1, N
! 242: IF( DIMAG( P( J, J ) ).NE.ZERO )
! 243: $ ILBBAD = .TRUE.
! 244: 20 CONTINUE
! 245: *
! 246: IF( ILBBAD ) THEN
! 247: INFO = -7
! 248: ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
! 249: INFO = -10
! 250: ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
! 251: INFO = -12
! 252: ELSE IF( MM.LT.IM ) THEN
! 253: INFO = -13
! 254: END IF
! 255: IF( INFO.NE.0 ) THEN
! 256: CALL XERBLA( 'ZTGEVC', -INFO )
! 257: RETURN
! 258: END IF
! 259: *
! 260: * Quick return if possible
! 261: *
! 262: M = IM
! 263: IF( N.EQ.0 )
! 264: $ RETURN
! 265: *
! 266: * Machine Constants
! 267: *
! 268: SAFMIN = DLAMCH( 'Safe minimum' )
! 269: BIG = ONE / SAFMIN
! 270: CALL DLABAD( SAFMIN, BIG )
! 271: ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
! 272: SMALL = SAFMIN*N / ULP
! 273: BIG = ONE / SMALL
! 274: BIGNUM = ONE / ( SAFMIN*N )
! 275: *
! 276: * Compute the 1-norm of each column of the strictly upper triangular
! 277: * part of A and B to check for possible overflow in the triangular
! 278: * solver.
! 279: *
! 280: ANORM = ABS1( S( 1, 1 ) )
! 281: BNORM = ABS1( P( 1, 1 ) )
! 282: RWORK( 1 ) = ZERO
! 283: RWORK( N+1 ) = ZERO
! 284: DO 40 J = 2, N
! 285: RWORK( J ) = ZERO
! 286: RWORK( N+J ) = ZERO
! 287: DO 30 I = 1, J - 1
! 288: RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
! 289: RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
! 290: 30 CONTINUE
! 291: ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
! 292: BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
! 293: 40 CONTINUE
! 294: *
! 295: ASCALE = ONE / MAX( ANORM, SAFMIN )
! 296: BSCALE = ONE / MAX( BNORM, SAFMIN )
! 297: *
! 298: * Left eigenvectors
! 299: *
! 300: IF( COMPL ) THEN
! 301: IEIG = 0
! 302: *
! 303: * Main loop over eigenvalues
! 304: *
! 305: DO 140 JE = 1, N
! 306: IF( ILALL ) THEN
! 307: ILCOMP = .TRUE.
! 308: ELSE
! 309: ILCOMP = SELECT( JE )
! 310: END IF
! 311: IF( ILCOMP ) THEN
! 312: IEIG = IEIG + 1
! 313: *
! 314: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
! 315: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
! 316: *
! 317: * Singular matrix pencil -- return unit eigenvector
! 318: *
! 319: DO 50 JR = 1, N
! 320: VL( JR, IEIG ) = CZERO
! 321: 50 CONTINUE
! 322: VL( IEIG, IEIG ) = CONE
! 323: GO TO 140
! 324: END IF
! 325: *
! 326: * Non-singular eigenvalue:
! 327: * Compute coefficients a and b in
! 328: * H
! 329: * y ( a A - b B ) = 0
! 330: *
! 331: TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
! 332: $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
! 333: SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
! 334: SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
! 335: ACOEFF = SBETA*ASCALE
! 336: BCOEFF = SALPHA*BSCALE
! 337: *
! 338: * Scale to avoid underflow
! 339: *
! 340: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
! 341: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
! 342: $ SMALL
! 343: *
! 344: SCALE = ONE
! 345: IF( LSA )
! 346: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
! 347: IF( LSB )
! 348: $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
! 349: $ MIN( BNORM, BIG ) )
! 350: IF( LSA .OR. LSB ) THEN
! 351: SCALE = MIN( SCALE, ONE /
! 352: $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
! 353: $ ABS1( BCOEFF ) ) ) )
! 354: IF( LSA ) THEN
! 355: ACOEFF = ASCALE*( SCALE*SBETA )
! 356: ELSE
! 357: ACOEFF = SCALE*ACOEFF
! 358: END IF
! 359: IF( LSB ) THEN
! 360: BCOEFF = BSCALE*( SCALE*SALPHA )
! 361: ELSE
! 362: BCOEFF = SCALE*BCOEFF
! 363: END IF
! 364: END IF
! 365: *
! 366: ACOEFA = ABS( ACOEFF )
! 367: BCOEFA = ABS1( BCOEFF )
! 368: XMAX = ONE
! 369: DO 60 JR = 1, N
! 370: WORK( JR ) = CZERO
! 371: 60 CONTINUE
! 372: WORK( JE ) = CONE
! 373: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
! 374: *
! 375: * H
! 376: * Triangular solve of (a A - b B) y = 0
! 377: *
! 378: * H
! 379: * (rowwise in (a A - b B) , or columnwise in a A - b B)
! 380: *
! 381: DO 100 J = JE + 1, N
! 382: *
! 383: * Compute
! 384: * j-1
! 385: * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
! 386: * k=je
! 387: * (Scale if necessary)
! 388: *
! 389: TEMP = ONE / XMAX
! 390: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
! 391: $ TEMP ) THEN
! 392: DO 70 JR = JE, J - 1
! 393: WORK( JR ) = TEMP*WORK( JR )
! 394: 70 CONTINUE
! 395: XMAX = ONE
! 396: END IF
! 397: SUMA = CZERO
! 398: SUMB = CZERO
! 399: *
! 400: DO 80 JR = JE, J - 1
! 401: SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
! 402: SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
! 403: 80 CONTINUE
! 404: SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
! 405: *
! 406: * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
! 407: *
! 408: * with scaling and perturbation of the denominator
! 409: *
! 410: D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
! 411: IF( ABS1( D ).LE.DMIN )
! 412: $ D = DCMPLX( DMIN )
! 413: *
! 414: IF( ABS1( D ).LT.ONE ) THEN
! 415: IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
! 416: TEMP = ONE / ABS1( SUM )
! 417: DO 90 JR = JE, J - 1
! 418: WORK( JR ) = TEMP*WORK( JR )
! 419: 90 CONTINUE
! 420: XMAX = TEMP*XMAX
! 421: SUM = TEMP*SUM
! 422: END IF
! 423: END IF
! 424: WORK( J ) = ZLADIV( -SUM, D )
! 425: XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
! 426: 100 CONTINUE
! 427: *
! 428: * Back transform eigenvector if HOWMNY='B'.
! 429: *
! 430: IF( ILBACK ) THEN
! 431: CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
! 432: $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
! 433: ISRC = 2
! 434: IBEG = 1
! 435: ELSE
! 436: ISRC = 1
! 437: IBEG = JE
! 438: END IF
! 439: *
! 440: * Copy and scale eigenvector into column of VL
! 441: *
! 442: XMAX = ZERO
! 443: DO 110 JR = IBEG, N
! 444: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
! 445: 110 CONTINUE
! 446: *
! 447: IF( XMAX.GT.SAFMIN ) THEN
! 448: TEMP = ONE / XMAX
! 449: DO 120 JR = IBEG, N
! 450: VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
! 451: 120 CONTINUE
! 452: ELSE
! 453: IBEG = N + 1
! 454: END IF
! 455: *
! 456: DO 130 JR = 1, IBEG - 1
! 457: VL( JR, IEIG ) = CZERO
! 458: 130 CONTINUE
! 459: *
! 460: END IF
! 461: 140 CONTINUE
! 462: END IF
! 463: *
! 464: * Right eigenvectors
! 465: *
! 466: IF( COMPR ) THEN
! 467: IEIG = IM + 1
! 468: *
! 469: * Main loop over eigenvalues
! 470: *
! 471: DO 250 JE = N, 1, -1
! 472: IF( ILALL ) THEN
! 473: ILCOMP = .TRUE.
! 474: ELSE
! 475: ILCOMP = SELECT( JE )
! 476: END IF
! 477: IF( ILCOMP ) THEN
! 478: IEIG = IEIG - 1
! 479: *
! 480: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
! 481: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
! 482: *
! 483: * Singular matrix pencil -- return unit eigenvector
! 484: *
! 485: DO 150 JR = 1, N
! 486: VR( JR, IEIG ) = CZERO
! 487: 150 CONTINUE
! 488: VR( IEIG, IEIG ) = CONE
! 489: GO TO 250
! 490: END IF
! 491: *
! 492: * Non-singular eigenvalue:
! 493: * Compute coefficients a and b in
! 494: *
! 495: * ( a A - b B ) x = 0
! 496: *
! 497: TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
! 498: $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
! 499: SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
! 500: SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
! 501: ACOEFF = SBETA*ASCALE
! 502: BCOEFF = SALPHA*BSCALE
! 503: *
! 504: * Scale to avoid underflow
! 505: *
! 506: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
! 507: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
! 508: $ SMALL
! 509: *
! 510: SCALE = ONE
! 511: IF( LSA )
! 512: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
! 513: IF( LSB )
! 514: $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
! 515: $ MIN( BNORM, BIG ) )
! 516: IF( LSA .OR. LSB ) THEN
! 517: SCALE = MIN( SCALE, ONE /
! 518: $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
! 519: $ ABS1( BCOEFF ) ) ) )
! 520: IF( LSA ) THEN
! 521: ACOEFF = ASCALE*( SCALE*SBETA )
! 522: ELSE
! 523: ACOEFF = SCALE*ACOEFF
! 524: END IF
! 525: IF( LSB ) THEN
! 526: BCOEFF = BSCALE*( SCALE*SALPHA )
! 527: ELSE
! 528: BCOEFF = SCALE*BCOEFF
! 529: END IF
! 530: END IF
! 531: *
! 532: ACOEFA = ABS( ACOEFF )
! 533: BCOEFA = ABS1( BCOEFF )
! 534: XMAX = ONE
! 535: DO 160 JR = 1, N
! 536: WORK( JR ) = CZERO
! 537: 160 CONTINUE
! 538: WORK( JE ) = CONE
! 539: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
! 540: *
! 541: * Triangular solve of (a A - b B) x = 0 (columnwise)
! 542: *
! 543: * WORK(1:j-1) contains sums w,
! 544: * WORK(j+1:JE) contains x
! 545: *
! 546: DO 170 JR = 1, JE - 1
! 547: WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
! 548: 170 CONTINUE
! 549: WORK( JE ) = CONE
! 550: *
! 551: DO 210 J = JE - 1, 1, -1
! 552: *
! 553: * Form x(j) := - w(j) / d
! 554: * with scaling and perturbation of the denominator
! 555: *
! 556: D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
! 557: IF( ABS1( D ).LE.DMIN )
! 558: $ D = DCMPLX( DMIN )
! 559: *
! 560: IF( ABS1( D ).LT.ONE ) THEN
! 561: IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
! 562: TEMP = ONE / ABS1( WORK( J ) )
! 563: DO 180 JR = 1, JE
! 564: WORK( JR ) = TEMP*WORK( JR )
! 565: 180 CONTINUE
! 566: END IF
! 567: END IF
! 568: *
! 569: WORK( J ) = ZLADIV( -WORK( J ), D )
! 570: *
! 571: IF( J.GT.1 ) THEN
! 572: *
! 573: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
! 574: *
! 575: IF( ABS1( WORK( J ) ).GT.ONE ) THEN
! 576: TEMP = ONE / ABS1( WORK( J ) )
! 577: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
! 578: $ BIGNUM*TEMP ) THEN
! 579: DO 190 JR = 1, JE
! 580: WORK( JR ) = TEMP*WORK( JR )
! 581: 190 CONTINUE
! 582: END IF
! 583: END IF
! 584: *
! 585: CA = ACOEFF*WORK( J )
! 586: CB = BCOEFF*WORK( J )
! 587: DO 200 JR = 1, J - 1
! 588: WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
! 589: $ CB*P( JR, J )
! 590: 200 CONTINUE
! 591: END IF
! 592: 210 CONTINUE
! 593: *
! 594: * Back transform eigenvector if HOWMNY='B'.
! 595: *
! 596: IF( ILBACK ) THEN
! 597: CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
! 598: $ CZERO, WORK( N+1 ), 1 )
! 599: ISRC = 2
! 600: IEND = N
! 601: ELSE
! 602: ISRC = 1
! 603: IEND = JE
! 604: END IF
! 605: *
! 606: * Copy and scale eigenvector into column of VR
! 607: *
! 608: XMAX = ZERO
! 609: DO 220 JR = 1, IEND
! 610: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
! 611: 220 CONTINUE
! 612: *
! 613: IF( XMAX.GT.SAFMIN ) THEN
! 614: TEMP = ONE / XMAX
! 615: DO 230 JR = 1, IEND
! 616: VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
! 617: 230 CONTINUE
! 618: ELSE
! 619: IEND = 0
! 620: END IF
! 621: *
! 622: DO 240 JR = IEND + 1, N
! 623: VR( JR, IEIG ) = CZERO
! 624: 240 CONTINUE
! 625: *
! 626: END IF
! 627: 250 CONTINUE
! 628: END IF
! 629: *
! 630: RETURN
! 631: *
! 632: * End of ZTGEVC
! 633: *
! 634: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgevc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack