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

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: *
        !           221: *> \date November 2011
        !           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.8     ! bertrand  229: *  -- LAPACK driver routine (version 3.4.0) --
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.8     ! bertrand  232: *     November 2011
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: *
                    576:       IF( ILASCL ) THEN
                    577:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
                    578:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
                    579:       END IF
                    580: *
                    581:       IF( ILBSCL ) THEN
                    582:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
                    583:       END IF
                    584: *
                    585:   110 CONTINUE
                    586: *
                    587:       WORK( 1 ) = MAXWRK
                    588: *
                    589:       RETURN
                    590: *
                    591: *     End of DGGEV
                    592: *
                    593:       END

CVSweb interface <joel.bertrand@systella.fr>