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