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*  =====================================================================
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,
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
490