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

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: *
                    185: *  @precisions fortran d -> s
1.9       bertrand  186: *
                    187: *> \ingroup doubleGEeigen
                    188: *
                    189: *  =====================================================================
1.1       bertrand  190:       SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
                    191:      $                  LDVR, WORK, LWORK, INFO )
1.15      bertrand  192:       implicit none
1.1       bertrand  193: *
1.20    ! bertrand  194: *  -- LAPACK driver routine --
1.1       bertrand  195: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    196: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    197: *
                    198: *     .. Scalar Arguments ..
                    199:       CHARACTER          JOBVL, JOBVR
                    200:       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
                    201: *     ..
                    202: *     .. Array Arguments ..
                    203:       DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
                    204:      $                   WI( * ), WORK( * ), WR( * )
                    205: *     ..
                    206: *
                    207: *  =====================================================================
                    208: *
                    209: *     .. Parameters ..
                    210:       DOUBLE PRECISION   ZERO, ONE
                    211:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
                    212: *     ..
                    213: *     .. Local Scalars ..
                    214:       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
                    215:       CHARACTER          SIDE
                    216:       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
1.15      bertrand  217:      $                   LWORK_TREVC, MAXWRK, MINWRK, NOUT
1.1       bertrand  218:       DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
                    219:      $                   SN
                    220: *     ..
                    221: *     .. Local Arrays ..
                    222:       LOGICAL            SELECT( 1 )
                    223:       DOUBLE PRECISION   DUM( 1 )
                    224: *     ..
                    225: *     .. External Subroutines ..
                    226:       EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
1.15      bertrand  227:      $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
1.1       bertrand  228:      $                   XERBLA
                    229: *     ..
                    230: *     .. External Functions ..
                    231:       LOGICAL            LSAME
                    232:       INTEGER            IDAMAX, ILAENV
                    233:       DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
                    234:       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
                    235:      $                   DNRM2
                    236: *     ..
                    237: *     .. Intrinsic Functions ..
                    238:       INTRINSIC          MAX, SQRT
                    239: *     ..
                    240: *     .. Executable Statements ..
                    241: *
                    242: *     Test the input arguments
                    243: *
                    244:       INFO = 0
                    245:       LQUERY = ( LWORK.EQ.-1 )
                    246:       WANTVL = LSAME( JOBVL, 'V' )
                    247:       WANTVR = LSAME( JOBVR, 'V' )
                    248:       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
                    249:          INFO = -1
                    250:       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
                    251:          INFO = -2
                    252:       ELSE IF( N.LT.0 ) THEN
                    253:          INFO = -3
                    254:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    255:          INFO = -5
                    256:       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
                    257:          INFO = -9
                    258:       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
                    259:          INFO = -11
                    260:       END IF
                    261: *
                    262: *     Compute workspace
                    263: *      (Note: Comments in the code beginning "Workspace:" describe the
                    264: *       minimal amount of workspace needed at that point in the code,
                    265: *       as well as the preferred amount for good performance.
                    266: *       NB refers to the optimal block size for the immediately
                    267: *       following subroutine, as returned by ILAENV.
                    268: *       HSWORK refers to the workspace preferred by DHSEQR, as
                    269: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
                    270: *       the worst case.)
                    271: *
                    272:       IF( INFO.EQ.0 ) THEN
                    273:          IF( N.EQ.0 ) THEN
                    274:             MINWRK = 1
                    275:             MAXWRK = 1
                    276:          ELSE
                    277:             MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
                    278:             IF( WANTVL ) THEN
                    279:                MINWRK = 4*N
                    280:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
                    281:      $                       'DORGHR', ' ', N, 1, N, -1 ) )
                    282:                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
1.15      bertrand  283:      $                      WORK, -1, INFO )
                    284:                HSWORK = INT( WORK(1) )
1.1       bertrand  285:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
1.15      bertrand  286:                CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
                    287:      $                       VL, LDVL, VR, LDVR, N, NOUT,
                    288:      $                       WORK, -1, IERR )
                    289:                LWORK_TREVC = INT( WORK(1) )
                    290:                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
1.1       bertrand  291:                MAXWRK = MAX( MAXWRK, 4*N )
                    292:             ELSE IF( WANTVR ) THEN
                    293:                MINWRK = 4*N
                    294:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
                    295:      $                       'DORGHR', ' ', N, 1, N, -1 ) )
                    296:                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
1.15      bertrand  297:      $                      WORK, -1, INFO )
                    298:                HSWORK = INT( WORK(1) )
1.1       bertrand  299:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
1.15      bertrand  300:                CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
                    301:      $                       VL, LDVL, VR, LDVR, N, NOUT,
                    302:      $                       WORK, -1, IERR )
                    303:                LWORK_TREVC = INT( WORK(1) )
                    304:                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
1.1       bertrand  305:                MAXWRK = MAX( MAXWRK, 4*N )
1.17      bertrand  306:             ELSE
1.1       bertrand  307:                MINWRK = 3*N
                    308:                CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
1.15      bertrand  309:      $                      WORK, -1, INFO )
                    310:                HSWORK = INT( WORK(1) )
1.1       bertrand  311:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
                    312:             END IF
                    313:             MAXWRK = MAX( MAXWRK, MINWRK )
                    314:          END IF
                    315:          WORK( 1 ) = MAXWRK
                    316: *
                    317:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
                    318:             INFO = -13
                    319:          END IF
                    320:       END IF
                    321: *
                    322:       IF( INFO.NE.0 ) THEN
                    323:          CALL XERBLA( 'DGEEV ', -INFO )
                    324:          RETURN
                    325:       ELSE IF( LQUERY ) THEN
                    326:          RETURN
                    327:       END IF
                    328: *
                    329: *     Quick return if possible
                    330: *
                    331:       IF( N.EQ.0 )
                    332:      $   RETURN
                    333: *
                    334: *     Get machine constants
                    335: *
                    336:       EPS = DLAMCH( 'P' )
                    337:       SMLNUM = DLAMCH( 'S' )
                    338:       BIGNUM = ONE / SMLNUM
                    339:       CALL DLABAD( SMLNUM, BIGNUM )
                    340:       SMLNUM = SQRT( SMLNUM ) / EPS
                    341:       BIGNUM = ONE / SMLNUM
                    342: *
                    343: *     Scale A if max element outside range [SMLNUM,BIGNUM]
                    344: *
                    345:       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
                    346:       SCALEA = .FALSE.
                    347:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
                    348:          SCALEA = .TRUE.
                    349:          CSCALE = SMLNUM
                    350:       ELSE IF( ANRM.GT.BIGNUM ) THEN
                    351:          SCALEA = .TRUE.
                    352:          CSCALE = BIGNUM
                    353:       END IF
                    354:       IF( SCALEA )
                    355:      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
                    356: *
                    357: *     Balance the matrix
                    358: *     (Workspace: need N)
                    359: *
                    360:       IBAL = 1
                    361:       CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
                    362: *
                    363: *     Reduce to upper Hessenberg form
                    364: *     (Workspace: need 3*N, prefer 2*N+N*NB)
                    365: *
                    366:       ITAU = IBAL + N
                    367:       IWRK = ITAU + N
                    368:       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
                    369:      $             LWORK-IWRK+1, IERR )
                    370: *
                    371:       IF( WANTVL ) THEN
                    372: *
                    373: *        Want left eigenvectors
                    374: *        Copy Householder vectors to VL
                    375: *
                    376:          SIDE = 'L'
                    377:          CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
                    378: *
                    379: *        Generate orthogonal matrix in VL
                    380: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
                    381: *
                    382:          CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
                    383:      $                LWORK-IWRK+1, IERR )
                    384: *
                    385: *        Perform QR iteration, accumulating Schur vectors in VL
                    386: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    387: *
                    388:          IWRK = ITAU
                    389:          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
                    390:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    391: *
                    392:          IF( WANTVR ) THEN
                    393: *
                    394: *           Want left and right eigenvectors
                    395: *           Copy Schur vectors to VR
                    396: *
                    397:             SIDE = 'B'
                    398:             CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
                    399:          END IF
                    400: *
                    401:       ELSE IF( WANTVR ) THEN
                    402: *
                    403: *        Want right eigenvectors
                    404: *        Copy Householder vectors to VR
                    405: *
                    406:          SIDE = 'R'
                    407:          CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
                    408: *
                    409: *        Generate orthogonal matrix in VR
                    410: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
                    411: *
                    412:          CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
                    413:      $                LWORK-IWRK+1, IERR )
                    414: *
                    415: *        Perform QR iteration, accumulating Schur vectors in VR
                    416: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    417: *
                    418:          IWRK = ITAU
                    419:          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
                    420:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    421: *
                    422:       ELSE
                    423: *
                    424: *        Compute eigenvalues only
                    425: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
                    426: *
                    427:          IWRK = ITAU
                    428:          CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
                    429:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
                    430:       END IF
                    431: *
1.15      bertrand  432: *     If INFO .NE. 0 from DHSEQR, then quit
1.1       bertrand  433: *
1.15      bertrand  434:       IF( INFO.NE.0 )
1.1       bertrand  435:      $   GO TO 50
                    436: *
                    437:       IF( WANTVL .OR. WANTVR ) THEN
                    438: *
                    439: *        Compute left and/or right eigenvectors
1.15      bertrand  440: *        (Workspace: need 4*N, prefer N + N + 2*N*NB)
1.1       bertrand  441: *
1.15      bertrand  442:          CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
                    443:      $                 N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
1.1       bertrand  444:       END IF
                    445: *
                    446:       IF( WANTVL ) THEN
                    447: *
                    448: *        Undo balancing of left eigenvectors
                    449: *        (Workspace: need N)
                    450: *
                    451:          CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
                    452:      $                IERR )
                    453: *
                    454: *        Normalize left eigenvectors and make largest component real
                    455: *
                    456:          DO 20 I = 1, N
                    457:             IF( WI( I ).EQ.ZERO ) THEN
                    458:                SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
                    459:                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
                    460:             ELSE IF( WI( I ).GT.ZERO ) THEN
                    461:                SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
                    462:      $               DNRM2( N, VL( 1, I+1 ), 1 ) )
                    463:                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
                    464:                CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
                    465:                DO 10 K = 1, N
                    466:                   WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
                    467:    10          CONTINUE
                    468:                K = IDAMAX( N, WORK( IWRK ), 1 )
                    469:                CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
                    470:                CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
                    471:                VL( K, I+1 ) = ZERO
                    472:             END IF
                    473:    20    CONTINUE
                    474:       END IF
                    475: *
                    476:       IF( WANTVR ) THEN
                    477: *
                    478: *        Undo balancing of right eigenvectors
                    479: *        (Workspace: need N)
                    480: *
                    481:          CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
                    482:      $                IERR )
                    483: *
                    484: *        Normalize right eigenvectors and make largest component real
                    485: *
                    486:          DO 40 I = 1, N
                    487:             IF( WI( I ).EQ.ZERO ) THEN
                    488:                SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
                    489:                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
                    490:             ELSE IF( WI( I ).GT.ZERO ) THEN
                    491:                SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
                    492:      $               DNRM2( N, VR( 1, I+1 ), 1 ) )
                    493:                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
                    494:                CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
                    495:                DO 30 K = 1, N
                    496:                   WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
                    497:    30          CONTINUE
                    498:                K = IDAMAX( N, WORK( IWRK ), 1 )
                    499:                CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
                    500:                CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
                    501:                VR( K, I+1 ) = ZERO
                    502:             END IF
                    503:    40    CONTINUE
                    504:       END IF
                    505: *
                    506: *     Undo scaling if necessary
                    507: *
                    508:    50 CONTINUE
                    509:       IF( SCALEA ) THEN
                    510:          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
                    511:      $                MAX( N-INFO, 1 ), IERR )
                    512:          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
                    513:      $                MAX( N-INFO, 1 ), IERR )
                    514:          IF( INFO.GT.0 ) THEN
                    515:             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
                    516:      $                   IERR )
                    517:             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
                    518:      $                   IERR )
                    519:          END IF
                    520:       END IF
                    521: *
                    522:       WORK( 1 ) = MAXWRK
                    523:       RETURN
                    524: *
                    525: *     End of DGEEV
                    526: *
                    527:       END

CVSweb interface <joel.bertrand@systella.fr>