Annotation of rpl/lapack/lapack/zggevx.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
! 2: $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
! 3: $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
! 4: $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
! 5: *
! 6: * -- LAPACK driver routine (version 3.2) --
! 7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 9: * November 2006
! 10: *
! 11: * .. Scalar Arguments ..
! 12: CHARACTER BALANC, JOBVL, JOBVR, SENSE
! 13: INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
! 14: DOUBLE PRECISION ABNRM, BBNRM
! 15: * ..
! 16: * .. Array Arguments ..
! 17: LOGICAL BWORK( * )
! 18: INTEGER IWORK( * )
! 19: DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
! 20: $ RSCALE( * ), RWORK( * )
! 21: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
! 22: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
! 23: $ WORK( * )
! 24: * ..
! 25: *
! 26: * Purpose
! 27: * =======
! 28: *
! 29: * ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
! 30: * (A,B) the generalized eigenvalues, and optionally, the left and/or
! 31: * right generalized eigenvectors.
! 32: *
! 33: * Optionally, it also computes a balancing transformation to improve
! 34: * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
! 35: * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
! 36: * the eigenvalues (RCONDE), and reciprocal condition numbers for the
! 37: * right eigenvectors (RCONDV).
! 38: *
! 39: * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
! 40: * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
! 41: * singular. It is usually represented as the pair (alpha,beta), as
! 42: * there is a reasonable interpretation for beta=0, and even for both
! 43: * being zero.
! 44: *
! 45: * The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
! 46: * of (A,B) satisfies
! 47: * A * v(j) = lambda(j) * B * v(j) .
! 48: * The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
! 49: * of (A,B) satisfies
! 50: * u(j)**H * A = lambda(j) * u(j)**H * B.
! 51: * where u(j)**H is the conjugate-transpose of u(j).
! 52: *
! 53: *
! 54: * Arguments
! 55: * =========
! 56: *
! 57: * BALANC (input) CHARACTER*1
! 58: * Specifies the balance option to be performed:
! 59: * = 'N': do not diagonally scale or permute;
! 60: * = 'P': permute only;
! 61: * = 'S': scale only;
! 62: * = 'B': both permute and scale.
! 63: * Computed reciprocal condition numbers will be for the
! 64: * matrices after permuting and/or balancing. Permuting does
! 65: * not change condition numbers (in exact arithmetic), but
! 66: * balancing does.
! 67: *
! 68: * JOBVL (input) CHARACTER*1
! 69: * = 'N': do not compute the left generalized eigenvectors;
! 70: * = 'V': compute the left generalized eigenvectors.
! 71: *
! 72: * JOBVR (input) CHARACTER*1
! 73: * = 'N': do not compute the right generalized eigenvectors;
! 74: * = 'V': compute the right generalized eigenvectors.
! 75: *
! 76: * SENSE (input) CHARACTER*1
! 77: * Determines which reciprocal condition numbers are computed.
! 78: * = 'N': none are computed;
! 79: * = 'E': computed for eigenvalues only;
! 80: * = 'V': computed for eigenvectors only;
! 81: * = 'B': computed for eigenvalues and eigenvectors.
! 82: *
! 83: * N (input) INTEGER
! 84: * The order of the matrices A, B, VL, and VR. N >= 0.
! 85: *
! 86: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
! 87: * On entry, the matrix A in the pair (A,B).
! 88: * On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
! 89: * or both, then A contains the first part of the complex Schur
! 90: * form of the "balanced" versions of the input A and B.
! 91: *
! 92: * LDA (input) INTEGER
! 93: * The leading dimension of A. LDA >= max(1,N).
! 94: *
! 95: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
! 96: * On entry, the matrix B in the pair (A,B).
! 97: * On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
! 98: * or both, then B contains the second part of the complex
! 99: * Schur form of the "balanced" versions of the input A and B.
! 100: *
! 101: * LDB (input) INTEGER
! 102: * The leading dimension of B. LDB >= max(1,N).
! 103: *
! 104: * ALPHA (output) COMPLEX*16 array, dimension (N)
! 105: * BETA (output) COMPLEX*16 array, dimension (N)
! 106: * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
! 107: * eigenvalues.
! 108: *
! 109: * Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
! 110: * underflow, and BETA(j) may even be zero. Thus, the user
! 111: * should avoid naively computing the ratio ALPHA/BETA.
! 112: * However, ALPHA will be always less than and usually
! 113: * comparable with norm(A) in magnitude, and BETA always less
! 114: * than and usually comparable with norm(B).
! 115: *
! 116: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
! 117: * If JOBVL = 'V', the left generalized eigenvectors u(j) are
! 118: * stored one after another in the columns of VL, in the same
! 119: * order as their eigenvalues.
! 120: * Each eigenvector will be scaled so the largest component
! 121: * will have abs(real part) + abs(imag. part) = 1.
! 122: * Not referenced if JOBVL = 'N'.
! 123: *
! 124: * LDVL (input) INTEGER
! 125: * The leading dimension of the matrix VL. LDVL >= 1, and
! 126: * if JOBVL = 'V', LDVL >= N.
! 127: *
! 128: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
! 129: * If JOBVR = 'V', the right generalized eigenvectors v(j) are
! 130: * stored one after another in the columns of VR, in the same
! 131: * order as their eigenvalues.
! 132: * Each eigenvector will be scaled so the largest component
! 133: * will have abs(real part) + abs(imag. part) = 1.
! 134: * Not referenced if JOBVR = 'N'.
! 135: *
! 136: * LDVR (input) INTEGER
! 137: * The leading dimension of the matrix VR. LDVR >= 1, and
! 138: * if JOBVR = 'V', LDVR >= N.
! 139: *
! 140: * ILO (output) INTEGER
! 141: * IHI (output) INTEGER
! 142: * ILO and IHI are integer values such that on exit
! 143: * A(i,j) = 0 and B(i,j) = 0 if i > j and
! 144: * j = 1,...,ILO-1 or i = IHI+1,...,N.
! 145: * If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
! 146: *
! 147: * LSCALE (output) DOUBLE PRECISION array, dimension (N)
! 148: * Details of the permutations and scaling factors applied
! 149: * to the left side of A and B. If PL(j) is the index of the
! 150: * row interchanged with row j, and DL(j) is the scaling
! 151: * factor applied to row j, then
! 152: * LSCALE(j) = PL(j) for j = 1,...,ILO-1
! 153: * = DL(j) for j = ILO,...,IHI
! 154: * = PL(j) for j = IHI+1,...,N.
! 155: * The order in which the interchanges are made is N to IHI+1,
! 156: * then 1 to ILO-1.
! 157: *
! 158: * RSCALE (output) DOUBLE PRECISION array, dimension (N)
! 159: * Details of the permutations and scaling factors applied
! 160: * to the right side of A and B. If PR(j) is the index of the
! 161: * column interchanged with column j, and DR(j) is the scaling
! 162: * factor applied to column j, then
! 163: * RSCALE(j) = PR(j) for j = 1,...,ILO-1
! 164: * = DR(j) for j = ILO,...,IHI
! 165: * = PR(j) for j = IHI+1,...,N
! 166: * The order in which the interchanges are made is N to IHI+1,
! 167: * then 1 to ILO-1.
! 168: *
! 169: * ABNRM (output) DOUBLE PRECISION
! 170: * The one-norm of the balanced matrix A.
! 171: *
! 172: * BBNRM (output) DOUBLE PRECISION
! 173: * The one-norm of the balanced matrix B.
! 174: *
! 175: * RCONDE (output) DOUBLE PRECISION array, dimension (N)
! 176: * If SENSE = 'E' or 'B', the reciprocal condition numbers of
! 177: * the eigenvalues, stored in consecutive elements of the array.
! 178: * If SENSE = 'N' or 'V', RCONDE is not referenced.
! 179: *
! 180: * RCONDV (output) DOUBLE PRECISION array, dimension (N)
! 181: * If JOB = 'V' or 'B', the estimated reciprocal condition
! 182: * numbers of the eigenvectors, stored in consecutive elements
! 183: * of the array. If the eigenvalues cannot be reordered to
! 184: * compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
! 185: * when the true value would be very small anyway.
! 186: * If SENSE = 'N' or 'E', RCONDV is not referenced.
! 187: *
! 188: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 189: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 190: *
! 191: * LWORK (input) INTEGER
! 192: * The dimension of the array WORK. LWORK >= max(1,2*N).
! 193: * If SENSE = 'E', LWORK >= max(1,4*N).
! 194: * If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
! 195: *
! 196: * If LWORK = -1, then a workspace query is assumed; the routine
! 197: * only calculates the optimal size of the WORK array, returns
! 198: * this value as the first entry of the WORK array, and no error
! 199: * message related to LWORK is issued by XERBLA.
! 200: *
! 201: * RWORK (workspace) REAL array, dimension (lrwork)
! 202: * lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
! 203: * and at least max(1,2*N) otherwise.
! 204: * Real workspace.
! 205: *
! 206: * IWORK (workspace) INTEGER array, dimension (N+2)
! 207: * If SENSE = 'E', IWORK is not referenced.
! 208: *
! 209: * BWORK (workspace) LOGICAL array, dimension (N)
! 210: * If SENSE = 'N', BWORK is not referenced.
! 211: *
! 212: * INFO (output) INTEGER
! 213: * = 0: successful exit
! 214: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 215: * = 1,...,N:
! 216: * The QZ iteration failed. No eigenvectors have been
! 217: * calculated, but ALPHA(j) and BETA(j) should be correct
! 218: * for j=INFO+1,...,N.
! 219: * > N: =N+1: other than QZ iteration failed in ZHGEQZ.
! 220: * =N+2: error return from ZTGEVC.
! 221: *
! 222: * Further Details
! 223: * ===============
! 224: *
! 225: * Balancing a matrix pair (A,B) includes, first, permuting rows and
! 226: * columns to isolate eigenvalues, second, applying diagonal similarity
! 227: * transformation to the rows and columns to make the rows and columns
! 228: * as close in norm as possible. The computed reciprocal condition
! 229: * numbers correspond to the balanced matrix. Permuting rows and columns
! 230: * will not change the condition numbers (in exact arithmetic) but
! 231: * diagonal scaling will. For further explanation of balancing, see
! 232: * section 4.11.1.2 of LAPACK Users' Guide.
! 233: *
! 234: * An approximate error bound on the chordal distance between the i-th
! 235: * computed generalized eigenvalue w and the corresponding exact
! 236: * eigenvalue lambda is
! 237: *
! 238: * chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
! 239: *
! 240: * An approximate error bound for the angle between the i-th computed
! 241: * eigenvector VL(i) or VR(i) is given by
! 242: *
! 243: * EPS * norm(ABNRM, BBNRM) / DIF(i).
! 244: *
! 245: * For further explanation of the reciprocal condition numbers RCONDE
! 246: * and RCONDV, see section 4.11 of LAPACK User's Guide.
! 247: *
! 248: * .. Parameters ..
! 249: DOUBLE PRECISION ZERO, ONE
! 250: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 251: COMPLEX*16 CZERO, CONE
! 252: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 253: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 254: * ..
! 255: * .. Local Scalars ..
! 256: LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
! 257: $ WANTSB, WANTSE, WANTSN, WANTSV
! 258: CHARACTER CHTEMP
! 259: INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
! 260: $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
! 261: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
! 262: $ SMLNUM, TEMP
! 263: COMPLEX*16 X
! 264: * ..
! 265: * .. Local Arrays ..
! 266: LOGICAL LDUMMA( 1 )
! 267: * ..
! 268: * .. External Subroutines ..
! 269: EXTERNAL DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
! 270: $ ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
! 271: $ ZTGSNA, ZUNGQR, ZUNMQR
! 272: * ..
! 273: * .. External Functions ..
! 274: LOGICAL LSAME
! 275: INTEGER ILAENV
! 276: DOUBLE PRECISION DLAMCH, ZLANGE
! 277: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 278: * ..
! 279: * .. Intrinsic Functions ..
! 280: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
! 281: * ..
! 282: * .. Statement Functions ..
! 283: DOUBLE PRECISION ABS1
! 284: * ..
! 285: * .. Statement Function definitions ..
! 286: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
! 287: * ..
! 288: * .. Executable Statements ..
! 289: *
! 290: * Decode the input arguments
! 291: *
! 292: IF( LSAME( JOBVL, 'N' ) ) THEN
! 293: IJOBVL = 1
! 294: ILVL = .FALSE.
! 295: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
! 296: IJOBVL = 2
! 297: ILVL = .TRUE.
! 298: ELSE
! 299: IJOBVL = -1
! 300: ILVL = .FALSE.
! 301: END IF
! 302: *
! 303: IF( LSAME( JOBVR, 'N' ) ) THEN
! 304: IJOBVR = 1
! 305: ILVR = .FALSE.
! 306: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
! 307: IJOBVR = 2
! 308: ILVR = .TRUE.
! 309: ELSE
! 310: IJOBVR = -1
! 311: ILVR = .FALSE.
! 312: END IF
! 313: ILV = ILVL .OR. ILVR
! 314: *
! 315: NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
! 316: WANTSN = LSAME( SENSE, 'N' )
! 317: WANTSE = LSAME( SENSE, 'E' )
! 318: WANTSV = LSAME( SENSE, 'V' )
! 319: WANTSB = LSAME( SENSE, 'B' )
! 320: *
! 321: * Test the input arguments
! 322: *
! 323: INFO = 0
! 324: LQUERY = ( LWORK.EQ.-1 )
! 325: IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
! 326: $ LSAME( BALANC, 'B' ) ) ) THEN
! 327: INFO = -1
! 328: ELSE IF( IJOBVL.LE.0 ) THEN
! 329: INFO = -2
! 330: ELSE IF( IJOBVR.LE.0 ) THEN
! 331: INFO = -3
! 332: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
! 333: $ THEN
! 334: INFO = -4
! 335: ELSE IF( N.LT.0 ) THEN
! 336: INFO = -5
! 337: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 338: INFO = -7
! 339: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 340: INFO = -9
! 341: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
! 342: INFO = -13
! 343: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
! 344: INFO = -15
! 345: END IF
! 346: *
! 347: * Compute workspace
! 348: * (Note: Comments in the code beginning "Workspace:" describe the
! 349: * minimal amount of workspace needed at that point in the code,
! 350: * as well as the preferred amount for good performance.
! 351: * NB refers to the optimal block size for the immediately
! 352: * following subroutine, as returned by ILAENV. The workspace is
! 353: * computed assuming ILO = 1 and IHI = N, the worst case.)
! 354: *
! 355: IF( INFO.EQ.0 ) THEN
! 356: IF( N.EQ.0 ) THEN
! 357: MINWRK = 1
! 358: MAXWRK = 1
! 359: ELSE
! 360: MINWRK = 2*N
! 361: IF( WANTSE ) THEN
! 362: MINWRK = 4*N
! 363: ELSE IF( WANTSV .OR. WANTSB ) THEN
! 364: MINWRK = 2*N*( N + 1)
! 365: END IF
! 366: MAXWRK = MINWRK
! 367: MAXWRK = MAX( MAXWRK,
! 368: $ N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
! 369: MAXWRK = MAX( MAXWRK,
! 370: $ N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
! 371: IF( ILVL ) THEN
! 372: MAXWRK = MAX( MAXWRK, N +
! 373: $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
! 374: END IF
! 375: END IF
! 376: WORK( 1 ) = MAXWRK
! 377: *
! 378: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 379: INFO = -25
! 380: END IF
! 381: END IF
! 382: *
! 383: IF( INFO.NE.0 ) THEN
! 384: CALL XERBLA( 'ZGGEVX', -INFO )
! 385: RETURN
! 386: ELSE IF( LQUERY ) THEN
! 387: RETURN
! 388: END IF
! 389: *
! 390: * Quick return if possible
! 391: *
! 392: IF( N.EQ.0 )
! 393: $ RETURN
! 394: *
! 395: * Get machine constants
! 396: *
! 397: EPS = DLAMCH( 'P' )
! 398: SMLNUM = DLAMCH( 'S' )
! 399: BIGNUM = ONE / SMLNUM
! 400: CALL DLABAD( SMLNUM, BIGNUM )
! 401: SMLNUM = SQRT( SMLNUM ) / EPS
! 402: BIGNUM = ONE / SMLNUM
! 403: *
! 404: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 405: *
! 406: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
! 407: ILASCL = .FALSE.
! 408: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 409: ANRMTO = SMLNUM
! 410: ILASCL = .TRUE.
! 411: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 412: ANRMTO = BIGNUM
! 413: ILASCL = .TRUE.
! 414: END IF
! 415: IF( ILASCL )
! 416: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
! 417: *
! 418: * Scale B if max element outside range [SMLNUM,BIGNUM]
! 419: *
! 420: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
! 421: ILBSCL = .FALSE.
! 422: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 423: BNRMTO = SMLNUM
! 424: ILBSCL = .TRUE.
! 425: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 426: BNRMTO = BIGNUM
! 427: ILBSCL = .TRUE.
! 428: END IF
! 429: IF( ILBSCL )
! 430: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
! 431: *
! 432: * Permute and/or balance the matrix pair (A,B)
! 433: * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
! 434: *
! 435: CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
! 436: $ RWORK, IERR )
! 437: *
! 438: * Compute ABNRM and BBNRM
! 439: *
! 440: ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
! 441: IF( ILASCL ) THEN
! 442: RWORK( 1 ) = ABNRM
! 443: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
! 444: $ IERR )
! 445: ABNRM = RWORK( 1 )
! 446: END IF
! 447: *
! 448: BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
! 449: IF( ILBSCL ) THEN
! 450: RWORK( 1 ) = BBNRM
! 451: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
! 452: $ IERR )
! 453: BBNRM = RWORK( 1 )
! 454: END IF
! 455: *
! 456: * Reduce B to triangular form (QR decomposition of B)
! 457: * (Complex Workspace: need N, prefer N*NB )
! 458: *
! 459: IROWS = IHI + 1 - ILO
! 460: IF( ILV .OR. .NOT.WANTSN ) THEN
! 461: ICOLS = N + 1 - ILO
! 462: ELSE
! 463: ICOLS = IROWS
! 464: END IF
! 465: ITAU = 1
! 466: IWRK = ITAU + IROWS
! 467: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
! 468: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
! 469: *
! 470: * Apply the unitary transformation to A
! 471: * (Complex Workspace: need N, prefer N*NB)
! 472: *
! 473: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
! 474: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
! 475: $ LWORK+1-IWRK, IERR )
! 476: *
! 477: * Initialize VL and/or VR
! 478: * (Workspace: need N, prefer N*NB)
! 479: *
! 480: IF( ILVL ) THEN
! 481: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
! 482: IF( IROWS.GT.1 ) THEN
! 483: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
! 484: $ VL( ILO+1, ILO ), LDVL )
! 485: END IF
! 486: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
! 487: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
! 488: END IF
! 489: *
! 490: IF( ILVR )
! 491: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
! 492: *
! 493: * Reduce to generalized Hessenberg form
! 494: * (Workspace: none needed)
! 495: *
! 496: IF( ILV .OR. .NOT.WANTSN ) THEN
! 497: *
! 498: * Eigenvectors requested -- work on whole matrix.
! 499: *
! 500: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
! 501: $ LDVL, VR, LDVR, IERR )
! 502: ELSE
! 503: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
! 504: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
! 505: END IF
! 506: *
! 507: * Perform QZ algorithm (Compute eigenvalues, and optionally, the
! 508: * Schur forms and Schur vectors)
! 509: * (Complex Workspace: need N)
! 510: * (Real Workspace: need N)
! 511: *
! 512: IWRK = ITAU
! 513: IF( ILV .OR. .NOT.WANTSN ) THEN
! 514: CHTEMP = 'S'
! 515: ELSE
! 516: CHTEMP = 'E'
! 517: END IF
! 518: *
! 519: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
! 520: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
! 521: $ LWORK+1-IWRK, RWORK, IERR )
! 522: IF( IERR.NE.0 ) THEN
! 523: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
! 524: INFO = IERR
! 525: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
! 526: INFO = IERR - N
! 527: ELSE
! 528: INFO = N + 1
! 529: END IF
! 530: GO TO 90
! 531: END IF
! 532: *
! 533: * Compute Eigenvectors and estimate condition numbers if desired
! 534: * ZTGEVC: (Complex Workspace: need 2*N )
! 535: * (Real Workspace: need 2*N )
! 536: * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
! 537: * (Integer Workspace: need N+2 )
! 538: *
! 539: IF( ILV .OR. .NOT.WANTSN ) THEN
! 540: IF( ILV ) THEN
! 541: IF( ILVL ) THEN
! 542: IF( ILVR ) THEN
! 543: CHTEMP = 'B'
! 544: ELSE
! 545: CHTEMP = 'L'
! 546: END IF
! 547: ELSE
! 548: CHTEMP = 'R'
! 549: END IF
! 550: *
! 551: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
! 552: $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
! 553: $ IERR )
! 554: IF( IERR.NE.0 ) THEN
! 555: INFO = N + 2
! 556: GO TO 90
! 557: END IF
! 558: END IF
! 559: *
! 560: IF( .NOT.WANTSN ) THEN
! 561: *
! 562: * compute eigenvectors (DTGEVC) and estimate condition
! 563: * numbers (DTGSNA). Note that the definition of the condition
! 564: * number is not invariant under transformation (u,v) to
! 565: * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
! 566: * Schur form (S,T), Q and Z are orthogonal matrices. In order
! 567: * to avoid using extra 2*N*N workspace, we have to
! 568: * re-calculate eigenvectors and estimate the condition numbers
! 569: * one at a time.
! 570: *
! 571: DO 20 I = 1, N
! 572: *
! 573: DO 10 J = 1, N
! 574: BWORK( J ) = .FALSE.
! 575: 10 CONTINUE
! 576: BWORK( I ) = .TRUE.
! 577: *
! 578: IWRK = N + 1
! 579: IWRK1 = IWRK + N
! 580: *
! 581: IF( WANTSE .OR. WANTSB ) THEN
! 582: CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
! 583: $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
! 584: $ WORK( IWRK1 ), RWORK, IERR )
! 585: IF( IERR.NE.0 ) THEN
! 586: INFO = N + 2
! 587: GO TO 90
! 588: END IF
! 589: END IF
! 590: *
! 591: CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
! 592: $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
! 593: $ RCONDV( I ), 1, M, WORK( IWRK1 ),
! 594: $ LWORK-IWRK1+1, IWORK, IERR )
! 595: *
! 596: 20 CONTINUE
! 597: END IF
! 598: END IF
! 599: *
! 600: * Undo balancing on VL and VR and normalization
! 601: * (Workspace: none needed)
! 602: *
! 603: IF( ILVL ) THEN
! 604: CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
! 605: $ LDVL, IERR )
! 606: *
! 607: DO 50 JC = 1, N
! 608: TEMP = ZERO
! 609: DO 30 JR = 1, N
! 610: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
! 611: 30 CONTINUE
! 612: IF( TEMP.LT.SMLNUM )
! 613: $ GO TO 50
! 614: TEMP = ONE / TEMP
! 615: DO 40 JR = 1, N
! 616: VL( JR, JC ) = VL( JR, JC )*TEMP
! 617: 40 CONTINUE
! 618: 50 CONTINUE
! 619: END IF
! 620: *
! 621: IF( ILVR ) THEN
! 622: CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
! 623: $ LDVR, IERR )
! 624: DO 80 JC = 1, N
! 625: TEMP = ZERO
! 626: DO 60 JR = 1, N
! 627: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
! 628: 60 CONTINUE
! 629: IF( TEMP.LT.SMLNUM )
! 630: $ GO TO 80
! 631: TEMP = ONE / TEMP
! 632: DO 70 JR = 1, N
! 633: VR( JR, JC ) = VR( JR, JC )*TEMP
! 634: 70 CONTINUE
! 635: 80 CONTINUE
! 636: END IF
! 637: *
! 638: * Undo scaling if necessary
! 639: *
! 640: IF( ILASCL )
! 641: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
! 642: *
! 643: IF( ILBSCL )
! 644: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
! 645: *
! 646: 90 CONTINUE
! 647: WORK( 1 ) = MAXWRK
! 648: *
! 649: RETURN
! 650: *
! 651: * End of ZGGEVX
! 652: *
! 653: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggevx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack