Annotation of rpl/lapack/lapack/dggev3.f, revision 1.3

1.1       bertrand    1: *> \brief <b> DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)</b>
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
                      7: *
                      8: *> \htmlonly
                      9: *> Download DGGEV3 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggev3.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggev3.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggev3.f">
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
                     22: *      $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
                     23: *      $                   INFO )
                     24: *
                     25: *       .. Scalar Arguments ..
                     26: *       CHARACTER          JOBVL, JOBVR
                     27: *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
                     28: *       ..
                     29: *       .. Array Arguments ..
                     30: *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
                     31: *      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
                     32: *      $                   VR( LDVR, * ), WORK( * )
                     33: *       ..
                     34: *
                     35: *
                     36: *> \par Purpose:
                     37: *  =============
                     38: *>
                     39: *> \verbatim
                     40: *>
                     41: *> DGGEV3 computes for a pair of N-by-N real nonsymmetric matrices (A,B)
                     42: *> the generalized eigenvalues, and optionally, the left and/or right
                     43: *> generalized eigenvectors.
                     44: *>
                     45: *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
                     46: *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
                     47: *> singular. It is usually represented as the pair (alpha,beta), as
                     48: *> there is a reasonable interpretation for beta=0, and even for both
                     49: *> being zero.
                     50: *>
                     51: *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
                     52: *> of (A,B) satisfies
                     53: *>
                     54: *>                  A * v(j) = lambda(j) * B * v(j).
                     55: *>
                     56: *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
                     57: *> of (A,B) satisfies
                     58: *>
                     59: *>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
                     60: *>
                     61: *> where u(j)**H is the conjugate-transpose of u(j).
                     62: *>
                     63: *> \endverbatim
                     64: *
                     65: *  Arguments:
                     66: *  ==========
                     67: *
                     68: *> \param[in] JOBVL
                     69: *> \verbatim
                     70: *>          JOBVL is CHARACTER*1
                     71: *>          = 'N':  do not compute the left generalized eigenvectors;
                     72: *>          = 'V':  compute the left generalized eigenvectors.
                     73: *> \endverbatim
                     74: *>
                     75: *> \param[in] JOBVR
                     76: *> \verbatim
                     77: *>          JOBVR is CHARACTER*1
                     78: *>          = 'N':  do not compute the right generalized eigenvectors;
                     79: *>          = 'V':  compute the right generalized eigenvectors.
                     80: *> \endverbatim
                     81: *>
                     82: *> \param[in] N
                     83: *> \verbatim
                     84: *>          N is INTEGER
                     85: *>          The order of the matrices A, B, VL, and VR.  N >= 0.
                     86: *> \endverbatim
                     87: *>
                     88: *> \param[in,out] A
                     89: *> \verbatim
                     90: *>          A is DOUBLE PRECISION array, dimension (LDA, N)
                     91: *>          On entry, the matrix A in the pair (A,B).
                     92: *>          On exit, A has been overwritten.
                     93: *> \endverbatim
                     94: *>
                     95: *> \param[in] LDA
                     96: *> \verbatim
                     97: *>          LDA is INTEGER
                     98: *>          The leading dimension of A.  LDA >= max(1,N).
                     99: *> \endverbatim
                    100: *>
                    101: *> \param[in,out] B
                    102: *> \verbatim
                    103: *>          B is DOUBLE PRECISION array, dimension (LDB, N)
                    104: *>          On entry, the matrix B in the pair (A,B).
                    105: *>          On exit, B has been overwritten.
                    106: *> \endverbatim
                    107: *>
                    108: *> \param[in] LDB
                    109: *> \verbatim
                    110: *>          LDB is INTEGER
                    111: *>          The leading dimension of B.  LDB >= max(1,N).
                    112: *> \endverbatim
                    113: *>
                    114: *> \param[out] ALPHAR
                    115: *> \verbatim
                    116: *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
                    117: *> \endverbatim
                    118: *>
                    119: *> \param[out] ALPHAI
                    120: *> \verbatim
                    121: *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
                    122: *> \endverbatim
                    123: *>
                    124: *> \param[out] BETA
                    125: *> \verbatim
                    126: *>          BETA is DOUBLE PRECISION array, dimension (N)
                    127: *>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
                    128: *>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
                    129: *>          the j-th eigenvalue is real; if positive, then the j-th and
                    130: *>          (j+1)-st eigenvalues are a complex conjugate pair, with
                    131: *>          ALPHAI(j+1) negative.
                    132: *>
                    133: *>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
                    134: *>          may easily over- or underflow, and BETA(j) may even be zero.
                    135: *>          Thus, the user should avoid naively computing the ratio
                    136: *>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
                    137: *>          than and usually comparable with norm(A) in magnitude, and
                    138: *>          BETA always less than and usually comparable with norm(B).
                    139: *> \endverbatim
                    140: *>
                    141: *> \param[out] VL
                    142: *> \verbatim
                    143: *>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
                    144: *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
                    145: *>          after another in the columns of VL, in the same order as
                    146: *>          their eigenvalues. If the j-th eigenvalue is real, then
                    147: *>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
                    148: *>          (j+1)-th eigenvalues form a complex conjugate pair, then
                    149: *>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
                    150: *>          Each eigenvector is scaled so the largest component has
                    151: *>          abs(real part)+abs(imag. part)=1.
                    152: *>          Not referenced if JOBVL = 'N'.
                    153: *> \endverbatim
                    154: *>
                    155: *> \param[in] LDVL
                    156: *> \verbatim
                    157: *>          LDVL is INTEGER
                    158: *>          The leading dimension of the matrix VL. LDVL >= 1, and
                    159: *>          if JOBVL = 'V', LDVL >= N.
                    160: *> \endverbatim
                    161: *>
                    162: *> \param[out] VR
                    163: *> \verbatim
                    164: *>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
                    165: *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
                    166: *>          after another in the columns of VR, in the same order as
                    167: *>          their eigenvalues. If the j-th eigenvalue is real, then
                    168: *>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
                    169: *>          (j+1)-th eigenvalues form a complex conjugate pair, then
                    170: *>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
                    171: *>          Each eigenvector is scaled so the largest component has
                    172: *>          abs(real part)+abs(imag. part)=1.
                    173: *>          Not referenced if JOBVR = 'N'.
                    174: *> \endverbatim
                    175: *>
                    176: *> \param[in] LDVR
                    177: *> \verbatim
                    178: *>          LDVR is INTEGER
                    179: *>          The leading dimension of the matrix VR. LDVR >= 1, and
                    180: *>          if JOBVR = 'V', LDVR >= N.
                    181: *> \endverbatim
                    182: *>
                    183: *> \param[out] WORK
                    184: *> \verbatim
                    185: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                    186: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                    187: *> \endverbatim
                    188: *>
                    189: *> \param[in] LWORK
                    190: *> \verbatim
                    191: *>          LWORK is INTEGER
                    192: *>
                    193: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    194: *>          only calculates the optimal size of the WORK array, returns
                    195: *>          this value as the first entry of the WORK array, and no error
                    196: *>          message related to LWORK is issued by XERBLA.
                    197: *> \endverbatim
                    198: *>
                    199: *> \param[out] INFO
                    200: *> \verbatim
                    201: *>          INFO is INTEGER
                    202: *>          = 0:  successful exit
                    203: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    204: *>          = 1,...,N:
                    205: *>                The QZ iteration failed.  No eigenvectors have been
                    206: *>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
                    207: *>                should be correct for j=INFO+1,...,N.
                    208: *>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
                    209: *>                =N+2: error return from DTGEVC.
                    210: *> \endverbatim
                    211: *
                    212: *  Authors:
                    213: *  ========
                    214: *
                    215: *> \author Univ. of Tennessee
                    216: *> \author Univ. of California Berkeley
                    217: *> \author Univ. of Colorado Denver
                    218: *> \author NAG Ltd.
                    219: *
                    220: *> \date January 2015
                    221: *
                    222: *> \ingroup doubleGEeigen
                    223: *
                    224: *  =====================================================================
                    225:       SUBROUTINE DGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
                    226:      $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
                    227:      $                   INFO )
                    228: *
                    229: *  -- LAPACK driver routine (version 3.6.0) --
                    230: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    231: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    232: *     January 2015
                    233: *
                    234: *     .. Scalar Arguments ..
                    235:       CHARACTER          JOBVL, JOBVR
                    236:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
                    237: *     ..
                    238: *     .. Array Arguments ..
                    239:       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
                    240:      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
                    241:      $                   VR( LDVR, * ), WORK( * )
                    242: *     ..
                    243: *
                    244: *  =====================================================================
                    245: *
                    246: *     .. Parameters ..
                    247:       DOUBLE PRECISION   ZERO, ONE
                    248:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
                    249: *     ..
                    250: *     .. Local Scalars ..
                    251:       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
                    252:       CHARACTER          CHTEMP
                    253:       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
                    254:      $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, LWKOPT
                    255:       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
                    256:      $                   SMLNUM, TEMP
                    257: *     ..
                    258: *     .. Local Arrays ..
                    259:       LOGICAL            LDUMMA( 1 )
                    260: *     ..
                    261: *     .. External Subroutines ..
                    262:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHD3, DHGEQZ, DLABAD,
                    263:      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
                    264:      $                   XERBLA
                    265: *     ..
                    266: *     .. External Functions ..
                    267:       LOGICAL            LSAME
                    268:       DOUBLE PRECISION   DLAMCH, DLANGE
                    269:       EXTERNAL           LSAME, DLAMCH, DLANGE
                    270: *     ..
                    271: *     .. Intrinsic Functions ..
                    272:       INTRINSIC          ABS, MAX, SQRT
                    273: *     ..
                    274: *     .. Executable Statements ..
                    275: *
                    276: *     Decode the input arguments
                    277: *
                    278:       IF( LSAME( JOBVL, 'N' ) ) THEN
                    279:          IJOBVL = 1
                    280:          ILVL = .FALSE.
                    281:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
                    282:          IJOBVL = 2
                    283:          ILVL = .TRUE.
                    284:       ELSE
                    285:          IJOBVL = -1
                    286:          ILVL = .FALSE.
                    287:       END IF
                    288: *
                    289:       IF( LSAME( JOBVR, 'N' ) ) THEN
                    290:          IJOBVR = 1
                    291:          ILVR = .FALSE.
                    292:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
                    293:          IJOBVR = 2
                    294:          ILVR = .TRUE.
                    295:       ELSE
                    296:          IJOBVR = -1
                    297:          ILVR = .FALSE.
                    298:       END IF
                    299:       ILV = ILVL .OR. ILVR
                    300: *
                    301: *     Test the input arguments
                    302: *
                    303:       INFO = 0
                    304:       LQUERY = ( LWORK.EQ.-1 )
                    305:       IF( IJOBVL.LE.0 ) THEN
                    306:          INFO = -1
                    307:       ELSE IF( IJOBVR.LE.0 ) THEN
                    308:          INFO = -2
                    309:       ELSE IF( N.LT.0 ) THEN
                    310:          INFO = -3
                    311:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    312:          INFO = -5
                    313:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    314:          INFO = -7
                    315:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
                    316:          INFO = -12
                    317:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
                    318:          INFO = -14
                    319:       ELSE IF( LWORK.LT.MAX( 1, 8*N ) .AND. .NOT.LQUERY ) THEN
                    320:          INFO = -16
                    321:       END IF
                    322: *
                    323: *     Compute workspace
                    324: *
                    325:       IF( INFO.EQ.0 ) THEN
                    326:          CALL DGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
                    327:          LWKOPT = MAX(1, 8*N, 3*N+INT( WORK( 1 ) ) )
                    328:          CALL DORMQR( 'L', 'T', N, N, N, B, LDB, WORK, A, LDA, WORK, -1,
                    329:      $                IERR )
                    330:          LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
                    331:          IF( ILVL ) THEN
                    332:             CALL DORGQR( N, N, N, VL, LDVL, WORK, WORK, -1, IERR )
                    333:             LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
                    334:          END IF
                    335:          IF( ILV ) THEN
                    336:             CALL DGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
                    337:      $                   LDVL, VR, LDVR, WORK, -1, IERR )
                    338:             LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
                    339:             CALL DHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
                    340:      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
                    341:      $                   WORK, -1, IERR )
                    342:             LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
                    343:          ELSE
                    344:             CALL DGGHD3( 'N', 'N', N, 1, N, A, LDA, B, LDB, VL, LDVL,
                    345:      $                   VR, LDVR, WORK, -1, IERR )
                    346:             LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
                    347:             CALL DHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
                    348:      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
                    349:      $                   WORK, -1, IERR )
                    350:             LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
                    351:          END IF
                    352: 
                    353:          WORK( 1 ) = LWKOPT
                    354:       END IF
                    355: *
                    356:       IF( INFO.NE.0 ) THEN
                    357:          CALL XERBLA( 'DGGEV3 ', -INFO )
                    358:          RETURN
                    359:       ELSE IF( LQUERY ) THEN
                    360:          RETURN
                    361:       END IF
                    362: *
                    363: *     Quick return if possible
                    364: *
                    365:       IF( N.EQ.0 )
                    366:      $   RETURN
                    367: *
                    368: *     Get machine constants
                    369: *
                    370:       EPS = DLAMCH( 'P' )
                    371:       SMLNUM = DLAMCH( 'S' )
                    372:       BIGNUM = ONE / SMLNUM
                    373:       CALL DLABAD( SMLNUM, BIGNUM )
                    374:       SMLNUM = SQRT( SMLNUM ) / EPS
                    375:       BIGNUM = ONE / SMLNUM
                    376: *
                    377: *     Scale A if max element outside range [SMLNUM,BIGNUM]
                    378: *
                    379:       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
                    380:       ILASCL = .FALSE.
                    381:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
                    382:          ANRMTO = SMLNUM
                    383:          ILASCL = .TRUE.
                    384:       ELSE IF( ANRM.GT.BIGNUM ) THEN
                    385:          ANRMTO = BIGNUM
                    386:          ILASCL = .TRUE.
                    387:       END IF
                    388:       IF( ILASCL )
                    389:      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
                    390: *
                    391: *     Scale B if max element outside range [SMLNUM,BIGNUM]
                    392: *
                    393:       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
                    394:       ILBSCL = .FALSE.
                    395:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
                    396:          BNRMTO = SMLNUM
                    397:          ILBSCL = .TRUE.
                    398:       ELSE IF( BNRM.GT.BIGNUM ) THEN
                    399:          BNRMTO = BIGNUM
                    400:          ILBSCL = .TRUE.
                    401:       END IF
                    402:       IF( ILBSCL )
                    403:      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
                    404: *
                    405: *     Permute the matrices A, B to isolate eigenvalues if possible
                    406: *
                    407:       ILEFT = 1
                    408:       IRIGHT = N + 1
                    409:       IWRK = IRIGHT + N
                    410:       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
                    411:      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
                    412: *
                    413: *     Reduce B to triangular form (QR decomposition of B)
                    414: *
                    415:       IROWS = IHI + 1 - ILO
                    416:       IF( ILV ) THEN
                    417:          ICOLS = N + 1 - ILO
                    418:       ELSE
                    419:          ICOLS = IROWS
                    420:       END IF
                    421:       ITAU = IWRK
                    422:       IWRK = ITAU + IROWS
                    423:       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
                    424:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
                    425: *
                    426: *     Apply the orthogonal transformation to matrix A
                    427: *
                    428:       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
                    429:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
                    430:      $             LWORK+1-IWRK, IERR )
                    431: *
                    432: *     Initialize VL
                    433: *
                    434:       IF( ILVL ) THEN
                    435:          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
                    436:          IF( IROWS.GT.1 ) THEN
                    437:             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
                    438:      $                   VL( ILO+1, ILO ), LDVL )
                    439:          END IF
                    440:          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
                    441:      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
                    442:       END IF
                    443: *
                    444: *     Initialize VR
                    445: *
                    446:       IF( ILVR )
                    447:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
                    448: *
                    449: *     Reduce to generalized Hessenberg form
                    450: *
                    451:       IF( ILV ) THEN
                    452: *
                    453: *        Eigenvectors requested -- work on whole matrix.
                    454: *
                    455:          CALL DGGHD3( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
                    456:      $                LDVL, VR, LDVR, WORK( IWRK ), LWORK+1-IWRK, IERR )
                    457:       ELSE
                    458:          CALL DGGHD3( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
                    459:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR,
                    460:      $                WORK( IWRK ), LWORK+1-IWRK, IERR )
                    461:       END IF
                    462: *
                    463: *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
                    464: *     Schur forms and Schur vectors)
                    465: *
                    466:       IWRK = ITAU
                    467:       IF( ILV ) THEN
                    468:          CHTEMP = 'S'
                    469:       ELSE
                    470:          CHTEMP = 'E'
                    471:       END IF
                    472:       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
                    473:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
                    474:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
                    475:       IF( IERR.NE.0 ) THEN
                    476:          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
                    477:             INFO = IERR
                    478:          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
                    479:             INFO = IERR - N
                    480:          ELSE
                    481:             INFO = N + 1
                    482:          END IF
                    483:          GO TO 110
                    484:       END IF
                    485: *
                    486: *     Compute Eigenvectors
                    487: *
                    488:       IF( ILV ) THEN
                    489:          IF( ILVL ) THEN
                    490:             IF( ILVR ) THEN
                    491:                CHTEMP = 'B'
                    492:             ELSE
                    493:                CHTEMP = 'L'
                    494:             END IF
                    495:          ELSE
                    496:             CHTEMP = 'R'
                    497:          END IF
                    498:          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
                    499:      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
                    500:          IF( IERR.NE.0 ) THEN
                    501:             INFO = N + 2
                    502:             GO TO 110
                    503:          END IF
                    504: *
                    505: *        Undo balancing on VL and VR and normalization
                    506: *
                    507:          IF( ILVL ) THEN
                    508:             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
                    509:      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
                    510:             DO 50 JC = 1, N
                    511:                IF( ALPHAI( JC ).LT.ZERO )
                    512:      $            GO TO 50
                    513:                TEMP = ZERO
                    514:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    515:                   DO 10 JR = 1, N
                    516:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
                    517:    10             CONTINUE
                    518:                ELSE
                    519:                   DO 20 JR = 1, N
                    520:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
                    521:      $                      ABS( VL( JR, JC+1 ) ) )
                    522:    20             CONTINUE
                    523:                END IF
                    524:                IF( TEMP.LT.SMLNUM )
                    525:      $            GO TO 50
                    526:                TEMP = ONE / TEMP
                    527:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    528:                   DO 30 JR = 1, N
                    529:                      VL( JR, JC ) = VL( JR, JC )*TEMP
                    530:    30             CONTINUE
                    531:                ELSE
                    532:                   DO 40 JR = 1, N
                    533:                      VL( JR, JC ) = VL( JR, JC )*TEMP
                    534:                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
                    535:    40             CONTINUE
                    536:                END IF
                    537:    50       CONTINUE
                    538:          END IF
                    539:          IF( ILVR ) THEN
                    540:             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
                    541:      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
                    542:             DO 100 JC = 1, N
                    543:                IF( ALPHAI( JC ).LT.ZERO )
                    544:      $            GO TO 100
                    545:                TEMP = ZERO
                    546:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    547:                   DO 60 JR = 1, N
                    548:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
                    549:    60             CONTINUE
                    550:                ELSE
                    551:                   DO 70 JR = 1, N
                    552:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
                    553:      $                      ABS( VR( JR, JC+1 ) ) )
                    554:    70             CONTINUE
                    555:                END IF
                    556:                IF( TEMP.LT.SMLNUM )
                    557:      $            GO TO 100
                    558:                TEMP = ONE / TEMP
                    559:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    560:                   DO 80 JR = 1, N
                    561:                      VR( JR, JC ) = VR( JR, JC )*TEMP
                    562:    80             CONTINUE
                    563:                ELSE
                    564:                   DO 90 JR = 1, N
                    565:                      VR( JR, JC ) = VR( JR, JC )*TEMP
                    566:                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
                    567:    90             CONTINUE
                    568:                END IF
                    569:   100       CONTINUE
                    570:          END IF
                    571: *
                    572: *        End of eigenvector calculation
                    573: *
                    574:       END IF
                    575: *
                    576: *     Undo scaling if necessary
                    577: *
                    578:   110 CONTINUE
                    579: *
                    580:       IF( ILASCL ) THEN
                    581:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
                    582:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
                    583:       END IF
                    584: *
                    585:       IF( ILBSCL ) THEN
                    586:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
                    587:       END IF
                    588: *
                    589:       WORK( 1 ) = LWKOPT
                    590:       RETURN
                    591: *
                    592: *     End of DGGEV3
                    593: *
                    594:       END

CVSweb interface <joel.bertrand@systella.fr>