LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dgeev.f
Go to the documentation of this file.
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)**H * A = lambda(j) * u(j)**H
47 *> where u(j)**H denotes the conjugate-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 *
185 * @precisions fortran d -> s
186 *
187 *> \ingroup doubleGEeigen
188 *
189 * =====================================================================
190  SUBROUTINE dgeev( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
191  $ LDVR, WORK, LWORK, INFO )
192  implicit none
193 *
194 * -- LAPACK driver routine --
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,
217  $ lwork_trevc, maxwrk, minwrk, nout
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,
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,
283  $ work, -1, info )
284  hswork = int( work(1) )
285  maxwrk = max( maxwrk, n + 1, n + hswork )
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 )
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,
297  $ work, -1, info )
298  hswork = int( work(1) )
299  maxwrk = max( maxwrk, n + 1, n + hswork )
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 )
305  maxwrk = max( maxwrk, 4*n )
306  ELSE
307  minwrk = 3*n
308  CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
309  $ work, -1, info )
310  hswork = int( work(1) )
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 *
432 * If INFO .NE. 0 from DHSEQR, then quit
433 *
434  IF( info.NE.0 )
435  $ GO TO 50
436 *
437  IF( wantvl .OR. wantvr ) THEN
438 *
439 * Compute left and/or right eigenvectors
440 * (Workspace: need 4*N, prefer N + N + 2*N*NB)
441 *
442  CALL dtrevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
443  $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
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
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:160
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:130
subroutine dgeev(JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: dgeev.f:192
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
Definition: dtrevc3.f:237