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