File:  [local] / rpl / lapack / lapack / dggev3.f
Revision 1.5: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:51 2023 UTC (9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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 DLAQZ0.
  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: *
  227: *  -- LAPACK driver routine --
  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 ..
  259:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHD3, DLAQZ0, DLABAD,
  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 ) ) )
  336:             CALL DLAQZ0( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
  337:      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
  338:      $                   WORK, -1, 0, IERR )
  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 ) ) )
  344:             CALL DLAQZ0( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
  345:      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
  346:      $                   WORK, -1, 0, IERR )
  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
  469:       CALL DLAQZ0( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
  470:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
  471:      $             WORK( IWRK ), LWORK+1-IWRK, 0, IERR )
  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>