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*> \date September 2012
185*
186*> \ingroup doubleGEeigen
187*
188*  =====================================================================
189      SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
190     $                  LDVR, WORK, LWORK, INFO )
191*
192*  -- LAPACK driver routine (version 3.4.2) --
193*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
194*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*     September 2012
196*
197*     .. Scalar Arguments ..
198      CHARACTER          JOBVL, JOBVR
199      INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
200*     ..
201*     .. Array Arguments ..
202      DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
203     $                   WI( * ), WORK( * ), WR( * )
204*     ..
205*
206*  =====================================================================
207*
208*     .. Parameters ..
209      DOUBLE PRECISION   ZERO, ONE
210      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
211*     ..
212*     .. Local Scalars ..
213      LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
214      CHARACTER          SIDE
215      INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
216     $                   MAXWRK, MINWRK, NOUT
217      DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
218     $                   SN
219*     ..
220*     .. Local Arrays ..
221      LOGICAL            SELECT( 1 )
222      DOUBLE PRECISION   DUM( 1 )
223*     ..
224*     .. External Subroutines ..
225      EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
226     $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
227     $                   XERBLA
228*     ..
229*     .. External Functions ..
230      LOGICAL            LSAME
231      INTEGER            IDAMAX, ILAENV
232      DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
233      EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
234     $                   DNRM2
235*     ..
236*     .. Intrinsic Functions ..
237      INTRINSIC          MAX, SQRT
238*     ..
239*     .. Executable Statements ..
240*
241*     Test the input arguments
242*
243      INFO = 0
244      LQUERY = ( LWORK.EQ.-1 )
245      WANTVL = LSAME( JOBVL, 'V' )
246      WANTVR = LSAME( JOBVR, 'V' )
247      IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
248         INFO = -1
249      ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
250         INFO = -2
251      ELSE IF( N.LT.0 ) THEN
252         INFO = -3
253      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
254         INFO = -5
255      ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
256         INFO = -9
257      ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
258         INFO = -11
259      END IF
260*
261*     Compute workspace
262*      (Note: Comments in the code beginning "Workspace:" describe the
263*       minimal amount of workspace needed at that point in the code,
264*       as well as the preferred amount for good performance.
265*       NB refers to the optimal block size for the immediately
266*       following subroutine, as returned by ILAENV.
267*       HSWORK refers to the workspace preferred by DHSEQR, as
268*       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
269*       the worst case.)
270*
271      IF( INFO.EQ.0 ) THEN
272         IF( N.EQ.0 ) THEN
273            MINWRK = 1
274            MAXWRK = 1
275         ELSE
276            MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
277            IF( WANTVL ) THEN
278               MINWRK = 4*N
279               MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
280     $                       'DORGHR', ' ', N, 1, N, -1 ) )
281               CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
282     $                WORK, -1, INFO )
283               HSWORK = WORK( 1 )
284               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
285               MAXWRK = MAX( MAXWRK, 4*N )
286            ELSE IF( WANTVR ) THEN
287               MINWRK = 4*N
288               MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
289     $                       'DORGHR', ' ', N, 1, N, -1 ) )
290               CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
291     $                WORK, -1, INFO )
292               HSWORK = WORK( 1 )
293               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
294               MAXWRK = MAX( MAXWRK, 4*N )
295            ELSE
296               MINWRK = 3*N
297               CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
298     $                WORK, -1, INFO )
299               HSWORK = WORK( 1 )
300               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
301            END IF
302            MAXWRK = MAX( MAXWRK, MINWRK )
303         END IF
304         WORK( 1 ) = MAXWRK
305*
306         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
307            INFO = -13
308         END IF
309      END IF
310*
311      IF( INFO.NE.0 ) THEN
312         CALL XERBLA( 'DGEEV ', -INFO )
313         RETURN
314      ELSE IF( LQUERY ) THEN
315         RETURN
316      END IF
317*
318*     Quick return if possible
319*
320      IF( N.EQ.0 )
321     $   RETURN
322*
323*     Get machine constants
324*
325      EPS = DLAMCH( 'P' )
326      SMLNUM = DLAMCH( 'S' )
327      BIGNUM = ONE / SMLNUM
328      CALL DLABAD( SMLNUM, BIGNUM )
329      SMLNUM = SQRT( SMLNUM ) / EPS
330      BIGNUM = ONE / SMLNUM
331*
332*     Scale A if max element outside range [SMLNUM,BIGNUM]
333*
334      ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
335      SCALEA = .FALSE.
336      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
337         SCALEA = .TRUE.
338         CSCALE = SMLNUM
339      ELSE IF( ANRM.GT.BIGNUM ) THEN
340         SCALEA = .TRUE.
341         CSCALE = BIGNUM
342      END IF
343      IF( SCALEA )
344     $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
345*
346*     Balance the matrix
347*     (Workspace: need N)
348*
349      IBAL = 1
350      CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
351*
352*     Reduce to upper Hessenberg form
353*     (Workspace: need 3*N, prefer 2*N+N*NB)
354*
355      ITAU = IBAL + N
356      IWRK = ITAU + N
357      CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
358     $             LWORK-IWRK+1, IERR )
359*
360      IF( WANTVL ) THEN
361*
362*        Want left eigenvectors
363*        Copy Householder vectors to VL
364*
365         SIDE = 'L'
366         CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
367*
368*        Generate orthogonal matrix in VL
369*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
370*
371         CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
372     $                LWORK-IWRK+1, IERR )
373*
374*        Perform QR iteration, accumulating Schur vectors in VL
375*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
376*
377         IWRK = ITAU
378         CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
379     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
380*
381         IF( WANTVR ) THEN
382*
383*           Want left and right eigenvectors
384*           Copy Schur vectors to VR
385*
386            SIDE = 'B'
387            CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
388         END IF
389*
390      ELSE IF( WANTVR ) THEN
391*
392*        Want right eigenvectors
393*        Copy Householder vectors to VR
394*
395         SIDE = 'R'
396         CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
397*
398*        Generate orthogonal matrix in VR
399*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
400*
401         CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
402     $                LWORK-IWRK+1, IERR )
403*
404*        Perform QR iteration, accumulating Schur vectors in VR
405*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
406*
407         IWRK = ITAU
408         CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
409     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
410*
411      ELSE
412*
413*        Compute eigenvalues only
414*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
415*
416         IWRK = ITAU
417         CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
418     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
419      END IF
420*
421*     If INFO > 0 from DHSEQR, then quit
422*
423      IF( INFO.GT.0 )
424     $   GO TO 50
425*
426      IF( WANTVL .OR. WANTVR ) THEN
427*
428*        Compute left and/or right eigenvectors
429*        (Workspace: need 4*N)
430*
431         CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
432     $                N, NOUT, WORK( IWRK ), IERR )
433      END IF
434*
435      IF( WANTVL ) THEN
436*
437*        Undo balancing of left eigenvectors
438*        (Workspace: need N)
439*
440         CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
441     $                IERR )
442*
443*        Normalize left eigenvectors and make largest component real
444*
445         DO 20 I = 1, N
446            IF( WI( I ).EQ.ZERO ) THEN
447               SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
448               CALL DSCAL( N, SCL, VL( 1, I ), 1 )
449            ELSE IF( WI( I ).GT.ZERO ) THEN
450               SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
451     $               DNRM2( N, VL( 1, I+1 ), 1 ) )
452               CALL DSCAL( N, SCL, VL( 1, I ), 1 )
453               CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
454               DO 10 K = 1, N
455                  WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
456   10          CONTINUE
457               K = IDAMAX( N, WORK( IWRK ), 1 )
458               CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
459               CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
460               VL( K, I+1 ) = ZERO
461            END IF
462   20    CONTINUE
463      END IF
464*
465      IF( WANTVR ) THEN
466*
467*        Undo balancing of right eigenvectors
468*        (Workspace: need N)
469*
470         CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
471     $                IERR )
472*
473*        Normalize right eigenvectors and make largest component real
474*
475         DO 40 I = 1, N
476            IF( WI( I ).EQ.ZERO ) THEN
477               SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
478               CALL DSCAL( N, SCL, VR( 1, I ), 1 )
479            ELSE IF( WI( I ).GT.ZERO ) THEN
480               SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
481     $               DNRM2( N, VR( 1, I+1 ), 1 ) )
482               CALL DSCAL( N, SCL, VR( 1, I ), 1 )
483               CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
484               DO 30 K = 1, N
485                  WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
486   30          CONTINUE
487               K = IDAMAX( N, WORK( IWRK ), 1 )
488               CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
489               CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
490               VR( K, I+1 ) = ZERO
491            END IF
492   40    CONTINUE
493      END IF
494*
495*     Undo scaling if necessary
496*
497   50 CONTINUE
498      IF( SCALEA ) THEN
499         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
500     $                MAX( N-INFO, 1 ), IERR )
501         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
502     $                MAX( N-INFO, 1 ), IERR )
503         IF( INFO.GT.0 ) THEN
504            CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
505     $                   IERR )
506            CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
507     $                   IERR )
508         END IF
509      END IF
510*
511      WORK( 1 ) = MAXWRK
512      RETURN
513*
514*     End of DGEEV
515*
516      END
517c $Id$
518