Annotation of rpl/lapack/lapack/zgegv.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGEGV( 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: * This routine is deprecated and has been replaced by routine ZGGEV.
! 24: *
! 25: * ZGEGV computes the eigenvalues and, optionally, the left and/or right
! 26: * eigenvectors of a complex matrix pair (A,B).
! 27: * Given two square matrices A and B,
! 28: * the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
! 29: * eigenvalues lambda and corresponding (non-zero) eigenvectors x such
! 30: * that
! 31: * A*x = lambda*B*x.
! 32: *
! 33: * An alternate form is to find the eigenvalues mu and corresponding
! 34: * eigenvectors y such that
! 35: * mu*A*y = B*y.
! 36: *
! 37: * These two forms are equivalent with mu = 1/lambda and x = y if
! 38: * neither lambda nor mu is zero. In order to deal with the case that
! 39: * lambda or mu is zero or small, two values alpha and beta are returned
! 40: * for each eigenvalue, such that lambda = alpha/beta and
! 41: * mu = beta/alpha.
! 42: *
! 43: * The vectors x and y in the above equations are right eigenvectors of
! 44: * the matrix pair (A,B). Vectors u and v satisfying
! 45: * u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
! 46: * are left eigenvectors of (A,B).
! 47: *
! 48: * Note: this routine performs "full balancing" on A and B -- see
! 49: * "Further Details", below.
! 50: *
! 51: * Arguments
! 52: * =========
! 53: *
! 54: * JOBVL (input) CHARACTER*1
! 55: * = 'N': do not compute the left generalized eigenvectors;
! 56: * = 'V': compute the left generalized eigenvectors (returned
! 57: * in VL).
! 58: *
! 59: * JOBVR (input) CHARACTER*1
! 60: * = 'N': do not compute the right generalized eigenvectors;
! 61: * = 'V': compute the right generalized eigenvectors (returned
! 62: * in VR).
! 63: *
! 64: * N (input) INTEGER
! 65: * The order of the matrices A, B, VL, and VR. N >= 0.
! 66: *
! 67: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
! 68: * On entry, the matrix A.
! 69: * If JOBVL = 'V' or JOBVR = 'V', then on exit A
! 70: * contains the Schur form of A from the generalized Schur
! 71: * factorization of the pair (A,B) after balancing. If no
! 72: * eigenvectors were computed, then only the diagonal elements
! 73: * of the Schur form will be correct. See ZGGHRD and ZHGEQZ
! 74: * for details.
! 75: *
! 76: * LDA (input) INTEGER
! 77: * The leading dimension of A. LDA >= max(1,N).
! 78: *
! 79: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
! 80: * On entry, the matrix B.
! 81: * If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
! 82: * upper triangular matrix obtained from B in the generalized
! 83: * Schur factorization of the pair (A,B) after balancing.
! 84: * If no eigenvectors were computed, then only the diagonal
! 85: * elements of B will be correct. See ZGGHRD and ZHGEQZ for
! 86: * details.
! 87: *
! 88: * LDB (input) INTEGER
! 89: * The leading dimension of B. LDB >= max(1,N).
! 90: *
! 91: * ALPHA (output) COMPLEX*16 array, dimension (N)
! 92: * The complex scalars alpha that define the eigenvalues of
! 93: * GNEP.
! 94: *
! 95: * BETA (output) COMPLEX*16 array, dimension (N)
! 96: * The complex scalars beta that define the eigenvalues of GNEP.
! 97: *
! 98: * Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
! 99: * represent the j-th eigenvalue of the matrix pair (A,B), in
! 100: * one of the forms lambda = alpha/beta or mu = beta/alpha.
! 101: * Since either lambda or mu may overflow, they should not,
! 102: * in general, be computed.
! 103: *
! 104: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
! 105: * If JOBVL = 'V', the left eigenvectors u(j) are stored
! 106: * in the columns of VL, in the same order as their eigenvalues.
! 107: * Each eigenvector is scaled so that its largest component has
! 108: * abs(real part) + abs(imag. part) = 1, except for eigenvectors
! 109: * corresponding to an eigenvalue with alpha = beta = 0, which
! 110: * are set to zero.
! 111: * Not referenced if JOBVL = 'N'.
! 112: *
! 113: * LDVL (input) INTEGER
! 114: * The leading dimension of the matrix VL. LDVL >= 1, and
! 115: * if JOBVL = 'V', LDVL >= N.
! 116: *
! 117: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
! 118: * If JOBVR = 'V', the right eigenvectors x(j) are stored
! 119: * in the columns of VR, in the same order as their eigenvalues.
! 120: * Each eigenvector is scaled so that its largest component has
! 121: * abs(real part) + abs(imag. part) = 1, except for eigenvectors
! 122: * corresponding to an eigenvalue with alpha = beta = 0, which
! 123: * are set to zero.
! 124: * Not referenced if JOBVR = 'N'.
! 125: *
! 126: * LDVR (input) INTEGER
! 127: * The leading dimension of the matrix VR. LDVR >= 1, and
! 128: * if JOBVR = 'V', LDVR >= N.
! 129: *
! 130: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 131: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 132: *
! 133: * LWORK (input) INTEGER
! 134: * The dimension of the array WORK. LWORK >= max(1,2*N).
! 135: * For good performance, LWORK must generally be larger.
! 136: * To compute the optimal value of LWORK, call ILAENV to get
! 137: * blocksizes (for ZGEQRF, ZUNMQR, and ZUNGQR.) Then compute:
! 138: * NB -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and ZUNGQR;
! 139: * The optimal LWORK is MAX( 2*N, N*(NB+1) ).
! 140: *
! 141: * If LWORK = -1, then a workspace query is assumed; the routine
! 142: * only calculates the optimal size of the WORK array, returns
! 143: * this value as the first entry of the WORK array, and no error
! 144: * message related to LWORK is issued by XERBLA.
! 145: *
! 146: * RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N)
! 147: *
! 148: * INFO (output) INTEGER
! 149: * = 0: successful exit
! 150: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 151: * =1,...,N:
! 152: * The QZ iteration failed. No eigenvectors have been
! 153: * calculated, but ALPHA(j) and BETA(j) should be
! 154: * correct for j=INFO+1,...,N.
! 155: * > N: errors that usually indicate LAPACK problems:
! 156: * =N+1: error return from ZGGBAL
! 157: * =N+2: error return from ZGEQRF
! 158: * =N+3: error return from ZUNMQR
! 159: * =N+4: error return from ZUNGQR
! 160: * =N+5: error return from ZGGHRD
! 161: * =N+6: error return from ZHGEQZ (other than failed
! 162: * iteration)
! 163: * =N+7: error return from ZTGEVC
! 164: * =N+8: error return from ZGGBAK (computing VL)
! 165: * =N+9: error return from ZGGBAK (computing VR)
! 166: * =N+10: error return from ZLASCL (various calls)
! 167: *
! 168: * Further Details
! 169: * ===============
! 170: *
! 171: * Balancing
! 172: * ---------
! 173: *
! 174: * This driver calls ZGGBAL to both permute and scale rows and columns
! 175: * of A and B. The permutations PL and PR are chosen so that PL*A*PR
! 176: * and PL*B*R will be upper triangular except for the diagonal blocks
! 177: * A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
! 178: * possible. The diagonal scaling matrices DL and DR are chosen so
! 179: * that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
! 180: * one (except for the elements that start out zero.)
! 181: *
! 182: * After the eigenvalues and eigenvectors of the balanced matrices
! 183: * have been computed, ZGGBAK transforms the eigenvectors back to what
! 184: * they would have been (in perfect arithmetic) if they had not been
! 185: * balanced.
! 186: *
! 187: * Contents of A and B on Exit
! 188: * -------- -- - --- - -- ----
! 189: *
! 190: * If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
! 191: * both), then on exit the arrays A and B will contain the complex Schur
! 192: * form[*] of the "balanced" versions of A and B. If no eigenvectors
! 193: * are computed, then only the diagonal blocks will be correct.
! 194: *
! 195: * [*] In other words, upper triangular form.
! 196: *
! 197: * =====================================================================
! 198: *
! 199: * .. Parameters ..
! 200: DOUBLE PRECISION ZERO, ONE
! 201: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 202: COMPLEX*16 CZERO, CONE
! 203: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
! 204: $ CONE = ( 1.0D0, 0.0D0 ) )
! 205: * ..
! 206: * .. Local Scalars ..
! 207: LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
! 208: CHARACTER CHTEMP
! 209: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
! 210: $ IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR,
! 211: $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
! 212: DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
! 213: $ BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
! 214: $ SALFAR, SBETA, SCALE, TEMP
! 215: COMPLEX*16 X
! 216: * ..
! 217: * .. Local Arrays ..
! 218: LOGICAL LDUMMA( 1 )
! 219: * ..
! 220: * .. External Subroutines ..
! 221: EXTERNAL XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
! 222: $ ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR, ZUNMQR
! 223: * ..
! 224: * .. External Functions ..
! 225: LOGICAL LSAME
! 226: INTEGER ILAENV
! 227: DOUBLE PRECISION DLAMCH, ZLANGE
! 228: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 229: * ..
! 230: * .. Intrinsic Functions ..
! 231: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, MAX
! 232: * ..
! 233: * .. Statement Functions ..
! 234: DOUBLE PRECISION ABS1
! 235: * ..
! 236: * .. Statement Function definitions ..
! 237: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
! 238: * ..
! 239: * .. Executable Statements ..
! 240: *
! 241: * Decode the input arguments
! 242: *
! 243: IF( LSAME( JOBVL, 'N' ) ) THEN
! 244: IJOBVL = 1
! 245: ILVL = .FALSE.
! 246: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
! 247: IJOBVL = 2
! 248: ILVL = .TRUE.
! 249: ELSE
! 250: IJOBVL = -1
! 251: ILVL = .FALSE.
! 252: END IF
! 253: *
! 254: IF( LSAME( JOBVR, 'N' ) ) THEN
! 255: IJOBVR = 1
! 256: ILVR = .FALSE.
! 257: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
! 258: IJOBVR = 2
! 259: ILVR = .TRUE.
! 260: ELSE
! 261: IJOBVR = -1
! 262: ILVR = .FALSE.
! 263: END IF
! 264: ILV = ILVL .OR. ILVR
! 265: *
! 266: * Test the input arguments
! 267: *
! 268: LWKMIN = MAX( 2*N, 1 )
! 269: LWKOPT = LWKMIN
! 270: WORK( 1 ) = LWKOPT
! 271: LQUERY = ( LWORK.EQ.-1 )
! 272: INFO = 0
! 273: IF( IJOBVL.LE.0 ) THEN
! 274: INFO = -1
! 275: ELSE IF( IJOBVR.LE.0 ) THEN
! 276: INFO = -2
! 277: ELSE IF( N.LT.0 ) THEN
! 278: INFO = -3
! 279: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 280: INFO = -5
! 281: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 282: INFO = -7
! 283: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
! 284: INFO = -11
! 285: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
! 286: INFO = -13
! 287: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
! 288: INFO = -15
! 289: END IF
! 290: *
! 291: IF( INFO.EQ.0 ) THEN
! 292: NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
! 293: NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
! 294: NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
! 295: NB = MAX( NB1, NB2, NB3 )
! 296: LOPT = MAX( 2*N, N*( NB+1 ) )
! 297: WORK( 1 ) = LOPT
! 298: END IF
! 299: *
! 300: IF( INFO.NE.0 ) THEN
! 301: CALL XERBLA( 'ZGEGV ', -INFO )
! 302: RETURN
! 303: ELSE IF( LQUERY ) THEN
! 304: RETURN
! 305: END IF
! 306: *
! 307: * Quick return if possible
! 308: *
! 309: IF( N.EQ.0 )
! 310: $ RETURN
! 311: *
! 312: * Get machine constants
! 313: *
! 314: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
! 315: SAFMIN = DLAMCH( 'S' )
! 316: SAFMIN = SAFMIN + SAFMIN
! 317: SAFMAX = ONE / SAFMIN
! 318: *
! 319: * Scale A
! 320: *
! 321: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
! 322: ANRM1 = ANRM
! 323: ANRM2 = ONE
! 324: IF( ANRM.LT.ONE ) THEN
! 325: IF( SAFMAX*ANRM.LT.ONE ) THEN
! 326: ANRM1 = SAFMIN
! 327: ANRM2 = SAFMAX*ANRM
! 328: END IF
! 329: END IF
! 330: *
! 331: IF( ANRM.GT.ZERO ) THEN
! 332: CALL ZLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
! 333: IF( IINFO.NE.0 ) THEN
! 334: INFO = N + 10
! 335: RETURN
! 336: END IF
! 337: END IF
! 338: *
! 339: * Scale B
! 340: *
! 341: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
! 342: BNRM1 = BNRM
! 343: BNRM2 = ONE
! 344: IF( BNRM.LT.ONE ) THEN
! 345: IF( SAFMAX*BNRM.LT.ONE ) THEN
! 346: BNRM1 = SAFMIN
! 347: BNRM2 = SAFMAX*BNRM
! 348: END IF
! 349: END IF
! 350: *
! 351: IF( BNRM.GT.ZERO ) THEN
! 352: CALL ZLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
! 353: IF( IINFO.NE.0 ) THEN
! 354: INFO = N + 10
! 355: RETURN
! 356: END IF
! 357: END IF
! 358: *
! 359: * Permute the matrix to make it more nearly triangular
! 360: * Also "balance" the matrix.
! 361: *
! 362: ILEFT = 1
! 363: IRIGHT = N + 1
! 364: IRWORK = IRIGHT + N
! 365: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
! 366: $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
! 367: IF( IINFO.NE.0 ) THEN
! 368: INFO = N + 1
! 369: GO TO 80
! 370: END IF
! 371: *
! 372: * Reduce B to triangular form, and initialize VL and/or VR
! 373: *
! 374: IROWS = IHI + 1 - ILO
! 375: IF( ILV ) THEN
! 376: ICOLS = N + 1 - ILO
! 377: ELSE
! 378: ICOLS = IROWS
! 379: END IF
! 380: ITAU = 1
! 381: IWORK = ITAU + IROWS
! 382: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
! 383: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
! 384: IF( IINFO.GE.0 )
! 385: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 386: IF( IINFO.NE.0 ) THEN
! 387: INFO = N + 2
! 388: GO TO 80
! 389: END IF
! 390: *
! 391: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
! 392: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
! 393: $ LWORK+1-IWORK, IINFO )
! 394: IF( IINFO.GE.0 )
! 395: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 396: IF( IINFO.NE.0 ) THEN
! 397: INFO = N + 3
! 398: GO TO 80
! 399: END IF
! 400: *
! 401: IF( ILVL ) THEN
! 402: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
! 403: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
! 404: $ VL( ILO+1, ILO ), LDVL )
! 405: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
! 406: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
! 407: $ IINFO )
! 408: IF( IINFO.GE.0 )
! 409: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 410: IF( IINFO.NE.0 ) THEN
! 411: INFO = N + 4
! 412: GO TO 80
! 413: END IF
! 414: END IF
! 415: *
! 416: IF( ILVR )
! 417: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
! 418: *
! 419: * Reduce to generalized Hessenberg form
! 420: *
! 421: IF( ILV ) THEN
! 422: *
! 423: * Eigenvectors requested -- work on whole matrix.
! 424: *
! 425: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
! 426: $ LDVL, VR, LDVR, IINFO )
! 427: ELSE
! 428: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
! 429: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
! 430: END IF
! 431: IF( IINFO.NE.0 ) THEN
! 432: INFO = N + 5
! 433: GO TO 80
! 434: END IF
! 435: *
! 436: * Perform QZ algorithm
! 437: *
! 438: IWORK = ITAU
! 439: IF( ILV ) THEN
! 440: CHTEMP = 'S'
! 441: ELSE
! 442: CHTEMP = 'E'
! 443: END IF
! 444: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
! 445: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ),
! 446: $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
! 447: IF( IINFO.GE.0 )
! 448: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 449: IF( IINFO.NE.0 ) THEN
! 450: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
! 451: INFO = IINFO
! 452: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
! 453: INFO = IINFO - N
! 454: ELSE
! 455: INFO = N + 6
! 456: END IF
! 457: GO TO 80
! 458: END IF
! 459: *
! 460: IF( ILV ) THEN
! 461: *
! 462: * Compute Eigenvectors
! 463: *
! 464: IF( ILVL ) THEN
! 465: IF( ILVR ) THEN
! 466: CHTEMP = 'B'
! 467: ELSE
! 468: CHTEMP = 'L'
! 469: END IF
! 470: ELSE
! 471: CHTEMP = 'R'
! 472: END IF
! 473: *
! 474: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
! 475: $ VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ),
! 476: $ IINFO )
! 477: IF( IINFO.NE.0 ) THEN
! 478: INFO = N + 7
! 479: GO TO 80
! 480: END IF
! 481: *
! 482: * Undo balancing on VL and VR, rescale
! 483: *
! 484: IF( ILVL ) THEN
! 485: CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
! 486: $ RWORK( IRIGHT ), N, VL, LDVL, IINFO )
! 487: IF( IINFO.NE.0 ) THEN
! 488: INFO = N + 8
! 489: GO TO 80
! 490: END IF
! 491: DO 30 JC = 1, N
! 492: TEMP = ZERO
! 493: DO 10 JR = 1, N
! 494: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
! 495: 10 CONTINUE
! 496: IF( TEMP.LT.SAFMIN )
! 497: $ GO TO 30
! 498: TEMP = ONE / TEMP
! 499: DO 20 JR = 1, N
! 500: VL( JR, JC ) = VL( JR, JC )*TEMP
! 501: 20 CONTINUE
! 502: 30 CONTINUE
! 503: END IF
! 504: IF( ILVR ) THEN
! 505: CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
! 506: $ RWORK( IRIGHT ), N, VR, LDVR, IINFO )
! 507: IF( IINFO.NE.0 ) THEN
! 508: INFO = N + 9
! 509: GO TO 80
! 510: END IF
! 511: DO 60 JC = 1, N
! 512: TEMP = ZERO
! 513: DO 40 JR = 1, N
! 514: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
! 515: 40 CONTINUE
! 516: IF( TEMP.LT.SAFMIN )
! 517: $ GO TO 60
! 518: TEMP = ONE / TEMP
! 519: DO 50 JR = 1, N
! 520: VR( JR, JC ) = VR( JR, JC )*TEMP
! 521: 50 CONTINUE
! 522: 60 CONTINUE
! 523: END IF
! 524: *
! 525: * End of eigenvector calculation
! 526: *
! 527: END IF
! 528: *
! 529: * Undo scaling in alpha, beta
! 530: *
! 531: * Note: this does not give the alpha and beta for the unscaled
! 532: * problem.
! 533: *
! 534: * Un-scaling is limited to avoid underflow in alpha and beta
! 535: * if they are significant.
! 536: *
! 537: DO 70 JC = 1, N
! 538: ABSAR = ABS( DBLE( ALPHA( JC ) ) )
! 539: ABSAI = ABS( DIMAG( ALPHA( JC ) ) )
! 540: ABSB = ABS( DBLE( BETA( JC ) ) )
! 541: SALFAR = ANRM*DBLE( ALPHA( JC ) )
! 542: SALFAI = ANRM*DIMAG( ALPHA( JC ) )
! 543: SBETA = BNRM*DBLE( BETA( JC ) )
! 544: ILIMIT = .FALSE.
! 545: SCALE = ONE
! 546: *
! 547: * Check for significant underflow in imaginary part of ALPHA
! 548: *
! 549: IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
! 550: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
! 551: ILIMIT = .TRUE.
! 552: SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI )
! 553: END IF
! 554: *
! 555: * Check for significant underflow in real part of ALPHA
! 556: *
! 557: IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
! 558: $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
! 559: ILIMIT = .TRUE.
! 560: SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) /
! 561: $ MAX( SAFMIN, ANRM2*ABSAR ) )
! 562: END IF
! 563: *
! 564: * Check for significant underflow in BETA
! 565: *
! 566: IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
! 567: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
! 568: ILIMIT = .TRUE.
! 569: SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) /
! 570: $ MAX( SAFMIN, BNRM2*ABSB ) )
! 571: END IF
! 572: *
! 573: * Check for possible overflow when limiting scaling
! 574: *
! 575: IF( ILIMIT ) THEN
! 576: TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
! 577: $ ABS( SBETA ) )
! 578: IF( TEMP.GT.ONE )
! 579: $ SCALE = SCALE / TEMP
! 580: IF( SCALE.LT.ONE )
! 581: $ ILIMIT = .FALSE.
! 582: END IF
! 583: *
! 584: * Recompute un-scaled ALPHA, BETA if necessary.
! 585: *
! 586: IF( ILIMIT ) THEN
! 587: SALFAR = ( SCALE*DBLE( ALPHA( JC ) ) )*ANRM
! 588: SALFAI = ( SCALE*DIMAG( ALPHA( JC ) ) )*ANRM
! 589: SBETA = ( SCALE*BETA( JC ) )*BNRM
! 590: END IF
! 591: ALPHA( JC ) = DCMPLX( SALFAR, SALFAI )
! 592: BETA( JC ) = SBETA
! 593: 70 CONTINUE
! 594: *
! 595: 80 CONTINUE
! 596: WORK( 1 ) = LWKOPT
! 597: *
! 598: RETURN
! 599: *
! 600: * End of ZGEGV
! 601: *
! 602: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgegv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack