File:  [local] / rpl / lapack / lapack / dggev.f
Revision 1.9: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 22:19:29 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Cohérence

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

CVSweb interface <joel.bertrand@systella.fr>