Annotation of rpl/lapack/lapack/dgeev.f, revision 1.18

1.9       bertrand    1: *> \brief <b> DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.17      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.17      bertrand    9: *> Download DGEEV + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeev.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeev.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeev.f">
1.9       bertrand   15: *> [TXT]</a>
1.17      bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
                     22: *                         LDVR, WORK, LWORK, INFO )
1.17      bertrand   23: *
1.9       bertrand   24: *       .. Scalar Arguments ..
                     25: *       CHARACTER          JOBVL, JOBVR
                     26: *       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
                     27: *       ..
                     28: *       .. Array Arguments ..
                     29: *       DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
                     30: *      $                   WI( * ), WORK( * ), WR( * )
                     31: *       ..
1.17      bertrand   32: *
1.9       bertrand   33: *
                     34: *> \par Purpose:
                     35: *  =============
                     36: *>
                     37: *> \verbatim
                     38: *>
                     39: *> DGEEV computes for an N-by-N real nonsymmetric matrix A, the
                     40: *> eigenvalues and, optionally, the left and/or right eigenvectors.
                     41: *>
                     42: *> The right eigenvector v(j) of A satisfies
                     43: *>                  A * v(j) = lambda(j) * v(j)
                     44: *> where lambda(j) is its eigenvalue.
                     45: *> The left eigenvector u(j) of A satisfies
1.12      bertrand   46: *>               u(j)**H * A = lambda(j) * u(j)**H
                     47: *> where u(j)**H denotes the conjugate-transpose of u(j).
1.9       bertrand   48: *>
                     49: *> The computed eigenvectors are normalized to have Euclidean norm
                     50: *> equal to 1 and largest component real.
                     51: *> \endverbatim
                     52: *
                     53: *  Arguments:
                     54: *  ==========
                     55: *
                     56: *> \param[in] JOBVL
                     57: *> \verbatim
                     58: *>          JOBVL is CHARACTER*1
                     59: *>          = 'N': left eigenvectors of A are not computed;
                     60: *>          = 'V': left eigenvectors of A are computed.
                     61: *> \endverbatim
                     62: *>
                     63: *> \param[in] JOBVR
                     64: *> \verbatim
                     65: *>          JOBVR is CHARACTER*1
                     66: *>          = 'N': right eigenvectors of A are not computed;
                     67: *>          = 'V': right eigenvectors of A are computed.
                     68: *> \endverbatim
                     69: *>
                     70: *> \param[in] N
                     71: *> \verbatim
                     72: *>          N is INTEGER
                     73: *>          The order of the matrix A. N >= 0.
                     74: *> \endverbatim
                     75: *>
                     76: *> \param[in,out] A
                     77: *> \verbatim
                     78: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
                     79: *>          On entry, the N-by-N matrix A.
                     80: *>          On exit, A has been overwritten.
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[in] LDA
                     84: *> \verbatim
                     85: *>          LDA is INTEGER
                     86: *>          The leading dimension of the array A.  LDA >= max(1,N).
                     87: *> \endverbatim
                     88: *>
                     89: *> \param[out] WR
                     90: *> \verbatim
                     91: *>          WR is DOUBLE PRECISION array, dimension (N)
                     92: *> \endverbatim
                     93: *>
                     94: *> \param[out] WI
                     95: *> \verbatim
                     96: *>          WI is DOUBLE PRECISION array, dimension (N)
                     97: *>          WR and WI contain the real and imaginary parts,
                     98: *>          respectively, of the computed eigenvalues.  Complex
                     99: *>          conjugate pairs of eigenvalues appear consecutively
                    100: *>          with the eigenvalue having the positive imaginary part
                    101: *>          first.
                    102: *> \endverbatim
                    103: *>
                    104: *> \param[out] VL
                    105: *> \verbatim
                    106: *>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
                    107: *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
                    108: *>          after another in the columns of VL, in the same order
                    109: *>          as their eigenvalues.
                    110: *>          If JOBVL = 'N', VL is not referenced.
                    111: *>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
                    112: *>          the j-th column of VL.
                    113: *>          If the j-th and (j+1)-st eigenvalues form a complex
                    114: *>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
                    115: *>          u(j+1) = VL(:,j) - i*VL(:,j+1).
                    116: *> \endverbatim
                    117: *>
                    118: *> \param[in] LDVL
                    119: *> \verbatim
                    120: *>          LDVL is INTEGER
                    121: *>          The leading dimension of the array VL.  LDVL >= 1; if
                    122: *>          JOBVL = 'V', LDVL >= N.
                    123: *> \endverbatim
                    124: *>
                    125: *> \param[out] VR
                    126: *> \verbatim
                    127: *>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
                    128: *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
                    129: *>          after another in the columns of VR, in the same order
                    130: *>          as their eigenvalues.
                    131: *>          If JOBVR = 'N', VR is not referenced.
                    132: *>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
                    133: *>          the j-th column of VR.
                    134: *>          If the j-th and (j+1)-st eigenvalues form a complex
                    135: *>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
                    136: *>          v(j+1) = VR(:,j) - i*VR(:,j+1).
                    137: *> \endverbatim
                    138: *>
                    139: *> \param[in] LDVR
                    140: *> \verbatim
                    141: *>          LDVR is INTEGER
                    142: *>          The leading dimension of the array VR.  LDVR >= 1; if
                    143: *>          JOBVR = 'V', LDVR >= N.
                    144: *> \endverbatim
                    145: *>
                    146: *> \param[out] WORK
                    147: *> \verbatim
                    148: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                    149: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                    150: *> \endverbatim
                    151: *>
                    152: *> \param[in] LWORK
                    153: *> \verbatim
                    154: *>          LWORK is INTEGER
                    155: *>          The dimension of the array WORK.  LWORK >= max(1,3*N), and
                    156: *>          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
                    157: *>          performance, LWORK must generally be larger.
                    158: *>
                    159: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    160: *>          only calculates the optimal size of the WORK array, returns
                    161: *>          this value as the first entry of the WORK array, and no error
                    162: *>          message related to LWORK is issued by XERBLA.
                    163: *> \endverbatim
                    164: *>
                    165: *> \param[out] INFO
                    166: *> \verbatim
                    167: *>          INFO is INTEGER
                    168: *>          = 0:  successful exit
                    169: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    170: *>          > 0:  if INFO = i, the QR algorithm failed to compute all the
                    171: *>                eigenvalues, and no eigenvectors have been computed;
                    172: *>                elements i+1:N of WR and WI contain eigenvalues which
                    173: *>                have converged.
                    174: *> \endverbatim
                    175: *
                    176: *  Authors:
                    177: *  ========
                    178: *
1.17      bertrand  179: *> \author Univ. of Tennessee
                    180: *> \author Univ. of California Berkeley
                    181: *> \author Univ. of Colorado Denver
                    182: *> \author NAG Ltd.
1.9       bertrand  183: *
1.15      bertrand  184: *> \date June 2016
                    185: *
                    186: *  @precisions fortran d -> s
1.9       bertrand  187: *
                    188: *> \ingroup doubleGEeigen
                    189: *
                    190: *  =====================================================================
1.1       bertrand  191:       SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
                    192:      $                  LDVR, WORK, LWORK, INFO )
1.15      bertrand  193:       implicit none
1.1       bertrand  194: *
1.17      bertrand  195: *  -- LAPACK driver routine (version 3.7.0) --
1.1       bertrand  196: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    197: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15      bertrand  198: *     June 2016
1.1       bertrand  199: *
                    200: *     .. Scalar Arguments ..
                    201:       CHARACTER          JOBVL, JOBVR
                    202:       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
                    203: *     ..
                    204: *     .. Array Arguments ..
                    205:       DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
                    206:      $                   WI( * ), WORK( * ), WR( * )
                    207: *     ..
                    208: *
                    209: *  =====================================================================
                    210: *
                    211: *     .. Parameters ..
                    212:       DOUBLE PRECISION   ZERO, ONE
                    213:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
                    214: *     ..
                    215: *     .. Local Scalars ..
                    216:       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
                    217:       CHARACTER          SIDE
                    218:       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
1.15      bertrand  219:      $                   LWORK_TREVC, MAXWRK, MINWRK, NOUT
1.1       bertrand  220:       DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
                    221:      $                   SN
                    222: *     ..
                    223: *     .. Local Arrays ..
                    224:       LOGICAL            SELECT( 1 )
                    225:       DOUBLE PRECISION   DUM( 1 )
                    226: *     ..
                    227: *     .. External Subroutines ..
                    228:       EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
1.15      bertrand  229:      $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
1.1       bertrand  230:      $                   XERBLA
                    231: *     ..
                    232: *     .. External Functions ..
                    233:       LOGICAL            LSAME
                    234:       INTEGER            IDAMAX, ILAENV
                    235:       DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
                    236:       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
                    237:      $                   DNRM2
                    238: *     ..
                    239: *     .. Intrinsic Functions ..
                    240:       INTRINSIC          MAX, SQRT
                    241: *     ..
                    242: *     .. Executable Statements ..
                    243: *
                    244: *     Test the input arguments
                    245: *
                    246:       INFO = 0
                    247:       LQUERY = ( LWORK.EQ.-1 )
                    248:       WANTVL = LSAME( JOBVL, 'V' )
                    249:       WANTVR = LSAME( JOBVR, 'V' )
                    250:       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
                    251:          INFO = -1
                    252:       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
                    253:          INFO = -2
                    254:       ELSE IF( N.LT.0 ) THEN
                    255:          INFO = -3
                    256:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    257:          INFO = -5
                    258:       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
                    259:          INFO = -9
                    260:       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
                    261:          INFO = -11
                    262:       END IF
                    263: *
                    264: *     Compute workspace
                    265: *      (Note: Comments in the code beginning "Workspace:" describe the
                    266: *       minimal amount of workspace needed at that point in the code,
                    267: *       as well as the preferred amount for good performance.
                    268: *       NB refers to the optimal block size for the immediately
                    269: *       following subroutine, as returned by ILAENV.
                    270: *       HSWORK refers to the workspace preferred by DHSEQR, as
                    271: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
                    272: *       the worst case.)
                    273: *
                    274:       IF( INFO.EQ.0 ) THEN
                    275:          IF( N.EQ.0 ) THEN
                    276:             MINWRK = 1
                    277:             MAXWRK = 1
                    278:          ELSE
                    279:             MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
                    280:             IF( WANTVL ) THEN
                    281:                MINWRK = 4*N
                    282:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
                    283:      $                       'DORGHR', ' ', N, 1, N, -1 ) )
                    284:                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
1.15      bertrand  285:      $                      WORK, -1, INFO )
                    286:                HSWORK = INT( WORK(1) )
1.1       bertrand  287:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
1.15      bertrand  288:                CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
                    289:      $                       VL, LDVL, VR, LDVR, N, NOUT,
                    290:      $                       WORK, -1, IERR )
                    291:                LWORK_TREVC = INT( WORK(1) )
                    292:                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
1.1       bertrand  293:                MAXWRK = MAX( MAXWRK, 4*N )
                    294:             ELSE IF( WANTVR ) THEN
                    295:                MINWRK = 4*N
                    296:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
                    297:      $                       'DORGHR', ' ', N, 1, N, -1 ) )
                    298:                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
1.15      bertrand  299:      $                      WORK, -1, INFO )
                    300:                HSWORK = INT( WORK(1) )
1.1       bertrand  301:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
1.15      bertrand  302:                CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
                    303:      $                       VL, LDVL, VR, LDVR, N, NOUT,
                    304:      $                       WORK, -1, IERR )
                    305:                LWORK_TREVC = INT( WORK(1) )
                    306:                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
1.1       bertrand  307:                MAXWRK = MAX( MAXWRK, 4*N )
1.17      bertrand  308:             ELSE
1.1       bertrand  309:                MINWRK = 3*N
                    310:                CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
1.15      bertrand  311:      $                      WORK, -1, INFO )
                    312:                HSWORK = INT( WORK(1) )
1.1       bertrand  313:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
                    314:             END IF
                    315:             MAXWRK = MAX( MAXWRK, MINWRK )
                    316:          END IF
                    317:          WORK( 1 ) = MAXWRK
                    318: *
                    319:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
                    320:             INFO = -13
                    321:          END IF
                    322:       END IF
                    323: *
                    324:       IF( INFO.NE.0 ) THEN
                    325:          CALL XERBLA( 'DGEEV ', -INFO )
                    326:          RETURN
                    327:       ELSE IF( LQUERY ) THEN
                    328:          RETURN
                    329:       END IF
                    330: *
                    331: *     Quick return if possible
                    332: *
                    333:       IF( N.EQ.0 )
                    334:      $   RETURN
                    335: *
                    336: *     Get machine constants
                    337: *
                    338:       EPS = DLAMCH( 'P' )
                    339:       SMLNUM = DLAMCH( 'S' )
                    340:       BIGNUM = ONE / SMLNUM
                    341:       CALL DLABAD( SMLNUM, BIGNUM )
                    342:       SMLNUM = SQRT( SMLNUM ) / EPS
                    343:       BIGNUM = ONE / SMLNUM
                    344: *
                    345: *     Scale A if max element outside range [SMLNUM,BIGNUM]
                    346: *
                    347:       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
                    348:       SCALEA = .FALSE.
                    349:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
                    350:          SCALEA = .TRUE.
                    351:          CSCALE = SMLNUM
                    352:       ELSE IF( ANRM.GT.BIGNUM ) THEN
                    353:          SCALEA = .TRUE.
                    354:          CSCALE = BIGNUM
                    355:       END IF
                    356:       IF( SCALEA )
                    357:      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
                    358: *
                    359: *     Balance the matrix
                    360: *     (Workspace: need N)
                    361: *
                    362:       IBAL = 1
                    363:       CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
                    364: *
                    365: *     Reduce to upper Hessenberg form
                    366: *     (Workspace: need 3*N, prefer 2*N+N*NB)
                    367: *
                    368:       ITAU = IBAL + N
                    369:       IWRK = ITAU + N
                    370:       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
                    371:      $             LWORK-IWRK+1, IERR )
                    372: *
                    373:       IF( WANTVL ) THEN
                    374: *
                    375: *        Want left eigenvectors
                    376: *        Copy Householder vectors to VL
                    377: *
                    378:          SIDE = 'L'
                    379:          CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
                    380: *
                    381: *        Generate orthogonal matrix in VL
                    382: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
                    383: *
                    384:          CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
                    385:      $                LWORK-IWRK+1, IERR )
                    386: *
                    387: *        Perform QR iteration, accumulating Schur vectors in VL
                    388: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    389: *
                    390:          IWRK = ITAU
                    391:          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
                    392:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    393: *
                    394:          IF( WANTVR ) THEN
                    395: *
                    396: *           Want left and right eigenvectors
                    397: *           Copy Schur vectors to VR
                    398: *
                    399:             SIDE = 'B'
                    400:             CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
                    401:          END IF
                    402: *
                    403:       ELSE IF( WANTVR ) THEN
                    404: *
                    405: *        Want right eigenvectors
                    406: *        Copy Householder vectors to VR
                    407: *
                    408:          SIDE = 'R'
                    409:          CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
                    410: *
                    411: *        Generate orthogonal matrix in VR
                    412: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
                    413: *
                    414:          CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
                    415:      $                LWORK-IWRK+1, IERR )
                    416: *
                    417: *        Perform QR iteration, accumulating Schur vectors in VR
                    418: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    419: *
                    420:          IWRK = ITAU
                    421:          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
                    422:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    423: *
                    424:       ELSE
                    425: *
                    426: *        Compute eigenvalues only
                    427: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    428: *
                    429:          IWRK = ITAU
                    430:          CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
                    431:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    432:       END IF
                    433: *
1.15      bertrand  434: *     If INFO .NE. 0 from DHSEQR, then quit
1.1       bertrand  435: *
1.15      bertrand  436:       IF( INFO.NE.0 )
1.1       bertrand  437:      $   GO TO 50
                    438: *
                    439:       IF( WANTVL .OR. WANTVR ) THEN
                    440: *
                    441: *        Compute left and/or right eigenvectors
1.15      bertrand  442: *        (Workspace: need 4*N, prefer N + N + 2*N*NB)
1.1       bertrand  443: *
1.15      bertrand  444:          CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
                    445:      $                 N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
1.1       bertrand  446:       END IF
                    447: *
                    448:       IF( WANTVL ) THEN
                    449: *
                    450: *        Undo balancing of left eigenvectors
                    451: *        (Workspace: need N)
                    452: *
                    453:          CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
                    454:      $                IERR )
                    455: *
                    456: *        Normalize left eigenvectors and make largest component real
                    457: *
                    458:          DO 20 I = 1, N
                    459:             IF( WI( I ).EQ.ZERO ) THEN
                    460:                SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
                    461:                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
                    462:             ELSE IF( WI( I ).GT.ZERO ) THEN
                    463:                SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
                    464:      $               DNRM2( N, VL( 1, I+1 ), 1 ) )
                    465:                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
                    466:                CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
                    467:                DO 10 K = 1, N
                    468:                   WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
                    469:    10          CONTINUE
                    470:                K = IDAMAX( N, WORK( IWRK ), 1 )
                    471:                CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
                    472:                CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
                    473:                VL( K, I+1 ) = ZERO
                    474:             END IF
                    475:    20    CONTINUE
                    476:       END IF
                    477: *
                    478:       IF( WANTVR ) THEN
                    479: *
                    480: *        Undo balancing of right eigenvectors
                    481: *        (Workspace: need N)
                    482: *
                    483:          CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
                    484:      $                IERR )
                    485: *
                    486: *        Normalize right eigenvectors and make largest component real
                    487: *
                    488:          DO 40 I = 1, N
                    489:             IF( WI( I ).EQ.ZERO ) THEN
                    490:                SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
                    491:                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
                    492:             ELSE IF( WI( I ).GT.ZERO ) THEN
                    493:                SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
                    494:      $               DNRM2( N, VR( 1, I+1 ), 1 ) )
                    495:                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
                    496:                CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
                    497:                DO 30 K = 1, N
                    498:                   WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
                    499:    30          CONTINUE
                    500:                K = IDAMAX( N, WORK( IWRK ), 1 )
                    501:                CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
                    502:                CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
                    503:                VR( K, I+1 ) = ZERO
                    504:             END IF
                    505:    40    CONTINUE
                    506:       END IF
                    507: *
                    508: *     Undo scaling if necessary
                    509: *
                    510:    50 CONTINUE
                    511:       IF( SCALEA ) THEN
                    512:          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
                    513:      $                MAX( N-INFO, 1 ), IERR )
                    514:          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
                    515:      $                MAX( N-INFO, 1 ), IERR )
                    516:          IF( INFO.GT.0 ) THEN
                    517:             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
                    518:      $                   IERR )
                    519:             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
                    520:      $                   IERR )
                    521:          END IF
                    522:       END IF
                    523: *
                    524:       WORK( 1 ) = MAXWRK
                    525:       RETURN
                    526: *
                    527: *     End of DGEEV
                    528: *
                    529:       END

CVSweb interface <joel.bertrand@systella.fr>