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

CVSweb interface <joel.bertrand@systella.fr>