Annotation of rpl/lapack/lapack/zggev.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
! 2: $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK driver 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 JOBVL, JOBVR
! 11: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION RWORK( * )
! 15: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
! 16: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
! 17: $ WORK( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
! 24: * (A,B), the generalized eigenvalues, and optionally, the left and/or
! 25: * right generalized eigenvectors.
! 26: *
! 27: * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
! 28: * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
! 29: * singular. It is usually represented as the pair (alpha,beta), as
! 30: * there is a reasonable interpretation for beta=0, and even for both
! 31: * being zero.
! 32: *
! 33: * The right generalized eigenvector v(j) corresponding to the
! 34: * generalized eigenvalue lambda(j) of (A,B) satisfies
! 35: *
! 36: * A * v(j) = lambda(j) * B * v(j).
! 37: *
! 38: * The left generalized eigenvector u(j) corresponding to the
! 39: * generalized eigenvalues lambda(j) of (A,B) satisfies
! 40: *
! 41: * u(j)**H * A = lambda(j) * u(j)**H * B
! 42: *
! 43: * where u(j)**H is the conjugate-transpose of u(j).
! 44: *
! 45: * Arguments
! 46: * =========
! 47: *
! 48: * JOBVL (input) CHARACTER*1
! 49: * = 'N': do not compute the left generalized eigenvectors;
! 50: * = 'V': compute the left generalized eigenvectors.
! 51: *
! 52: * JOBVR (input) CHARACTER*1
! 53: * = 'N': do not compute the right generalized eigenvectors;
! 54: * = 'V': compute the right generalized eigenvectors.
! 55: *
! 56: * N (input) INTEGER
! 57: * The order of the matrices A, B, VL, and VR. N >= 0.
! 58: *
! 59: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
! 60: * On entry, the matrix A in the pair (A,B).
! 61: * On exit, A has been overwritten.
! 62: *
! 63: * LDA (input) INTEGER
! 64: * The leading dimension of A. LDA >= max(1,N).
! 65: *
! 66: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
! 67: * On entry, the matrix B in the pair (A,B).
! 68: * On exit, B has been overwritten.
! 69: *
! 70: * LDB (input) INTEGER
! 71: * The leading dimension of B. LDB >= max(1,N).
! 72: *
! 73: * ALPHA (output) COMPLEX*16 array, dimension (N)
! 74: * BETA (output) COMPLEX*16 array, dimension (N)
! 75: * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
! 76: * generalized eigenvalues.
! 77: *
! 78: * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
! 79: * underflow, and BETA(j) may even be zero. Thus, the user
! 80: * should avoid naively computing the ratio alpha/beta.
! 81: * However, ALPHA will be always less than and usually
! 82: * comparable with norm(A) in magnitude, and BETA always less
! 83: * than and usually comparable with norm(B).
! 84: *
! 85: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
! 86: * If JOBVL = 'V', the left generalized eigenvectors u(j) are
! 87: * stored one after another in the columns of VL, in the same
! 88: * order as their eigenvalues.
! 89: * Each eigenvector is scaled so the largest component has
! 90: * abs(real part) + abs(imag. part) = 1.
! 91: * Not referenced if JOBVL = 'N'.
! 92: *
! 93: * LDVL (input) INTEGER
! 94: * The leading dimension of the matrix VL. LDVL >= 1, and
! 95: * if JOBVL = 'V', LDVL >= N.
! 96: *
! 97: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
! 98: * If JOBVR = 'V', the right generalized eigenvectors v(j) are
! 99: * stored one after another in the columns of VR, in the same
! 100: * order as their eigenvalues.
! 101: * Each eigenvector is scaled so the largest component has
! 102: * abs(real part) + abs(imag. part) = 1.
! 103: * Not referenced if JOBVR = 'N'.
! 104: *
! 105: * LDVR (input) INTEGER
! 106: * The leading dimension of the matrix VR. LDVR >= 1, and
! 107: * if JOBVR = 'V', LDVR >= N.
! 108: *
! 109: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 110: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 111: *
! 112: * LWORK (input) INTEGER
! 113: * The dimension of the array WORK. LWORK >= max(1,2*N).
! 114: * For good performance, LWORK must generally be larger.
! 115: *
! 116: * If LWORK = -1, then a workspace query is assumed; the routine
! 117: * only calculates the optimal size of the WORK array, returns
! 118: * this value as the first entry of the WORK array, and no error
! 119: * message related to LWORK is issued by XERBLA.
! 120: *
! 121: * RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N)
! 122: *
! 123: * INFO (output) INTEGER
! 124: * = 0: successful exit
! 125: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 126: * =1,...,N:
! 127: * The QZ iteration failed. No eigenvectors have been
! 128: * calculated, but ALPHA(j) and BETA(j) should be
! 129: * correct for j=INFO+1,...,N.
! 130: * > N: =N+1: other then QZ iteration failed in DHGEQZ,
! 131: * =N+2: error return from DTGEVC.
! 132: *
! 133: * =====================================================================
! 134: *
! 135: * .. Parameters ..
! 136: DOUBLE PRECISION ZERO, ONE
! 137: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 138: COMPLEX*16 CZERO, CONE
! 139: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
! 140: $ CONE = ( 1.0D0, 0.0D0 ) )
! 141: * ..
! 142: * .. Local Scalars ..
! 143: LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
! 144: CHARACTER CHTEMP
! 145: INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
! 146: $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
! 147: $ LWKMIN, LWKOPT
! 148: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
! 149: $ SMLNUM, TEMP
! 150: COMPLEX*16 X
! 151: * ..
! 152: * .. Local Arrays ..
! 153: LOGICAL LDUMMA( 1 )
! 154: * ..
! 155: * .. External Subroutines ..
! 156: EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
! 157: $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
! 158: $ ZUNMQR
! 159: * ..
! 160: * .. External Functions ..
! 161: LOGICAL LSAME
! 162: INTEGER ILAENV
! 163: DOUBLE PRECISION DLAMCH, ZLANGE
! 164: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 165: * ..
! 166: * .. Intrinsic Functions ..
! 167: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
! 168: * ..
! 169: * .. Statement Functions ..
! 170: DOUBLE PRECISION ABS1
! 171: * ..
! 172: * .. Statement Function definitions ..
! 173: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
! 174: * ..
! 175: * .. Executable Statements ..
! 176: *
! 177: * Decode the input arguments
! 178: *
! 179: IF( LSAME( JOBVL, 'N' ) ) THEN
! 180: IJOBVL = 1
! 181: ILVL = .FALSE.
! 182: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
! 183: IJOBVL = 2
! 184: ILVL = .TRUE.
! 185: ELSE
! 186: IJOBVL = -1
! 187: ILVL = .FALSE.
! 188: END IF
! 189: *
! 190: IF( LSAME( JOBVR, 'N' ) ) THEN
! 191: IJOBVR = 1
! 192: ILVR = .FALSE.
! 193: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
! 194: IJOBVR = 2
! 195: ILVR = .TRUE.
! 196: ELSE
! 197: IJOBVR = -1
! 198: ILVR = .FALSE.
! 199: END IF
! 200: ILV = ILVL .OR. ILVR
! 201: *
! 202: * Test the input arguments
! 203: *
! 204: INFO = 0
! 205: LQUERY = ( LWORK.EQ.-1 )
! 206: IF( IJOBVL.LE.0 ) THEN
! 207: INFO = -1
! 208: ELSE IF( IJOBVR.LE.0 ) THEN
! 209: INFO = -2
! 210: ELSE IF( N.LT.0 ) THEN
! 211: INFO = -3
! 212: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 213: INFO = -5
! 214: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 215: INFO = -7
! 216: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
! 217: INFO = -11
! 218: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
! 219: INFO = -13
! 220: END IF
! 221: *
! 222: * Compute workspace
! 223: * (Note: Comments in the code beginning "Workspace:" describe the
! 224: * minimal amount of workspace needed at that point in the code,
! 225: * as well as the preferred amount for good performance.
! 226: * NB refers to the optimal block size for the immediately
! 227: * following subroutine, as returned by ILAENV. The workspace is
! 228: * computed assuming ILO = 1 and IHI = N, the worst case.)
! 229: *
! 230: IF( INFO.EQ.0 ) THEN
! 231: LWKMIN = MAX( 1, 2*N )
! 232: LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
! 233: LWKOPT = MAX( LWKOPT, N +
! 234: $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
! 235: IF( ILVL ) THEN
! 236: LWKOPT = MAX( LWKOPT, N +
! 237: $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
! 238: END IF
! 239: WORK( 1 ) = LWKOPT
! 240: *
! 241: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
! 242: $ INFO = -15
! 243: END IF
! 244: *
! 245: IF( INFO.NE.0 ) THEN
! 246: CALL XERBLA( 'ZGGEV ', -INFO )
! 247: RETURN
! 248: ELSE IF( LQUERY ) THEN
! 249: RETURN
! 250: END IF
! 251: *
! 252: * Quick return if possible
! 253: *
! 254: IF( N.EQ.0 )
! 255: $ RETURN
! 256: *
! 257: * Get machine constants
! 258: *
! 259: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
! 260: SMLNUM = DLAMCH( 'S' )
! 261: BIGNUM = ONE / SMLNUM
! 262: CALL DLABAD( SMLNUM, BIGNUM )
! 263: SMLNUM = SQRT( SMLNUM ) / EPS
! 264: BIGNUM = ONE / SMLNUM
! 265: *
! 266: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 267: *
! 268: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
! 269: ILASCL = .FALSE.
! 270: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 271: ANRMTO = SMLNUM
! 272: ILASCL = .TRUE.
! 273: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 274: ANRMTO = BIGNUM
! 275: ILASCL = .TRUE.
! 276: END IF
! 277: IF( ILASCL )
! 278: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
! 279: *
! 280: * Scale B if max element outside range [SMLNUM,BIGNUM]
! 281: *
! 282: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
! 283: ILBSCL = .FALSE.
! 284: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 285: BNRMTO = SMLNUM
! 286: ILBSCL = .TRUE.
! 287: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 288: BNRMTO = BIGNUM
! 289: ILBSCL = .TRUE.
! 290: END IF
! 291: IF( ILBSCL )
! 292: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
! 293: *
! 294: * Permute the matrices A, B to isolate eigenvalues if possible
! 295: * (Real Workspace: need 6*N)
! 296: *
! 297: ILEFT = 1
! 298: IRIGHT = N + 1
! 299: IRWRK = IRIGHT + N
! 300: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
! 301: $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
! 302: *
! 303: * Reduce B to triangular form (QR decomposition of B)
! 304: * (Complex Workspace: need N, prefer N*NB)
! 305: *
! 306: IROWS = IHI + 1 - ILO
! 307: IF( ILV ) THEN
! 308: ICOLS = N + 1 - ILO
! 309: ELSE
! 310: ICOLS = IROWS
! 311: END IF
! 312: ITAU = 1
! 313: IWRK = ITAU + IROWS
! 314: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
! 315: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
! 316: *
! 317: * Apply the orthogonal transformation to matrix A
! 318: * (Complex Workspace: need N, prefer N*NB)
! 319: *
! 320: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
! 321: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
! 322: $ LWORK+1-IWRK, IERR )
! 323: *
! 324: * Initialize VL
! 325: * (Complex Workspace: need N, prefer N*NB)
! 326: *
! 327: IF( ILVL ) THEN
! 328: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
! 329: IF( IROWS.GT.1 ) THEN
! 330: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
! 331: $ VL( ILO+1, ILO ), LDVL )
! 332: END IF
! 333: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
! 334: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
! 335: END IF
! 336: *
! 337: * Initialize VR
! 338: *
! 339: IF( ILVR )
! 340: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
! 341: *
! 342: * Reduce to generalized Hessenberg form
! 343: *
! 344: IF( ILV ) THEN
! 345: *
! 346: * Eigenvectors requested -- work on whole matrix.
! 347: *
! 348: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
! 349: $ LDVL, VR, LDVR, IERR )
! 350: ELSE
! 351: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
! 352: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
! 353: END IF
! 354: *
! 355: * Perform QZ algorithm (Compute eigenvalues, and optionally, the
! 356: * Schur form and Schur vectors)
! 357: * (Complex Workspace: need N)
! 358: * (Real Workspace: need N)
! 359: *
! 360: IWRK = ITAU
! 361: IF( ILV ) THEN
! 362: CHTEMP = 'S'
! 363: ELSE
! 364: CHTEMP = 'E'
! 365: END IF
! 366: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
! 367: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
! 368: $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
! 369: IF( IERR.NE.0 ) THEN
! 370: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
! 371: INFO = IERR
! 372: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
! 373: INFO = IERR - N
! 374: ELSE
! 375: INFO = N + 1
! 376: END IF
! 377: GO TO 70
! 378: END IF
! 379: *
! 380: * Compute Eigenvectors
! 381: * (Real Workspace: need 2*N)
! 382: * (Complex Workspace: need 2*N)
! 383: *
! 384: IF( ILV ) THEN
! 385: IF( ILVL ) THEN
! 386: IF( ILVR ) THEN
! 387: CHTEMP = 'B'
! 388: ELSE
! 389: CHTEMP = 'L'
! 390: END IF
! 391: ELSE
! 392: CHTEMP = 'R'
! 393: END IF
! 394: *
! 395: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
! 396: $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
! 397: $ IERR )
! 398: IF( IERR.NE.0 ) THEN
! 399: INFO = N + 2
! 400: GO TO 70
! 401: END IF
! 402: *
! 403: * Undo balancing on VL and VR and normalization
! 404: * (Workspace: none needed)
! 405: *
! 406: IF( ILVL ) THEN
! 407: CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
! 408: $ RWORK( IRIGHT ), N, VL, LDVL, IERR )
! 409: DO 30 JC = 1, N
! 410: TEMP = ZERO
! 411: DO 10 JR = 1, N
! 412: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
! 413: 10 CONTINUE
! 414: IF( TEMP.LT.SMLNUM )
! 415: $ GO TO 30
! 416: TEMP = ONE / TEMP
! 417: DO 20 JR = 1, N
! 418: VL( JR, JC ) = VL( JR, JC )*TEMP
! 419: 20 CONTINUE
! 420: 30 CONTINUE
! 421: END IF
! 422: IF( ILVR ) THEN
! 423: CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
! 424: $ RWORK( IRIGHT ), N, VR, LDVR, IERR )
! 425: DO 60 JC = 1, N
! 426: TEMP = ZERO
! 427: DO 40 JR = 1, N
! 428: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
! 429: 40 CONTINUE
! 430: IF( TEMP.LT.SMLNUM )
! 431: $ GO TO 60
! 432: TEMP = ONE / TEMP
! 433: DO 50 JR = 1, N
! 434: VR( JR, JC ) = VR( JR, JC )*TEMP
! 435: 50 CONTINUE
! 436: 60 CONTINUE
! 437: END IF
! 438: END IF
! 439: *
! 440: * Undo scaling if necessary
! 441: *
! 442: IF( ILASCL )
! 443: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
! 444: *
! 445: IF( ILBSCL )
! 446: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
! 447: *
! 448: 70 CONTINUE
! 449: WORK( 1 ) = LWKOPT
! 450: *
! 451: RETURN
! 452: *
! 453: * End of ZGGEV
! 454: *
! 455: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggev.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack