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