1*> \brief <b> CGGEV 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 CGGEV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggev.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggev.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggev.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
22*                         VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBVL, JOBVR
26*       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       REAL               RWORK( * )
30*       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
31*      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
32*      $                   WORK( * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
42*> (A,B), the generalized eigenvalues, and optionally, the left and/or
43*> right generalized eigenvectors.
44*>
45*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46*> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47*> singular. It is usually represented as the pair (alpha,beta), as
48*> there is a reasonable interpretation for beta=0, and even for both
49*> being zero.
50*>
51*> The right generalized eigenvector v(j) corresponding to the
52*> generalized eigenvalue lambda(j) of (A,B) satisfies
53*>
54*>              A * v(j) = lambda(j) * B * v(j).
55*>
56*> The left generalized eigenvector u(j) corresponding to the
57*> generalized eigenvalues lambda(j) of (A,B) satisfies
58*>
59*>              u(j)**H * A = lambda(j) * u(j)**H * B
60*>
61*> where u(j)**H is the conjugate-transpose of u(j).
62*> \endverbatim
63*
64*  Arguments:
65*  ==========
66*
67*> \param[in] JOBVL
68*> \verbatim
69*>          JOBVL is CHARACTER*1
70*>          = 'N':  do not compute the left generalized eigenvectors;
71*>          = 'V':  compute the left generalized eigenvectors.
72*> \endverbatim
73*>
74*> \param[in] JOBVR
75*> \verbatim
76*>          JOBVR is CHARACTER*1
77*>          = 'N':  do not compute the right generalized eigenvectors;
78*>          = 'V':  compute the right generalized eigenvectors.
79*> \endverbatim
80*>
81*> \param[in] N
82*> \verbatim
83*>          N is INTEGER
84*>          The order of the matrices A, B, VL, and VR.  N >= 0.
85*> \endverbatim
86*>
87*> \param[in,out] A
88*> \verbatim
89*>          A is COMPLEX array, dimension (LDA, N)
90*>          On entry, the matrix A in the pair (A,B).
91*>          On exit, A has been overwritten.
92*> \endverbatim
93*>
94*> \param[in] LDA
95*> \verbatim
96*>          LDA is INTEGER
97*>          The leading dimension of A.  LDA >= max(1,N).
98*> \endverbatim
99*>
100*> \param[in,out] B
101*> \verbatim
102*>          B is COMPLEX array, dimension (LDB, N)
103*>          On entry, the matrix B in the pair (A,B).
104*>          On exit, B has been overwritten.
105*> \endverbatim
106*>
107*> \param[in] LDB
108*> \verbatim
109*>          LDB is INTEGER
110*>          The leading dimension of B.  LDB >= max(1,N).
111*> \endverbatim
112*>
113*> \param[out] ALPHA
114*> \verbatim
115*>          ALPHA is COMPLEX array, dimension (N)
116*> \endverbatim
117*>
118*> \param[out] BETA
119*> \verbatim
120*>          BETA is COMPLEX array, dimension (N)
121*>          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
122*>          generalized eigenvalues.
123*>
124*>          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
125*>          underflow, and BETA(j) may even be zero.  Thus, the user
126*>          should avoid naively computing the ratio alpha/beta.
127*>          However, ALPHA will be always less than and usually
128*>          comparable with norm(A) in magnitude, and BETA always less
129*>          than and usually comparable with norm(B).
130*> \endverbatim
131*>
132*> \param[out] VL
133*> \verbatim
134*>          VL is COMPLEX array, dimension (LDVL,N)
135*>          If JOBVL = 'V', the left generalized eigenvectors u(j) are
136*>          stored one after another in the columns of VL, in the same
137*>          order as their eigenvalues.
138*>          Each eigenvector is scaled so the largest component has
139*>          abs(real part) + abs(imag. part) = 1.
140*>          Not referenced if JOBVL = 'N'.
141*> \endverbatim
142*>
143*> \param[in] LDVL
144*> \verbatim
145*>          LDVL is INTEGER
146*>          The leading dimension of the matrix VL. LDVL >= 1, and
147*>          if JOBVL = 'V', LDVL >= N.
148*> \endverbatim
149*>
150*> \param[out] VR
151*> \verbatim
152*>          VR is COMPLEX array, dimension (LDVR,N)
153*>          If JOBVR = 'V', the right generalized eigenvectors v(j) are
154*>          stored one after another in the columns of VR, in the same
155*>          order as their eigenvalues.
156*>          Each eigenvector is scaled so the largest component has
157*>          abs(real part) + abs(imag. part) = 1.
158*>          Not referenced if JOBVR = 'N'.
159*> \endverbatim
160*>
161*> \param[in] LDVR
162*> \verbatim
163*>          LDVR is INTEGER
164*>          The leading dimension of the matrix VR. LDVR >= 1, and
165*>          if JOBVR = 'V', LDVR >= N.
166*> \endverbatim
167*>
168*> \param[out] WORK
169*> \verbatim
170*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
171*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
172*> \endverbatim
173*>
174*> \param[in] LWORK
175*> \verbatim
176*>          LWORK is INTEGER
177*>          The dimension of the array WORK.  LWORK >= max(1,2*N).
178*>          For good performance, LWORK must generally be larger.
179*>
180*>          If LWORK = -1, then a workspace query is assumed; the routine
181*>          only calculates the optimal size of the WORK array, returns
182*>          this value as the first entry of the WORK array, and no error
183*>          message related to LWORK is issued by XERBLA.
184*> \endverbatim
185*>
186*> \param[out] RWORK
187*> \verbatim
188*>          RWORK is REAL array, dimension (8*N)
189*> \endverbatim
190*>
191*> \param[out] INFO
192*> \verbatim
193*>          INFO is INTEGER
194*>          = 0:  successful exit
195*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
196*>          =1,...,N:
197*>                The QZ iteration failed.  No eigenvectors have been
198*>                calculated, but ALPHA(j) and BETA(j) should be
199*>                correct for j=INFO+1,...,N.
200*>          > N:  =N+1: other then QZ iteration failed in CHGEQZ,
201*>                =N+2: error return from CTGEVC.
202*> \endverbatim
203*
204*  Authors:
205*  ========
206*
207*> \author Univ. of Tennessee
208*> \author Univ. of California Berkeley
209*> \author Univ. of Colorado Denver
210*> \author NAG Ltd.
211*
212*> \ingroup complexGEeigen
213*
214*  =====================================================================
215      SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
216     $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
217*
218*  -- LAPACK driver routine --
219*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
220*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221*
222*     .. Scalar Arguments ..
223      CHARACTER          JOBVL, JOBVR
224      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
225*     ..
226*     .. Array Arguments ..
227      REAL               RWORK( * )
228      COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
229     $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
230     $                   WORK( * )
231*     ..
232*
233*  =====================================================================
234*
235*     .. Parameters ..
236      REAL               ZERO, ONE
237      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
238      COMPLEX            CZERO, CONE
239      PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
240     $                   CONE = ( 1.0E0, 0.0E0 ) )
241*     ..
242*     .. Local Scalars ..
243      LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
244      CHARACTER          CHTEMP
245      INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
246     $                   IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
247     $                   LWKMIN, LWKOPT
248      REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
249     $                   SMLNUM, TEMP
250      COMPLEX            X
251*     ..
252*     .. Local Arrays ..
253      LOGICAL            LDUMMA( 1 )
254*     ..
255*     .. External Subroutines ..
256      EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
257     $                   CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, SLABAD,
258     $                   XERBLA
259*     ..
260*     .. External Functions ..
261      LOGICAL            LSAME
262      INTEGER            ILAENV
263      REAL               CLANGE, SLAMCH
264      EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
265*     ..
266*     .. Intrinsic Functions ..
267      INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
268*     ..
269*     .. Statement Functions ..
270      REAL               ABS1
271*     ..
272*     .. Statement Function definitions ..
273      ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
274*     ..
275*     .. Executable Statements ..
276*
277*     Decode the input arguments
278*
279      IF( LSAME( JOBVL, 'N' ) ) THEN
280         IJOBVL = 1
281         ILVL = .FALSE.
282      ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
283         IJOBVL = 2
284         ILVL = .TRUE.
285      ELSE
286         IJOBVL = -1
287         ILVL = .FALSE.
288      END IF
289*
290      IF( LSAME( JOBVR, 'N' ) ) THEN
291         IJOBVR = 1
292         ILVR = .FALSE.
293      ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
294         IJOBVR = 2
295         ILVR = .TRUE.
296      ELSE
297         IJOBVR = -1
298         ILVR = .FALSE.
299      END IF
300      ILV = ILVL .OR. ILVR
301*
302*     Test the input arguments
303*
304      INFO = 0
305      LQUERY = ( LWORK.EQ.-1 )
306      IF( IJOBVL.LE.0 ) THEN
307         INFO = -1
308      ELSE IF( IJOBVR.LE.0 ) THEN
309         INFO = -2
310      ELSE IF( N.LT.0 ) THEN
311         INFO = -3
312      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
313         INFO = -5
314      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
315         INFO = -7
316      ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
317         INFO = -11
318      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
319         INFO = -13
320      END IF
321*
322*     Compute workspace
323*      (Note: Comments in the code beginning "Workspace:" describe the
324*       minimal amount of workspace needed at that point in the code,
325*       as well as the preferred amount for good performance.
326*       NB refers to the optimal block size for the immediately
327*       following subroutine, as returned by ILAENV. The workspace is
328*       computed assuming ILO = 1 and IHI = N, the worst case.)
329*
330      IF( INFO.EQ.0 ) THEN
331         LWKMIN = MAX( 1, 2*N )
332         LWKOPT = MAX( 1, N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) )
333         LWKOPT = MAX( LWKOPT, N +
334     $                 N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) )
335         IF( ILVL ) THEN
336            LWKOPT = MAX( LWKOPT, N +
337     $                 N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) )
338         END IF
339         WORK( 1 ) = LWKOPT
340*
341         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
342     $      INFO = -15
343      END IF
344*
345      IF( INFO.NE.0 ) THEN
346         CALL XERBLA( 'CGGEV ', -INFO )
347         RETURN
348      ELSE IF( LQUERY ) THEN
349         RETURN
350      END IF
351*
352*     Quick return if possible
353*
354      IF( N.EQ.0 )
355     $   RETURN
356*
357*     Get machine constants
358*
359      EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
360      SMLNUM = SLAMCH( 'S' )
361      BIGNUM = ONE / SMLNUM
362      CALL SLABAD( SMLNUM, BIGNUM )
363      SMLNUM = SQRT( SMLNUM ) / EPS
364      BIGNUM = ONE / SMLNUM
365*
366*     Scale A if max element outside range [SMLNUM,BIGNUM]
367*
368      ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
369      ILASCL = .FALSE.
370      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
371         ANRMTO = SMLNUM
372         ILASCL = .TRUE.
373      ELSE IF( ANRM.GT.BIGNUM ) THEN
374         ANRMTO = BIGNUM
375         ILASCL = .TRUE.
376      END IF
377      IF( ILASCL )
378     $   CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
379*
380*     Scale B if max element outside range [SMLNUM,BIGNUM]
381*
382      BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
383      ILBSCL = .FALSE.
384      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
385         BNRMTO = SMLNUM
386         ILBSCL = .TRUE.
387      ELSE IF( BNRM.GT.BIGNUM ) THEN
388         BNRMTO = BIGNUM
389         ILBSCL = .TRUE.
390      END IF
391      IF( ILBSCL )
392     $   CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
393*
394*     Permute the matrices A, B to isolate eigenvalues if possible
395*     (Real Workspace: need 6*N)
396*
397      ILEFT = 1
398      IRIGHT = N + 1
399      IRWRK = IRIGHT + N
400      CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
401     $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
402*
403*     Reduce B to triangular form (QR decomposition of B)
404*     (Complex Workspace: need N, prefer N*NB)
405*
406      IROWS = IHI + 1 - ILO
407      IF( ILV ) THEN
408         ICOLS = N + 1 - ILO
409      ELSE
410         ICOLS = IROWS
411      END IF
412      ITAU = 1
413      IWRK = ITAU + IROWS
414      CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
415     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
416*
417*     Apply the orthogonal transformation to matrix A
418*     (Complex Workspace: need N, prefer N*NB)
419*
420      CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
421     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
422     $             LWORK+1-IWRK, IERR )
423*
424*     Initialize VL
425*     (Complex Workspace: need N, prefer N*NB)
426*
427      IF( ILVL ) THEN
428         CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
429         IF( IROWS.GT.1 ) THEN
430            CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
431     $                   VL( ILO+1, ILO ), LDVL )
432         END IF
433         CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
434     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
435      END IF
436*
437*     Initialize VR
438*
439      IF( ILVR )
440     $   CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
441*
442*     Reduce to generalized Hessenberg form
443*
444      IF( ILV ) THEN
445*
446*        Eigenvectors requested -- work on whole matrix.
447*
448         CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
449     $                LDVL, VR, LDVR, IERR )
450      ELSE
451         CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
452     $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
453      END IF
454*
455*     Perform QZ algorithm (Compute eigenvalues, and optionally, the
456*     Schur form and Schur vectors)
457*     (Complex Workspace: need N)
458*     (Real Workspace: need N)
459*
460      IWRK = ITAU
461      IF( ILV ) THEN
462         CHTEMP = 'S'
463      ELSE
464         CHTEMP = 'E'
465      END IF
466      CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
467     $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
468     $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
469      IF( IERR.NE.0 ) THEN
470         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
471            INFO = IERR
472         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
473            INFO = IERR - N
474         ELSE
475            INFO = N + 1
476         END IF
477         GO TO 70
478      END IF
479*
480*     Compute Eigenvectors
481*     (Real Workspace: need 2*N)
482*     (Complex Workspace: need 2*N)
483*
484      IF( ILV ) THEN
485         IF( ILVL ) THEN
486            IF( ILVR ) THEN
487               CHTEMP = 'B'
488            ELSE
489               CHTEMP = 'L'
490            END IF
491         ELSE
492            CHTEMP = 'R'
493         END IF
494*
495         CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
496     $                VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
497     $                IERR )
498         IF( IERR.NE.0 ) THEN
499            INFO = N + 2
500            GO TO 70
501         END IF
502*
503*        Undo balancing on VL and VR and normalization
504*        (Workspace: none needed)
505*
506         IF( ILVL ) THEN
507            CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
508     $                   RWORK( IRIGHT ), N, VL, LDVL, IERR )
509            DO 30 JC = 1, N
510               TEMP = ZERO
511               DO 10 JR = 1, N
512                  TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
513   10          CONTINUE
514               IF( TEMP.LT.SMLNUM )
515     $            GO TO 30
516               TEMP = ONE / TEMP
517               DO 20 JR = 1, N
518                  VL( JR, JC ) = VL( JR, JC )*TEMP
519   20          CONTINUE
520   30       CONTINUE
521         END IF
522         IF( ILVR ) THEN
523            CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
524     $                   RWORK( IRIGHT ), N, VR, LDVR, IERR )
525            DO 60 JC = 1, N
526               TEMP = ZERO
527               DO 40 JR = 1, N
528                  TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
529   40          CONTINUE
530               IF( TEMP.LT.SMLNUM )
531     $            GO TO 60
532               TEMP = ONE / TEMP
533               DO 50 JR = 1, N
534                  VR( JR, JC ) = VR( JR, JC )*TEMP
535   50          CONTINUE
536   60       CONTINUE
537         END IF
538      END IF
539*
540*     Undo scaling if necessary
541*
542   70 CONTINUE
543*
544      IF( ILASCL )
545     $   CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
546*
547      IF( ILBSCL )
548     $   CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
549*
550      WORK( 1 ) = LWKOPT
551      RETURN
552*
553*     End of CGGEV
554*
555      END
556