1 *> \brief <b> DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DGEEV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeev.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeev.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeev.f">
21 * SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
22 * LDVR, WORK, LWORK, INFO )
24 * .. Scalar Arguments ..
25 * CHARACTER JOBVL, JOBVR
26 * INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
28 * .. Array Arguments ..
29 * DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
30 * $ WI( * ), WORK( * ), WR( * )
39 *> DGEEV computes for an N-by-N real nonsymmetric matrix A, the
40 *> eigenvalues and, optionally, the left and/or right eigenvectors.
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)**H * A = lambda(j) * u(j)**H
47 *> where u(j)**H denotes the conjugate-transpose of u(j).
49 *> The computed eigenvectors are normalized to have Euclidean norm
50 *> equal to 1 and largest component real.
58 *> JOBVL is CHARACTER*1
59 *> = 'N': left eigenvectors of A are not computed;
60 *> = 'V': left eigenvectors of A are computed.
65 *> JOBVR is CHARACTER*1
66 *> = 'N': right eigenvectors of A are not computed;
67 *> = 'V': right eigenvectors of A are computed.
73 *> The order of the matrix A. N >= 0.
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.
86 *> The leading dimension of the array A. LDA >= max(1,N).
91 *> WR is DOUBLE PRECISION array, dimension (N)
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
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).
121 *> The leading dimension of the array VL. LDVL >= 1; if
122 *> JOBVL = 'V', LDVL >= N.
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).
142 *> The leading dimension of the array VR. LDVR >= 1; if
143 *> JOBVR = 'V', LDVR >= N.
148 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
149 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
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.
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.
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
179 *> \author Univ. of Tennessee
180 *> \author Univ. of California Berkeley
181 *> \author Univ. of Colorado Denver
186 * @precisions fortran d -> s
188 *> \ingroup doubleGEeigen
190 * =====================================================================
191 SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
192 $ LDVR, WORK, LWORK, INFO )
195 * -- LAPACK driver routine (version 3.6.1) --
196 * -- LAPACK is a software package provided by Univ. of Tennessee, --
197 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200 * .. Scalar Arguments ..
201 CHARACTER JOBVL, JOBVR
202 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
204 * .. Array Arguments ..
205 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
206 $ WI( * ), WORK( * ), WR( * )
209 * =====================================================================
212 DOUBLE PRECISION ZERO, ONE
213 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
215 * .. Local Scalars ..
216 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
218 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
219 $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
220 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
225 DOUBLE PRECISION DUM( 1 )
227 * .. External Subroutines ..
228 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
229 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
232 * .. External Functions ..
234 INTEGER IDAMAX, ILAENV
235 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
236 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
239 * .. Intrinsic Functions ..
242 * .. Executable Statements ..
244 * Test the input arguments
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
252 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
254 ELSE IF( N.LT.0 ) THEN
256 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
258 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
260 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
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,
279 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
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,
286 HSWORK = INT( WORK(1) )
287 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
288 CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
289 $ VL, LDVL, VR, LDVR, N, NOUT,
291 LWORK_TREVC = INT( WORK(1) )
292 MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
293 MAXWRK = MAX( MAXWRK, 4*N )
294 ELSE IF( WANTVR ) THEN
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,
300 HSWORK = INT( WORK(1) )
301 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
302 CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
303 $ VL, LDVL, VR, LDVR, N, NOUT,
305 LWORK_TREVC = INT( WORK(1) )
306 MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
307 MAXWRK = MAX( MAXWRK, 4*N )
310 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
312 HSWORK = INT( WORK(1) )
313 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
315 MAXWRK = MAX( MAXWRK, MINWRK )
319 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
325 CALL XERBLA( 'DGEEV ', -INFO )
327 ELSE IF( LQUERY ) THEN
331 * Quick return if possible
336 * Get machine constants
339 SMLNUM = DLAMCH( 'S' )
340 BIGNUM = ONE / SMLNUM
341 CALL DLABAD( SMLNUM, BIGNUM )
342 SMLNUM = SQRT( SMLNUM ) / EPS
343 BIGNUM = ONE / SMLNUM
345 * Scale A if max element outside range [SMLNUM,BIGNUM]
347 ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
349 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
352 ELSE IF( ANRM.GT.BIGNUM ) THEN
357 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
360 * (Workspace: need N)
363 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
365 * Reduce to upper Hessenberg form
366 * (Workspace: need 3*N, prefer 2*N+N*NB)
370 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
371 $ LWORK-IWRK+1, IERR )
375 * Want left eigenvectors
376 * Copy Householder vectors to VL
379 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
381 * Generate orthogonal matrix in VL
382 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
384 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
385 $ LWORK-IWRK+1, IERR )
387 * Perform QR iteration, accumulating Schur vectors in VL
388 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
391 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
392 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
396 * Want left and right eigenvectors
397 * Copy Schur vectors to VR
400 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
403 ELSE IF( WANTVR ) THEN
405 * Want right eigenvectors
406 * Copy Householder vectors to VR
409 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
411 * Generate orthogonal matrix in VR
412 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
414 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
415 $ LWORK-IWRK+1, IERR )
417 * Perform QR iteration, accumulating Schur vectors in VR
418 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
421 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
422 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
426 * Compute eigenvalues only
427 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
430 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
431 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
434 * If INFO .NE. 0 from DHSEQR, then quit
439 IF( WANTVL .OR. WANTVR ) THEN
441 * Compute left and/or right eigenvectors
442 * (Workspace: need 4*N, prefer N + N + 2*N*NB)
444 CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
445 $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
450 * Undo balancing of left eigenvectors
451 * (Workspace: need N)
453 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
456 * Normalize left eigenvectors and make largest component real
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 )
468 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
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 )
480 * Undo balancing of right eigenvectors
481 * (Workspace: need N)
483 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
486 * Normalize right eigenvectors and make largest component real
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 )
498 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
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 )
508 * Undo scaling if necessary
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 )
517 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
519 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,