001:       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
002:      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, 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:       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
015:      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
016:      $                   VR( LDVR, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)
023: *  the generalized eigenvalues, and optionally, the left and/or right
024: *  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 eigenvector v(j) corresponding to the eigenvalue lambda(j)
033: *  of (A,B) satisfies
034: *
035: *                   A * v(j) = lambda(j) * B * v(j).
036: *
037: *  The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
038: *  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: *
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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: *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
074: *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
075: *  BETA    (output) DOUBLE PRECISION array, dimension (N)
076: *          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
077: *          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
078: *          the j-th eigenvalue is real; if positive, then the j-th and
079: *          (j+1)-st eigenvalues are a complex conjugate pair, with
080: *          ALPHAI(j+1) negative.
081: *
082: *          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
083: *          may easily over- or underflow, and BETA(j) may even be zero.
084: *          Thus, the user should avoid naively computing the ratio
085: *          alpha/beta.  However, ALPHAR and ALPHAI will be always less
086: *          than and usually comparable with norm(A) in magnitude, and
087: *          BETA always less than and usually comparable with norm(B).
088: *
089: *  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
090: *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
091: *          after another in the columns of VL, in the same order as
092: *          their eigenvalues. If the j-th eigenvalue is real, then
093: *          u(j) = VL(:,j), the j-th column of VL. If the j-th and
094: *          (j+1)-th eigenvalues form a complex conjugate pair, then
095: *          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
096: *          Each eigenvector is scaled so the largest component has
097: *          abs(real part)+abs(imag. part)=1.
098: *          Not referenced if JOBVL = 'N'.
099: *
100: *  LDVL    (input) INTEGER
101: *          The leading dimension of the matrix VL. LDVL >= 1, and
102: *          if JOBVL = 'V', LDVL >= N.
103: *
104: *  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
105: *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
106: *          after another in the columns of VR, in the same order as
107: *          their eigenvalues. If the j-th eigenvalue is real, then
108: *          v(j) = VR(:,j), the j-th column of VR. If the j-th and
109: *          (j+1)-th eigenvalues form a complex conjugate pair, then
110: *          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
111: *          Each eigenvector is scaled so the largest component has
112: *          abs(real part)+abs(imag. part)=1.
113: *          Not referenced if JOBVR = 'N'.
114: *
115: *  LDVR    (input) INTEGER
116: *          The leading dimension of the matrix VR. LDVR >= 1, and
117: *          if JOBVR = 'V', LDVR >= N.
118: *
119: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
120: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
121: *
122: *  LWORK   (input) INTEGER
123: *          The dimension of the array WORK.  LWORK >= max(1,8*N).
124: *          For good performance, LWORK must generally be larger.
125: *
126: *          If LWORK = -1, then a workspace query is assumed; the routine
127: *          only calculates the optimal size of the WORK array, returns
128: *          this value as the first entry of the WORK array, and no error
129: *          message related to LWORK is issued by XERBLA.
130: *
131: *  INFO    (output) INTEGER
132: *          = 0:  successful exit
133: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
134: *          = 1,...,N:
135: *                The QZ iteration failed.  No eigenvectors have been
136: *                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
137: *                should be correct for j=INFO+1,...,N.
138: *          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
139: *                =N+2: error return from DTGEVC.
140: *
141: *  =====================================================================
142: *
143: *     .. Parameters ..
144:       DOUBLE PRECISION   ZERO, ONE
145:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
146: *     ..
147: *     .. Local Scalars ..
148:       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
149:       CHARACTER          CHTEMP
150:       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
151:      $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
152:      $                   MINWRK
153:       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
154:      $                   SMLNUM, TEMP
155: *     ..
156: *     .. Local Arrays ..
157:       LOGICAL            LDUMMA( 1 )
158: *     ..
159: *     .. External Subroutines ..
160:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
161:      $                   DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
162:      $                   XERBLA
163: *     ..
164: *     .. External Functions ..
165:       LOGICAL            LSAME
166:       INTEGER            ILAENV
167:       DOUBLE PRECISION   DLAMCH, DLANGE
168:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
169: *     ..
170: *     .. Intrinsic Functions ..
171:       INTRINSIC          ABS, MAX, SQRT
172: *     ..
173: *     .. Executable Statements ..
174: *
175: *     Decode the input arguments
176: *
177:       IF( LSAME( JOBVL, 'N' ) ) THEN
178:          IJOBVL = 1
179:          ILVL = .FALSE.
180:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
181:          IJOBVL = 2
182:          ILVL = .TRUE.
183:       ELSE
184:          IJOBVL = -1
185:          ILVL = .FALSE.
186:       END IF
187: *
188:       IF( LSAME( JOBVR, 'N' ) ) THEN
189:          IJOBVR = 1
190:          ILVR = .FALSE.
191:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
192:          IJOBVR = 2
193:          ILVR = .TRUE.
194:       ELSE
195:          IJOBVR = -1
196:          ILVR = .FALSE.
197:       END IF
198:       ILV = ILVL .OR. ILVR
199: *
200: *     Test the input arguments
201: *
202:       INFO = 0
203:       LQUERY = ( LWORK.EQ.-1 )
204:       IF( IJOBVL.LE.0 ) THEN
205:          INFO = -1
206:       ELSE IF( IJOBVR.LE.0 ) THEN
207:          INFO = -2
208:       ELSE IF( N.LT.0 ) THEN
209:          INFO = -3
210:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
211:          INFO = -5
212:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
213:          INFO = -7
214:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
215:          INFO = -12
216:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
217:          INFO = -14
218:       END IF
219: *
220: *     Compute workspace
221: *      (Note: Comments in the code beginning "Workspace:" describe the
222: *       minimal amount of workspace needed at that point in the code,
223: *       as well as the preferred amount for good performance.
224: *       NB refers to the optimal block size for the immediately
225: *       following subroutine, as returned by ILAENV. The workspace is
226: *       computed assuming ILO = 1 and IHI = N, the worst case.)
227: *
228:       IF( INFO.EQ.0 ) THEN
229:          MINWRK = MAX( 1, 8*N )
230:          MAXWRK = MAX( 1, N*( 7 +
231:      $                 ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) )
232:          MAXWRK = MAX( MAXWRK, N*( 7 +
233:      $                 ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) )
234:          IF( ILVL ) THEN
235:             MAXWRK = MAX( MAXWRK, N*( 7 +
236:      $                 ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) )
237:          END IF
238:          WORK( 1 ) = MAXWRK
239: *
240:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
241:      $      INFO = -16
242:       END IF
243: *
244:       IF( INFO.NE.0 ) THEN
245:          CALL XERBLA( 'DGGEV ', -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 = DLAMCH( 'P' )
259:       SMLNUM = DLAMCH( 'S' )
260:       BIGNUM = ONE / SMLNUM
261:       CALL DLABAD( 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 = DLANGE( 'M', N, N, A, LDA, WORK )
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 DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
278: *
279: *     Scale B if max element outside range [SMLNUM,BIGNUM]
280: *
281:       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
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 DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
292: *
293: *     Permute the matrices A, B to isolate eigenvalues if possible
294: *     (Workspace: need 6*N)
295: *
296:       ILEFT = 1
297:       IRIGHT = N + 1
298:       IWRK = IRIGHT + N
299:       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
300:      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
301: *
302: *     Reduce B to triangular form (QR decomposition of B)
303: *     (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 = IWRK
312:       IWRK = ITAU + IROWS
313:       CALL DGEQRF( 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: *     (Workspace: need N, prefer N*NB)
318: *
319:       CALL DORMQR( 'L', 'T', 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: *     (Workspace: need N, prefer N*NB)
325: *
326:       IF( ILVL ) THEN
327:          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
328:          IF( IROWS.GT.1 ) THEN
329:             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
330:      $                   VL( ILO+1, ILO ), LDVL )
331:          END IF
332:          CALL DORGQR( 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 DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
340: *
341: *     Reduce to generalized Hessenberg form
342: *     (Workspace: none needed)
343: *
344:       IF( ILV ) THEN
345: *
346: *        Eigenvectors requested -- work on whole matrix.
347: *
348:          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
349:      $                LDVL, VR, LDVR, IERR )
350:       ELSE
351:          CALL DGGHRD( '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 forms and Schur vectors)
357: *     (Workspace: need N)
358: *
359:       IWRK = ITAU
360:       IF( ILV ) THEN
361:          CHTEMP = 'S'
362:       ELSE
363:          CHTEMP = 'E'
364:       END IF
365:       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
366:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
367:      $             WORK( IWRK ), LWORK+1-IWRK, 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 110
377:       END IF
378: *
379: *     Compute Eigenvectors
380: *     (Workspace: need 6*N)
381: *
382:       IF( ILV ) THEN
383:          IF( ILVL ) THEN
384:             IF( ILVR ) THEN
385:                CHTEMP = 'B'
386:             ELSE
387:                CHTEMP = 'L'
388:             END IF
389:          ELSE
390:             CHTEMP = 'R'
391:          END IF
392:          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
393:      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
394:          IF( IERR.NE.0 ) THEN
395:             INFO = N + 2
396:             GO TO 110
397:          END IF
398: *
399: *        Undo balancing on VL and VR and normalization
400: *        (Workspace: none needed)
401: *
402:          IF( ILVL ) THEN
403:             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
404:      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
405:             DO 50 JC = 1, N
406:                IF( ALPHAI( JC ).LT.ZERO )
407:      $            GO TO 50
408:                TEMP = ZERO
409:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
410:                   DO 10 JR = 1, N
411:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
412:    10             CONTINUE
413:                ELSE
414:                   DO 20 JR = 1, N
415:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
416:      $                      ABS( VL( JR, JC+1 ) ) )
417:    20             CONTINUE
418:                END IF
419:                IF( TEMP.LT.SMLNUM )
420:      $            GO TO 50
421:                TEMP = ONE / TEMP
422:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
423:                   DO 30 JR = 1, N
424:                      VL( JR, JC ) = VL( JR, JC )*TEMP
425:    30             CONTINUE
426:                ELSE
427:                   DO 40 JR = 1, N
428:                      VL( JR, JC ) = VL( JR, JC )*TEMP
429:                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
430:    40             CONTINUE
431:                END IF
432:    50       CONTINUE
433:          END IF
434:          IF( ILVR ) THEN
435:             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
436:      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
437:             DO 100 JC = 1, N
438:                IF( ALPHAI( JC ).LT.ZERO )
439:      $            GO TO 100
440:                TEMP = ZERO
441:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
442:                   DO 60 JR = 1, N
443:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
444:    60             CONTINUE
445:                ELSE
446:                   DO 70 JR = 1, N
447:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
448:      $                      ABS( VR( JR, JC+1 ) ) )
449:    70             CONTINUE
450:                END IF
451:                IF( TEMP.LT.SMLNUM )
452:      $            GO TO 100
453:                TEMP = ONE / TEMP
454:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
455:                   DO 80 JR = 1, N
456:                      VR( JR, JC ) = VR( JR, JC )*TEMP
457:    80             CONTINUE
458:                ELSE
459:                   DO 90 JR = 1, N
460:                      VR( JR, JC ) = VR( JR, JC )*TEMP
461:                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
462:    90             CONTINUE
463:                END IF
464:   100       CONTINUE
465:          END IF
466: *
467: *        End of eigenvector calculation
468: *
469:       END IF
470: *
471: *     Undo scaling if necessary
472: *
473:       IF( ILASCL ) THEN
474:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
475:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
476:       END IF
477: *
478:       IF( ILBSCL ) THEN
479:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
480:       END IF
481: *
482:   110 CONTINUE
483: *
484:       WORK( 1 ) = MAXWRK
485: *
486:       RETURN
487: *
488: *     End of DGGEV
489: *
490:       END
491: