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