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