Annotation of rpl/lapack/lapack/dggev.f, revision 1.18

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

CVSweb interface <joel.bertrand@systella.fr>