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

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

CVSweb interface <joel.bertrand@systella.fr>