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

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: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      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"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     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 )
                     23: * 
                     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: *       ..
                     33: *  
                     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: *
                    216: *> \author Univ. of Tennessee 
                    217: *> \author Univ. of California Berkeley 
                    218: *> \author Univ. of Colorado Denver 
                    219: *> \author NAG Ltd. 
                    220: *
1.10      bertrand  221: *> \date April 2012
1.8       bertrand  222: *
                    223: *> \ingroup doubleGEeigen
                    224: *
                    225: *  =====================================================================
1.1       bertrand  226:       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
                    227:      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
                    228: *
1.10      bertrand  229: *  -- LAPACK driver routine (version 3.4.1) --
1.1       bertrand  230: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    231: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10      bertrand  232: *     April 2012
1.1       bertrand  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, MAXWRK,
                    255:      $                   MINWRK
                    256:       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
                    257:      $                   SMLNUM, TEMP
                    258: *     ..
                    259: *     .. Local Arrays ..
                    260:       LOGICAL            LDUMMA( 1 )
                    261: *     ..
                    262: *     .. External Subroutines ..
                    263:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
                    264:      $                   DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
                    265:      $                   XERBLA
                    266: *     ..
                    267: *     .. External Functions ..
                    268:       LOGICAL            LSAME
                    269:       INTEGER            ILAENV
                    270:       DOUBLE PRECISION   DLAMCH, DLANGE
                    271:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
                    272: *     ..
                    273: *     .. Intrinsic Functions ..
                    274:       INTRINSIC          ABS, MAX, SQRT
                    275: *     ..
                    276: *     .. Executable Statements ..
                    277: *
                    278: *     Decode the input arguments
                    279: *
                    280:       IF( LSAME( JOBVL, 'N' ) ) THEN
                    281:          IJOBVL = 1
                    282:          ILVL = .FALSE.
                    283:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
                    284:          IJOBVL = 2
                    285:          ILVL = .TRUE.
                    286:       ELSE
                    287:          IJOBVL = -1
                    288:          ILVL = .FALSE.
                    289:       END IF
                    290: *
                    291:       IF( LSAME( JOBVR, 'N' ) ) THEN
                    292:          IJOBVR = 1
                    293:          ILVR = .FALSE.
                    294:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
                    295:          IJOBVR = 2
                    296:          ILVR = .TRUE.
                    297:       ELSE
                    298:          IJOBVR = -1
                    299:          ILVR = .FALSE.
                    300:       END IF
                    301:       ILV = ILVL .OR. ILVR
                    302: *
                    303: *     Test the input arguments
                    304: *
                    305:       INFO = 0
                    306:       LQUERY = ( LWORK.EQ.-1 )
                    307:       IF( IJOBVL.LE.0 ) THEN
                    308:          INFO = -1
                    309:       ELSE IF( IJOBVR.LE.0 ) THEN
                    310:          INFO = -2
                    311:       ELSE IF( N.LT.0 ) THEN
                    312:          INFO = -3
                    313:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    314:          INFO = -5
                    315:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    316:          INFO = -7
                    317:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
                    318:          INFO = -12
                    319:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
                    320:          INFO = -14
                    321:       END IF
                    322: *
                    323: *     Compute workspace
                    324: *      (Note: Comments in the code beginning "Workspace:" describe the
                    325: *       minimal amount of workspace needed at that point in the code,
                    326: *       as well as the preferred amount for good performance.
                    327: *       NB refers to the optimal block size for the immediately
                    328: *       following subroutine, as returned by ILAENV. The workspace is
                    329: *       computed assuming ILO = 1 and IHI = N, the worst case.)
                    330: *
                    331:       IF( INFO.EQ.0 ) THEN
                    332:          MINWRK = MAX( 1, 8*N )
                    333:          MAXWRK = MAX( 1, N*( 7 +
                    334:      $                 ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) )
                    335:          MAXWRK = MAX( MAXWRK, N*( 7 +
                    336:      $                 ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) )
                    337:          IF( ILVL ) THEN
                    338:             MAXWRK = MAX( MAXWRK, N*( 7 +
                    339:      $                 ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) )
                    340:          END IF
                    341:          WORK( 1 ) = MAXWRK
                    342: *
                    343:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
                    344:      $      INFO = -16
                    345:       END IF
                    346: *
                    347:       IF( INFO.NE.0 ) THEN
                    348:          CALL XERBLA( 'DGGEV ', -INFO )
                    349:          RETURN
                    350:       ELSE IF( LQUERY ) THEN
                    351:          RETURN
                    352:       END IF
                    353: *
                    354: *     Quick return if possible
                    355: *
                    356:       IF( N.EQ.0 )
                    357:      $   RETURN
                    358: *
                    359: *     Get machine constants
                    360: *
                    361:       EPS = DLAMCH( 'P' )
                    362:       SMLNUM = DLAMCH( 'S' )
                    363:       BIGNUM = ONE / SMLNUM
                    364:       CALL DLABAD( SMLNUM, BIGNUM )
                    365:       SMLNUM = SQRT( SMLNUM ) / EPS
                    366:       BIGNUM = ONE / SMLNUM
                    367: *
                    368: *     Scale A if max element outside range [SMLNUM,BIGNUM]
                    369: *
                    370:       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
                    371:       ILASCL = .FALSE.
                    372:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
                    373:          ANRMTO = SMLNUM
                    374:          ILASCL = .TRUE.
                    375:       ELSE IF( ANRM.GT.BIGNUM ) THEN
                    376:          ANRMTO = BIGNUM
                    377:          ILASCL = .TRUE.
                    378:       END IF
                    379:       IF( ILASCL )
                    380:      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
                    381: *
                    382: *     Scale B if max element outside range [SMLNUM,BIGNUM]
                    383: *
                    384:       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
                    385:       ILBSCL = .FALSE.
                    386:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
                    387:          BNRMTO = SMLNUM
                    388:          ILBSCL = .TRUE.
                    389:       ELSE IF( BNRM.GT.BIGNUM ) THEN
                    390:          BNRMTO = BIGNUM
                    391:          ILBSCL = .TRUE.
                    392:       END IF
                    393:       IF( ILBSCL )
                    394:      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
                    395: *
                    396: *     Permute the matrices A, B to isolate eigenvalues if possible
                    397: *     (Workspace: need 6*N)
                    398: *
                    399:       ILEFT = 1
                    400:       IRIGHT = N + 1
                    401:       IWRK = IRIGHT + N
                    402:       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
                    403:      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
                    404: *
                    405: *     Reduce B to triangular form (QR decomposition of B)
                    406: *     (Workspace: need N, prefer N*NB)
                    407: *
                    408:       IROWS = IHI + 1 - ILO
                    409:       IF( ILV ) THEN
                    410:          ICOLS = N + 1 - ILO
                    411:       ELSE
                    412:          ICOLS = IROWS
                    413:       END IF
                    414:       ITAU = IWRK
                    415:       IWRK = ITAU + IROWS
                    416:       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
                    417:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
                    418: *
                    419: *     Apply the orthogonal transformation to matrix A
                    420: *     (Workspace: need N, prefer N*NB)
                    421: *
                    422:       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
                    423:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
                    424:      $             LWORK+1-IWRK, IERR )
                    425: *
                    426: *     Initialize VL
                    427: *     (Workspace: need N, prefer N*NB)
                    428: *
                    429:       IF( ILVL ) THEN
                    430:          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
                    431:          IF( IROWS.GT.1 ) THEN
                    432:             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
                    433:      $                   VL( ILO+1, ILO ), LDVL )
                    434:          END IF
                    435:          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
                    436:      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
                    437:       END IF
                    438: *
                    439: *     Initialize VR
                    440: *
                    441:       IF( ILVR )
                    442:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
                    443: *
                    444: *     Reduce to generalized Hessenberg form
                    445: *     (Workspace: none needed)
                    446: *
                    447:       IF( ILV ) THEN
                    448: *
                    449: *        Eigenvectors requested -- work on whole matrix.
                    450: *
                    451:          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
                    452:      $                LDVL, VR, LDVR, IERR )
                    453:       ELSE
                    454:          CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
                    455:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
                    456:       END IF
                    457: *
                    458: *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
                    459: *     Schur forms and Schur vectors)
                    460: *     (Workspace: need N)
                    461: *
                    462:       IWRK = ITAU
                    463:       IF( ILV ) THEN
                    464:          CHTEMP = 'S'
                    465:       ELSE
                    466:          CHTEMP = 'E'
                    467:       END IF
                    468:       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
                    469:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
                    470:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
                    471:       IF( IERR.NE.0 ) THEN
                    472:          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
                    473:             INFO = IERR
                    474:          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
                    475:             INFO = IERR - N
                    476:          ELSE
                    477:             INFO = N + 1
                    478:          END IF
                    479:          GO TO 110
                    480:       END IF
                    481: *
                    482: *     Compute Eigenvectors
                    483: *     (Workspace: need 6*N)
                    484: *
                    485:       IF( ILV ) THEN
                    486:          IF( ILVL ) THEN
                    487:             IF( ILVR ) THEN
                    488:                CHTEMP = 'B'
                    489:             ELSE
                    490:                CHTEMP = 'L'
                    491:             END IF
                    492:          ELSE
                    493:             CHTEMP = 'R'
                    494:          END IF
                    495:          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
                    496:      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
                    497:          IF( IERR.NE.0 ) THEN
                    498:             INFO = N + 2
                    499:             GO TO 110
                    500:          END IF
                    501: *
                    502: *        Undo balancing on VL and VR and normalization
                    503: *        (Workspace: none needed)
                    504: *
                    505:          IF( ILVL ) THEN
                    506:             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
                    507:      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
                    508:             DO 50 JC = 1, N
                    509:                IF( ALPHAI( JC ).LT.ZERO )
                    510:      $            GO TO 50
                    511:                TEMP = ZERO
                    512:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    513:                   DO 10 JR = 1, N
                    514:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
                    515:    10             CONTINUE
                    516:                ELSE
                    517:                   DO 20 JR = 1, N
                    518:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
                    519:      $                      ABS( VL( JR, JC+1 ) ) )
                    520:    20             CONTINUE
                    521:                END IF
                    522:                IF( TEMP.LT.SMLNUM )
                    523:      $            GO TO 50
                    524:                TEMP = ONE / TEMP
                    525:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    526:                   DO 30 JR = 1, N
                    527:                      VL( JR, JC ) = VL( JR, JC )*TEMP
                    528:    30             CONTINUE
                    529:                ELSE
                    530:                   DO 40 JR = 1, N
                    531:                      VL( JR, JC ) = VL( JR, JC )*TEMP
                    532:                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
                    533:    40             CONTINUE
                    534:                END IF
                    535:    50       CONTINUE
                    536:          END IF
                    537:          IF( ILVR ) THEN
                    538:             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
                    539:      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
                    540:             DO 100 JC = 1, N
                    541:                IF( ALPHAI( JC ).LT.ZERO )
                    542:      $            GO TO 100
                    543:                TEMP = ZERO
                    544:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    545:                   DO 60 JR = 1, N
                    546:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
                    547:    60             CONTINUE
                    548:                ELSE
                    549:                   DO 70 JR = 1, N
                    550:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
                    551:      $                      ABS( VR( JR, JC+1 ) ) )
                    552:    70             CONTINUE
                    553:                END IF
                    554:                IF( TEMP.LT.SMLNUM )
                    555:      $            GO TO 100
                    556:                TEMP = ONE / TEMP
                    557:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
                    558:                   DO 80 JR = 1, N
                    559:                      VR( JR, JC ) = VR( JR, JC )*TEMP
                    560:    80             CONTINUE
                    561:                ELSE
                    562:                   DO 90 JR = 1, N
                    563:                      VR( JR, JC ) = VR( JR, JC )*TEMP
                    564:                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
                    565:    90             CONTINUE
                    566:                END IF
                    567:   100       CONTINUE
                    568:          END IF
                    569: *
                    570: *        End of eigenvector calculation
                    571: *
                    572:       END IF
                    573: *
                    574: *     Undo scaling if necessary
                    575: *
1.10      bertrand  576:   110 CONTINUE
                    577: *
1.1       bertrand  578:       IF( ILASCL ) THEN
                    579:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
                    580:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
                    581:       END IF
                    582: *
                    583:       IF( ILBSCL ) THEN
                    584:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
                    585:       END IF
                    586: *
                    587:       WORK( 1 ) = MAXWRK
                    588:       RETURN
                    589: *
                    590: *     End of DGGEV
                    591: *
                    592:       END

CVSweb interface <joel.bertrand@systella.fr>