LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cgeev.f
Go to the documentation of this file.
1 *> \brief <b> CGEEV 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 CGEEV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeev.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeev.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeev.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGEEV( 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 * REAL RWORK( * )
30 * COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
31 * $ W( * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CGEEV 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 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 array, dimension (N)
93 *> W contains the computed eigenvalues.
94 *> \endverbatim
95 *>
96 *> \param[out] VL
97 *> \verbatim
98 *> VL is COMPLEX 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 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 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 REAL 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 complexGEeigen
175 *
176 * =====================================================================
177  SUBROUTINE cgeev( 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  REAL rwork( * )
191  COMPLEX a( lda, * ), vl( ldvl, * ), vr( ldvr, * ),
192  $ w( * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  REAL zero, one
199  parameter( zero = 0.0e0, one = 1.0e0 )
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  REAL anrm, bignum, cscale, eps, scl, smlnum
207  COMPLEX tmp
208 * ..
209 * .. Local Arrays ..
210  LOGICAL select( 1 )
211  REAL dum( 1 )
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL cgebak, cgebal, cgehrd, chseqr, clacpy, clascl,
216 * ..
217 * .. External Functions ..
218  LOGICAL lsame
219  INTEGER ilaenv, isamax
220  REAL clange, scnrm2, slamch
221  EXTERNAL lsame, ilaenv, isamax, clange, scnrm2, slamch
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC aimag, cmplx, conjg, max, REAL, 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 *
249 * Compute workspace
250 * (Note: Comments in the code beginning "Workspace:" describe the
251 * minimal amount of workspace needed at that point in the code,
252 * as well as the preferred amount for good performance.
253 * CWorkspace refers to complex workspace, and RWorkspace to real
254 * workspace. NB refers to the optimal block size for the
255 * immediately following subroutine, as returned by ILAENV.
256 * HSWORK refers to the workspace preferred by CHSEQR, as
257 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
258 * the worst case.)
259 *
260  IF( info.EQ.0 ) THEN
261  IF( n.EQ.0 ) THEN
262  minwrk = 1
263  maxwrk = 1
264  ELSE
265  maxwrk = n + n*ilaenv( 1, 'CGEHRD', ' ', n, 1, n, 0 )
266  minwrk = 2*n
267  IF( wantvl ) THEN
268  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
269  $ ' ', n, 1, n, -1 ) )
270  CALL chseqr( 'S', 'V', n, 1, n, a, lda, w, vl, ldvl,
271  $ work, -1, info )
272  ELSE IF( wantvr ) THEN
273  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
274  $ ' ', n, 1, n, -1 ) )
275  CALL chseqr( 'S', 'V', n, 1, n, a, lda, w, vr, ldvr,
276  $ work, -1, info )
277  ELSE
278  CALL chseqr( 'E', 'N', n, 1, n, a, lda, w, vr, ldvr,
279  $ work, -1, info )
280  END IF
281  hswork = work( 1 )
282  maxwrk = max( maxwrk, hswork, minwrk )
283  END IF
284  work( 1 ) = maxwrk
285 *
286  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
287  info = -12
288  END IF
289  END IF
290 *
291  IF( info.NE.0 ) THEN
292  CALL xerbla( 'CGEEV ', -info )
293  RETURN
294  ELSE IF( lquery ) THEN
295  RETURN
296  END IF
297 *
298 * Quick return if possible
299 *
300  IF( n.EQ.0 )
301  $ RETURN
302 *
303 * Get machine constants
304 *
305  eps = slamch( 'P' )
306  smlnum = slamch( 'S' )
307  bignum = one / smlnum
308  CALL slabad( smlnum, bignum )
309  smlnum = sqrt( smlnum ) / eps
310  bignum = one / smlnum
311 *
312 * Scale A if max element outside range [SMLNUM,BIGNUM]
313 *
314  anrm = clange( 'M', n, n, a, lda, dum )
315  scalea = .false.
316  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
317  scalea = .true.
318  cscale = smlnum
319  ELSE IF( anrm.GT.bignum ) THEN
320  scalea = .true.
321  cscale = bignum
322  END IF
323  IF( scalea )
324  $ CALL clascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
325 *
326 * Balance the matrix
327 * (CWorkspace: none)
328 * (RWorkspace: need N)
329 *
330  ibal = 1
331  CALL cgebal( 'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
332 *
333 * Reduce to upper Hessenberg form
334 * (CWorkspace: need 2*N, prefer N+N*NB)
335 * (RWorkspace: none)
336 *
337  itau = 1
338  iwrk = itau + n
339  CALL cgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
340  $ lwork-iwrk+1, ierr )
341 *
342  IF( wantvl ) THEN
343 *
344 * Want left eigenvectors
345 * Copy Householder vectors to VL
346 *
347  side = 'L'
348  CALL clacpy( 'L', n, n, a, lda, vl, ldvl )
349 *
350 * Generate unitary matrix in VL
351 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
352 * (RWorkspace: none)
353 *
354  CALL cunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
355  $ lwork-iwrk+1, ierr )
356 *
357 * Perform QR iteration, accumulating Schur vectors in VL
358 * (CWorkspace: need 1, prefer HSWORK (see comments) )
359 * (RWorkspace: none)
360 *
361  iwrk = itau
362  CALL chseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,
363  $ work( iwrk ), lwork-iwrk+1, info )
364 *
365  IF( wantvr ) THEN
366 *
367 * Want left and right eigenvectors
368 * Copy Schur vectors to VR
369 *
370  side = 'B'
371  CALL clacpy( 'F', n, n, vl, ldvl, vr, ldvr )
372  END IF
373 *
374  ELSE IF( wantvr ) THEN
375 *
376 * Want right eigenvectors
377 * Copy Householder vectors to VR
378 *
379  side = 'R'
380  CALL clacpy( 'L', n, n, a, lda, vr, ldvr )
381 *
382 * Generate unitary matrix in VR
383 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
384 * (RWorkspace: none)
385 *
386  CALL cunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
387  $ lwork-iwrk+1, ierr )
388 *
389 * Perform QR iteration, accumulating Schur vectors in VR
390 * (CWorkspace: need 1, prefer HSWORK (see comments) )
391 * (RWorkspace: none)
392 *
393  iwrk = itau
394  CALL chseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,
395  $ work( iwrk ), lwork-iwrk+1, info )
396 *
397  ELSE
398 *
399 * Compute eigenvalues only
400 * (CWorkspace: need 1, prefer HSWORK (see comments) )
401 * (RWorkspace: none)
402 *
403  iwrk = itau
404  CALL chseqr( 'E', 'N', n, ilo, ihi, a, lda, w, vr, ldvr,
405  $ work( iwrk ), lwork-iwrk+1, info )
406  END IF
407 *
408 * If INFO > 0 from CHSEQR, then quit
409 *
410  IF( info.GT.0 )
411  $ go to 50
412 *
413  IF( wantvl .OR. wantvr ) THEN
414 *
415 * Compute left and/or right eigenvectors
416 * (CWorkspace: need 2*N)
417 * (RWorkspace: need 2*N)
418 *
419  irwork = ibal + n
420  CALL ctrevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
421  $ n, nout, work( iwrk ), rwork( irwork ), ierr )
422  END IF
423 *
424  IF( wantvl ) THEN
425 *
426 * Undo balancing of left eigenvectors
427 * (CWorkspace: none)
428 * (RWorkspace: need N)
429 *
430  CALL cgebak( 'B', 'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,
431  $ ierr )
432 *
433 * Normalize left eigenvectors and make largest component real
434 *
435  DO 20 i = 1, n
436  scl = one / scnrm2( n, vl( 1, i ), 1 )
437  CALL csscal( n, scl, vl( 1, i ), 1 )
438  DO 10 k = 1, n
439  rwork( irwork+k-1 ) = REAL( VL( K, I ) )**2 +
440  $ aimag( vl( k, i ) )**2
441  10 CONTINUE
442  k = isamax( n, rwork( irwork ), 1 )
443  tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
444  CALL cscal( n, tmp, vl( 1, i ), 1 )
445  vl( k, i ) = cmplx( REAL( VL( K, I ) ), zero )
446  20 CONTINUE
447  END IF
448 *
449  IF( wantvr ) THEN
450 *
451 * Undo balancing of right eigenvectors
452 * (CWorkspace: none)
453 * (RWorkspace: need N)
454 *
455  CALL cgebak( 'B', 'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,
456  $ ierr )
457 *
458 * Normalize right eigenvectors and make largest component real
459 *
460  DO 40 i = 1, n
461  scl = one / scnrm2( n, vr( 1, i ), 1 )
462  CALL csscal( n, scl, vr( 1, i ), 1 )
463  DO 30 k = 1, n
464  rwork( irwork+k-1 ) = REAL( VR( K, I ) )**2 +
465  $ aimag( vr( k, i ) )**2
466  30 CONTINUE
467  k = isamax( n, rwork( irwork ), 1 )
468  tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
469  CALL cscal( n, tmp, vr( 1, i ), 1 )
470  vr( k, i ) = cmplx( REAL( VR( K, I ) ), zero )
471  40 CONTINUE
472  END IF
473 *
474 * Undo scaling if necessary
475 *
476  50 CONTINUE
477  IF( scalea ) THEN
478  CALL clascl( 'G', 0, 0, cscale, anrm, n-info, 1, w( info+1 ),
479  $ max( n-info, 1 ), ierr )
480  IF( info.GT.0 ) THEN
481  CALL clascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, w, n, ierr )
482  END IF
483  END IF
484 *
485  work( 1 ) = maxwrk
486  RETURN
487 *
488 * End of CGEEV
489 *
490  END