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