Annotation of rpl/lapack/lapack/zgegv.f, revision 1.16

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

CVSweb interface <joel.bertrand@systella.fr>