001:       SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
002:      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOBVL, JOBVR
010:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               RWORK( * )
014:       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
015:      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
016:      $                   WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
023: *  (A,B), the generalized eigenvalues, and optionally, the left and/or
024: *  right generalized eigenvectors.
025: *
026: *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar
027: *  lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
028: *  singular. It is usually represented as the pair (alpha,beta), as
029: *  there is a reasonable interpretation for beta=0, and even for both
030: *  being zero.
031: *
032: *  The right generalized eigenvector v(j) corresponding to the
033: *  generalized eigenvalue lambda(j) of (A,B) satisfies
034: *
035: *               A * v(j) = lambda(j) * B * v(j).
036: *
037: *  The left generalized eigenvector u(j) corresponding to the
038: *  generalized eigenvalues lambda(j) of (A,B) satisfies
039: *
040: *               u(j)**H * A = lambda(j) * u(j)**H * B
041: *
042: *  where u(j)**H is the conjugate-transpose of u(j).
043: *
044: *  Arguments
045: *  =========
046: *
047: *  JOBVL   (input) CHARACTER*1
048: *          = 'N':  do not compute the left generalized eigenvectors;
049: *          = 'V':  compute the left generalized eigenvectors.
050: *
051: *  JOBVR   (input) CHARACTER*1
052: *          = 'N':  do not compute the right generalized eigenvectors;
053: *          = 'V':  compute the right generalized eigenvectors.
054: *
055: *  N       (input) INTEGER
056: *          The order of the matrices A, B, VL, and VR.  N >= 0.
057: *
058: *  A       (input/output) COMPLEX array, dimension (LDA, N)
059: *          On entry, the matrix A in the pair (A,B).
060: *          On exit, A has been overwritten.
061: *
062: *  LDA     (input) INTEGER
063: *          The leading dimension of A.  LDA >= max(1,N).
064: *
065: *  B       (input/output) COMPLEX array, dimension (LDB, N)
066: *          On entry, the matrix B in the pair (A,B).
067: *          On exit, B has been overwritten.
068: *
069: *  LDB     (input) INTEGER
070: *          The leading dimension of B.  LDB >= max(1,N).
071: *
072: *  ALPHA   (output) COMPLEX array, dimension (N)
073: *  BETA    (output) COMPLEX array, dimension (N)
074: *          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
075: *          generalized eigenvalues.
076: *
077: *          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
078: *          underflow, and BETA(j) may even be zero.  Thus, the user
079: *          should avoid naively computing the ratio alpha/beta.
080: *          However, ALPHA will be always less than and usually
081: *          comparable with norm(A) in magnitude, and BETA always less
082: *          than and usually comparable with norm(B).
083: *
084: *  VL      (output) COMPLEX array, dimension (LDVL,N)
085: *          If JOBVL = 'V', the left generalized eigenvectors u(j) are
086: *          stored one after another in the columns of VL, in the same
087: *          order as their eigenvalues.
088: *          Each eigenvector is scaled so the largest component has
089: *          abs(real part) + abs(imag. part) = 1.
090: *          Not referenced if JOBVL = 'N'.
091: *
092: *  LDVL    (input) INTEGER
093: *          The leading dimension of the matrix VL. LDVL >= 1, and
094: *          if JOBVL = 'V', LDVL >= N.
095: *
096: *  VR      (output) COMPLEX array, dimension (LDVR,N)
097: *          If JOBVR = 'V', the right generalized eigenvectors v(j) are
098: *          stored one after another in the columns of VR, in the same
099: *          order as their eigenvalues.
100: *          Each eigenvector is scaled so the largest component has
101: *          abs(real part) + abs(imag. part) = 1.
102: *          Not referenced if JOBVR = 'N'.
103: *
104: *  LDVR    (input) INTEGER
105: *          The leading dimension of the matrix VR. LDVR >= 1, and
106: *          if JOBVR = 'V', LDVR >= N.
107: *
108: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
109: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
110: *
111: *  LWORK   (input) INTEGER
112: *          The dimension of the array WORK.  LWORK >= max(1,2*N).
113: *          For good performance, LWORK must generally be larger.
114: *
115: *          If LWORK = -1, then a workspace query is assumed; the routine
116: *          only calculates the optimal size of the WORK array, returns
117: *          this value as the first entry of the WORK array, and no error
118: *          message related to LWORK is issued by XERBLA.
119: *
120: *  RWORK   (workspace/output) REAL array, dimension (8*N)
121: *
122: *  INFO    (output) INTEGER
123: *          = 0:  successful exit
124: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
125: *          =1,...,N:
126: *                The QZ iteration failed.  No eigenvectors have been
127: *                calculated, but ALPHA(j) and BETA(j) should be
128: *                correct for j=INFO+1,...,N.
129: *          > N:  =N+1: other then QZ iteration failed in SHGEQZ,
130: *                =N+2: error return from STGEVC.
131: *
132: *  =====================================================================
133: *
134: *     .. Parameters ..
135:       REAL               ZERO, ONE
136:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
137:       COMPLEX            CZERO, CONE
138:       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
139:      $                   CONE = ( 1.0E0, 0.0E0 ) )
140: *     ..
141: *     .. Local Scalars ..
142:       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
143:       CHARACTER          CHTEMP
144:       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
145:      $                   IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
146:      $                   LWKMIN, LWKOPT
147:       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
148:      $                   SMLNUM, TEMP
149:       COMPLEX            X
150: *     ..
151: *     .. Local Arrays ..
152:       LOGICAL            LDUMMA( 1 )
153: *     ..
154: *     .. External Subroutines ..
155:       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
156:      $                   CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, SLABAD,
157:      $                   XERBLA
158: *     ..
159: *     .. External Functions ..
160:       LOGICAL            LSAME
161:       INTEGER            ILAENV
162:       REAL               CLANGE, SLAMCH
163:       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
164: *     ..
165: *     .. Intrinsic Functions ..
166:       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
167: *     ..
168: *     .. Statement Functions ..
169:       REAL               ABS1
170: *     ..
171: *     .. Statement Function definitions ..
172:       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
173: *     ..
174: *     .. Executable Statements ..
175: *
176: *     Decode the input arguments
177: *
178:       IF( LSAME( JOBVL, 'N' ) ) THEN
179:          IJOBVL = 1
180:          ILVL = .FALSE.
181:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
182:          IJOBVL = 2
183:          ILVL = .TRUE.
184:       ELSE
185:          IJOBVL = -1
186:          ILVL = .FALSE.
187:       END IF
188: *
189:       IF( LSAME( JOBVR, 'N' ) ) THEN
190:          IJOBVR = 1
191:          ILVR = .FALSE.
192:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
193:          IJOBVR = 2
194:          ILVR = .TRUE.
195:       ELSE
196:          IJOBVR = -1
197:          ILVR = .FALSE.
198:       END IF
199:       ILV = ILVL .OR. ILVR
200: *
201: *     Test the input arguments
202: *
203:       INFO = 0
204:       LQUERY = ( LWORK.EQ.-1 )
205:       IF( IJOBVL.LE.0 ) THEN
206:          INFO = -1
207:       ELSE IF( IJOBVR.LE.0 ) THEN
208:          INFO = -2
209:       ELSE IF( N.LT.0 ) THEN
210:          INFO = -3
211:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
212:          INFO = -5
213:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
214:          INFO = -7
215:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
216:          INFO = -11
217:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
218:          INFO = -13
219:       END IF
220: *
221: *     Compute workspace
222: *      (Note: Comments in the code beginning "Workspace:" describe the
223: *       minimal amount of workspace needed at that point in the code,
224: *       as well as the preferred amount for good performance.
225: *       NB refers to the optimal block size for the immediately
226: *       following subroutine, as returned by ILAENV. The workspace is
227: *       computed assuming ILO = 1 and IHI = N, the worst case.)
228: *
229:       IF( INFO.EQ.0 ) THEN
230:          LWKMIN = MAX( 1, 2*N )
231:          LWKOPT = MAX( 1, N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) )
232:          LWKOPT = MAX( LWKOPT, N +
233:      $                 N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) ) 
234:          IF( ILVL ) THEN
235:             LWKOPT = MAX( LWKOPT, N +
236:      $                 N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) )
237:          END IF
238:          WORK( 1 ) = LWKOPT
239: *
240:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
241:      $      INFO = -15
242:       END IF
243: *
244:       IF( INFO.NE.0 ) THEN
245:          CALL XERBLA( 'CGGEV ', -INFO )
246:          RETURN
247:       ELSE IF( LQUERY ) THEN
248:          RETURN
249:       END IF
250: *
251: *     Quick return if possible
252: *
253:       IF( N.EQ.0 )
254:      $   RETURN
255: *
256: *     Get machine constants
257: *
258:       EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
259:       SMLNUM = SLAMCH( 'S' )
260:       BIGNUM = ONE / SMLNUM
261:       CALL SLABAD( SMLNUM, BIGNUM )
262:       SMLNUM = SQRT( SMLNUM ) / EPS
263:       BIGNUM = ONE / SMLNUM
264: *
265: *     Scale A if max element outside range [SMLNUM,BIGNUM]
266: *
267:       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
268:       ILASCL = .FALSE.
269:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
270:          ANRMTO = SMLNUM
271:          ILASCL = .TRUE.
272:       ELSE IF( ANRM.GT.BIGNUM ) THEN
273:          ANRMTO = BIGNUM
274:          ILASCL = .TRUE.
275:       END IF
276:       IF( ILASCL )
277:      $   CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
278: *
279: *     Scale B if max element outside range [SMLNUM,BIGNUM]
280: *
281:       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
282:       ILBSCL = .FALSE.
283:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
284:          BNRMTO = SMLNUM
285:          ILBSCL = .TRUE.
286:       ELSE IF( BNRM.GT.BIGNUM ) THEN
287:          BNRMTO = BIGNUM
288:          ILBSCL = .TRUE.
289:       END IF
290:       IF( ILBSCL )
291:      $   CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
292: *
293: *     Permute the matrices A, B to isolate eigenvalues if possible
294: *     (Real Workspace: need 6*N)
295: *
296:       ILEFT = 1
297:       IRIGHT = N + 1
298:       IRWRK = IRIGHT + N
299:       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
300:      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
301: *
302: *     Reduce B to triangular form (QR decomposition of B)
303: *     (Complex Workspace: need N, prefer N*NB)
304: *
305:       IROWS = IHI + 1 - ILO
306:       IF( ILV ) THEN
307:          ICOLS = N + 1 - ILO
308:       ELSE
309:          ICOLS = IROWS
310:       END IF
311:       ITAU = 1
312:       IWRK = ITAU + IROWS
313:       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
314:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
315: *
316: *     Apply the orthogonal transformation to matrix A
317: *     (Complex Workspace: need N, prefer N*NB)
318: *
319:       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
320:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
321:      $             LWORK+1-IWRK, IERR )
322: *
323: *     Initialize VL
324: *     (Complex Workspace: need N, prefer N*NB)
325: *
326:       IF( ILVL ) THEN
327:          CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
328:          IF( IROWS.GT.1 ) THEN
329:             CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
330:      $                   VL( ILO+1, ILO ), LDVL )
331:          END IF
332:          CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
333:      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
334:       END IF
335: *
336: *     Initialize VR
337: *
338:       IF( ILVR )
339:      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
340: *
341: *     Reduce to generalized Hessenberg form
342: *
343:       IF( ILV ) THEN
344: *
345: *        Eigenvectors requested -- work on whole matrix.
346: *
347:          CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
348:      $                LDVL, VR, LDVR, IERR )
349:       ELSE
350:          CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
351:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
352:       END IF
353: *
354: *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
355: *     Schur form and Schur vectors)
356: *     (Complex Workspace: need N)
357: *     (Real Workspace: need N)
358: *
359:       IWRK = ITAU
360:       IF( ILV ) THEN
361:          CHTEMP = 'S'
362:       ELSE
363:          CHTEMP = 'E'
364:       END IF
365:       CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
366:      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
367:      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
368:       IF( IERR.NE.0 ) THEN
369:          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
370:             INFO = IERR
371:          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
372:             INFO = IERR - N
373:          ELSE
374:             INFO = N + 1
375:          END IF
376:          GO TO 70
377:       END IF
378: *
379: *     Compute Eigenvectors
380: *     (Real Workspace: need 2*N)
381: *     (Complex Workspace: need 2*N)
382: *
383:       IF( ILV ) THEN
384:          IF( ILVL ) THEN
385:             IF( ILVR ) THEN
386:                CHTEMP = 'B'
387:             ELSE
388:                CHTEMP = 'L'
389:             END IF
390:          ELSE
391:             CHTEMP = 'R'
392:          END IF
393: *
394:          CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
395:      $                VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
396:      $                IERR )
397:          IF( IERR.NE.0 ) THEN
398:             INFO = N + 2
399:             GO TO 70
400:          END IF
401: *
402: *        Undo balancing on VL and VR and normalization
403: *        (Workspace: none needed)
404: *
405:          IF( ILVL ) THEN
406:             CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
407:      $                   RWORK( IRIGHT ), N, VL, LDVL, IERR )
408:             DO 30 JC = 1, N
409:                TEMP = ZERO
410:                DO 10 JR = 1, N
411:                   TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
412:    10          CONTINUE
413:                IF( TEMP.LT.SMLNUM )
414:      $            GO TO 30
415:                TEMP = ONE / TEMP
416:                DO 20 JR = 1, N
417:                   VL( JR, JC ) = VL( JR, JC )*TEMP
418:    20          CONTINUE
419:    30       CONTINUE
420:          END IF
421:          IF( ILVR ) THEN
422:             CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
423:      $                   RWORK( IRIGHT ), N, VR, LDVR, IERR )
424:             DO 60 JC = 1, N
425:                TEMP = ZERO
426:                DO 40 JR = 1, N
427:                   TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
428:    40          CONTINUE
429:                IF( TEMP.LT.SMLNUM )
430:      $            GO TO 60
431:                TEMP = ONE / TEMP
432:                DO 50 JR = 1, N
433:                   VR( JR, JC ) = VR( JR, JC )*TEMP
434:    50          CONTINUE
435:    60       CONTINUE
436:          END IF
437:       END IF
438: *
439: *     Undo scaling if necessary
440: *
441:       IF( ILASCL )
442:      $   CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
443: *
444:       IF( ILBSCL )
445:      $   CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
446: *
447:    70 CONTINUE
448:       WORK( 1 ) = LWKOPT
449: *
450:       RETURN
451: *
452: *     End of CGGEV
453: *
454:       END
455: