1*> \brief <b> ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZGGEV3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggev3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggev3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggev3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZGGEV3( 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*       DOUBLE PRECISION   RWORK( * )
30*       COMPLEX*16         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*> ZGGEV3 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*16 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*16 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*16 array, dimension (N)
116*> \endverbatim
117*>
118*> \param[out] BETA
119*> \verbatim
120*>          BETA is COMPLEX*16 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*16 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*16 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*16 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.
178*>
179*>          If LWORK = -1, then a workspace query is assumed; the routine
180*>          only calculates the optimal size of the WORK array, returns
181*>          this value as the first entry of the WORK array, and no error
182*>          message related to LWORK is issued by XERBLA.
183*> \endverbatim
184*>
185*> \param[out] RWORK
186*> \verbatim
187*>          RWORK is DOUBLE PRECISION array, dimension (8*N)
188*> \endverbatim
189*>
190*> \param[out] INFO
191*> \verbatim
192*>          INFO is INTEGER
193*>          = 0:  successful exit
194*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
195*>          =1,...,N:
196*>                The QZ iteration failed.  No eigenvectors have been
197*>                calculated, but ALPHA(j) and BETA(j) should be
198*>                correct for j=INFO+1,...,N.
199*>          > N:  =N+1: other then QZ iteration failed in ZHGEQZ,
200*>                =N+2: error return from ZTGEVC.
201*> \endverbatim
202*
203*  Authors:
204*  ========
205*
206*> \author Univ. of Tennessee
207*> \author Univ. of California Berkeley
208*> \author Univ. of Colorado Denver
209*> \author NAG Ltd.
210*
211*> \ingroup complex16GEeigen
212*
213*  =====================================================================
214      SUBROUTINE ZGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
215     $                   VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
216*
217*  -- LAPACK driver routine --
218*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
219*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221*     .. Scalar Arguments ..
222      CHARACTER          JOBVL, JOBVR
223      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
224*     ..
225*     .. Array Arguments ..
226      DOUBLE PRECISION   RWORK( * )
227      COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
228     $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229     $                   WORK( * )
230*     ..
231*
232*  =====================================================================
233*
234*     .. Parameters ..
235      DOUBLE PRECISION   ZERO, ONE
236      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
237      COMPLEX*16         CZERO, CONE
238      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
239     $                   CONE = ( 1.0D0, 0.0D0 ) )
240*     ..
241*     .. Local Scalars ..
242      LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
243      CHARACTER          CHTEMP
244      INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
245     $                   IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
246     $                   LWKOPT
247      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248     $                   SMLNUM, TEMP
249      COMPLEX*16         X
250*     ..
251*     .. Local Arrays ..
252      LOGICAL            LDUMMA( 1 )
253*     ..
254*     .. External Subroutines ..
255      EXTERNAL           DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHD3,
256     $                   ZLAQZ0, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
257     $                   ZUNMQR
258*     ..
259*     .. External Functions ..
260      LOGICAL            LSAME
261      DOUBLE PRECISION   DLAMCH, ZLANGE
262      EXTERNAL           LSAME, DLAMCH, ZLANGE
263*     ..
264*     .. Intrinsic Functions ..
265      INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
266*     ..
267*     .. Statement Functions ..
268      DOUBLE PRECISION   ABS1
269*     ..
270*     .. Statement Function definitions ..
271      ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
272*     ..
273*     .. Executable Statements ..
274*
275*     Decode the input arguments
276*
277      IF( LSAME( JOBVL, 'N' ) ) THEN
278         IJOBVL = 1
279         ILVL = .FALSE.
280      ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
281         IJOBVL = 2
282         ILVL = .TRUE.
283      ELSE
284         IJOBVL = -1
285         ILVL = .FALSE.
286      END IF
287*
288      IF( LSAME( JOBVR, 'N' ) ) THEN
289         IJOBVR = 1
290         ILVR = .FALSE.
291      ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
292         IJOBVR = 2
293         ILVR = .TRUE.
294      ELSE
295         IJOBVR = -1
296         ILVR = .FALSE.
297      END IF
298      ILV = ILVL .OR. ILVR
299*
300*     Test the input arguments
301*
302      INFO = 0
303      LQUERY = ( LWORK.EQ.-1 )
304      IF( IJOBVL.LE.0 ) THEN
305         INFO = -1
306      ELSE IF( IJOBVR.LE.0 ) THEN
307         INFO = -2
308      ELSE IF( N.LT.0 ) THEN
309         INFO = -3
310      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
311         INFO = -5
312      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
313         INFO = -7
314      ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
315         INFO = -11
316      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
317         INFO = -13
318      ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
319         INFO = -15
320      END IF
321*
322*     Compute workspace
323*
324      IF( INFO.EQ.0 ) THEN
325         CALL ZGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
326         LWKOPT = MAX( 1,  N+INT( WORK( 1 ) ) )
327         CALL ZUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
328     $                -1, IERR )
329         LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
330         IF( ILVL ) THEN
331            CALL ZUNGQR( N, N, N, VL, LDVL, WORK, WORK, -1, IERR )
332            LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
333         END IF
334         IF( ILV ) THEN
335            CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
336     $                   LDVL, VR, LDVR, WORK, -1, IERR )
337            LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
338            CALL ZLAQZ0( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
339     $                   ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
340     $                   RWORK, 0, IERR )
341            LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
342         ELSE
343            CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
344     $                   LDVL, VR, LDVR, WORK, -1, IERR )
345            LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
346            CALL ZLAQZ0( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
347     $                   ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
348     $                   RWORK, 0, IERR )
349            LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
350         END IF
351         WORK( 1 ) = DCMPLX( LWKOPT )
352      END IF
353*
354      IF( INFO.NE.0 ) THEN
355         CALL XERBLA( 'ZGGEV3 ', -INFO )
356         RETURN
357      ELSE IF( LQUERY ) THEN
358         RETURN
359      END IF
360*
361*     Quick return if possible
362*
363      IF( N.EQ.0 )
364     $   RETURN
365*
366*     Get machine constants
367*
368      EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
369      SMLNUM = DLAMCH( 'S' )
370      BIGNUM = ONE / SMLNUM
371      CALL DLABAD( SMLNUM, BIGNUM )
372      SMLNUM = SQRT( SMLNUM ) / EPS
373      BIGNUM = ONE / SMLNUM
374*
375*     Scale A if max element outside range [SMLNUM,BIGNUM]
376*
377      ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
378      ILASCL = .FALSE.
379      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
380         ANRMTO = SMLNUM
381         ILASCL = .TRUE.
382      ELSE IF( ANRM.GT.BIGNUM ) THEN
383         ANRMTO = BIGNUM
384         ILASCL = .TRUE.
385      END IF
386      IF( ILASCL )
387     $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
388*
389*     Scale B if max element outside range [SMLNUM,BIGNUM]
390*
391      BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
392      ILBSCL = .FALSE.
393      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
394         BNRMTO = SMLNUM
395         ILBSCL = .TRUE.
396      ELSE IF( BNRM.GT.BIGNUM ) THEN
397         BNRMTO = BIGNUM
398         ILBSCL = .TRUE.
399      END IF
400      IF( ILBSCL )
401     $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
402*
403*     Permute the matrices A, B to isolate eigenvalues if possible
404*
405      ILEFT = 1
406      IRIGHT = N + 1
407      IRWRK = IRIGHT + N
408      CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
409     $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
410*
411*     Reduce B to triangular form (QR decomposition of B)
412*
413      IROWS = IHI + 1 - ILO
414      IF( ILV ) THEN
415         ICOLS = N + 1 - ILO
416      ELSE
417         ICOLS = IROWS
418      END IF
419      ITAU = 1
420      IWRK = ITAU + IROWS
421      CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
422     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
423*
424*     Apply the orthogonal transformation to matrix A
425*
426      CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
427     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
428     $             LWORK+1-IWRK, IERR )
429*
430*     Initialize VL
431*
432      IF( ILVL ) THEN
433         CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
434         IF( IROWS.GT.1 ) THEN
435            CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
436     $                   VL( ILO+1, ILO ), LDVL )
437         END IF
438         CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
439     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
440      END IF
441*
442*     Initialize VR
443*
444      IF( ILVR )
445     $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
446*
447*     Reduce to generalized Hessenberg form
448*
449      IF( ILV ) THEN
450*
451*        Eigenvectors requested -- work on whole matrix.
452*
453         CALL ZGGHD3( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
454     $                LDVL, VR, LDVR, WORK( IWRK ), LWORK+1-IWRK, IERR )
455      ELSE
456         CALL ZGGHD3( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
457     $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR,
458     $                WORK( IWRK ), LWORK+1-IWRK, IERR )
459      END IF
460*
461*     Perform QZ algorithm (Compute eigenvalues, and optionally, the
462*     Schur form and Schur vectors)
463*
464      IWRK = ITAU
465      IF( ILV ) THEN
466         CHTEMP = 'S'
467      ELSE
468         CHTEMP = 'E'
469      END IF
470      CALL ZLAQZ0( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
471     $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
472     $             LWORK+1-IWRK, RWORK( IRWRK ), 0, IERR )
473      IF( IERR.NE.0 ) THEN
474         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
475            INFO = IERR
476         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
477            INFO = IERR - N
478         ELSE
479            INFO = N + 1
480         END IF
481         GO TO 70
482      END IF
483*
484*     Compute Eigenvectors
485*
486      IF( ILV ) THEN
487         IF( ILVL ) THEN
488            IF( ILVR ) THEN
489               CHTEMP = 'B'
490            ELSE
491               CHTEMP = 'L'
492            END IF
493         ELSE
494            CHTEMP = 'R'
495         END IF
496*
497         CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
498     $                VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
499     $                IERR )
500         IF( IERR.NE.0 ) THEN
501            INFO = N + 2
502            GO TO 70
503         END IF
504*
505*        Undo balancing on VL and VR and normalization
506*
507         IF( ILVL ) THEN
508            CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
509     $                   RWORK( IRIGHT ), N, VL, LDVL, IERR )
510            DO 30 JC = 1, N
511               TEMP = ZERO
512               DO 10 JR = 1, N
513                  TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
514   10          CONTINUE
515               IF( TEMP.LT.SMLNUM )
516     $            GO TO 30
517               TEMP = ONE / TEMP
518               DO 20 JR = 1, N
519                  VL( JR, JC ) = VL( JR, JC )*TEMP
520   20          CONTINUE
521   30       CONTINUE
522         END IF
523         IF( ILVR ) THEN
524            CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
525     $                   RWORK( IRIGHT ), N, VR, LDVR, IERR )
526            DO 60 JC = 1, N
527               TEMP = ZERO
528               DO 40 JR = 1, N
529                  TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
530   40          CONTINUE
531               IF( TEMP.LT.SMLNUM )
532     $            GO TO 60
533               TEMP = ONE / TEMP
534               DO 50 JR = 1, N
535                  VR( JR, JC ) = VR( JR, JC )*TEMP
536   50          CONTINUE
537   60       CONTINUE
538         END IF
539      END IF
540*
541*     Undo scaling if necessary
542*
543   70 CONTINUE
544*
545      IF( ILASCL )
546     $   CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
547*
548      IF( ILBSCL )
549     $   CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
550*
551      WORK( 1 ) = DCMPLX( LWKOPT )
552      RETURN
553*
554*     End of ZGGEV3
555*
556      END
557