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

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

CVSweb interface <joel.bertrand@systella.fr>