LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
zgeev.f
Go to the documentation of this file.
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
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 * =====================================================================
177  SUBROUTINE zgeev( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
178  \$ work, lwork, rwork, info )
179 *
180 * -- LAPACK driver routine (version 3.4.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * November 2011
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,
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