1 *> \brief <b> DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DGGES + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgges.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgges.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgges.f">
21 * SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
22 * SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
23 * LDVSR, WORK, LWORK, BWORK, INFO )
25 * .. Scalar Arguments ..
26 * CHARACTER JOBVSL, JOBVSR, SORT
27 * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
29 * .. Array Arguments ..
31 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
32 * $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
33 * $ VSR( LDVSR, * ), WORK( * )
35 * .. Function Arguments ..
46 *> DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
47 *> the generalized eigenvalues, the generalized real Schur form (S,T),
48 *> optionally, the left and/or right matrices of Schur vectors (VSL and
49 *> VSR). This gives the generalized Schur factorization
51 *> (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
53 *> Optionally, it also orders the eigenvalues so that a selected cluster
54 *> of eigenvalues appears in the leading diagonal blocks of the upper
55 *> quasi-triangular matrix S and the upper triangular matrix T.The
56 *> leading columns of VSL and VSR then form an orthonormal basis for the
57 *> corresponding left and right eigenspaces (deflating subspaces).
59 *> (If only the generalized eigenvalues are needed, use the driver
60 *> DGGEV instead, which is faster.)
62 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
63 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
64 *> usually represented as the pair (alpha,beta), as there is a
65 *> reasonable interpretation for beta=0 or both being zero.
67 *> A pair of matrices (S,T) is in generalized real Schur form if T is
68 *> upper triangular with non-negative diagonal and S is block upper
69 *> triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
70 *> to real generalized eigenvalues, while 2-by-2 blocks of S will be
71 *> "standardized" by making the corresponding elements of T have the
76 *> and the pair of corresponding 2-by-2 blocks in S and T will have a
77 *> complex conjugate pair of generalized eigenvalues.
86 *> JOBVSL is CHARACTER*1
87 *> = 'N': do not compute the left Schur vectors;
88 *> = 'V': compute the left Schur vectors.
93 *> JOBVSR is CHARACTER*1
94 *> = 'N': do not compute the right Schur vectors;
95 *> = 'V': compute the right Schur vectors.
100 *> SORT is CHARACTER*1
101 *> Specifies whether or not to order the eigenvalues on the
102 *> diagonal of the generalized Schur form.
103 *> = 'N': Eigenvalues are not ordered;
104 *> = 'S': Eigenvalues are ordered (see SELCTG);
109 *> SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
110 *> SELCTG must be declared EXTERNAL in the calling subroutine.
111 *> If SORT = 'N', SELCTG is not referenced.
112 *> If SORT = 'S', SELCTG is used to select eigenvalues to sort
113 *> to the top left of the Schur form.
114 *> An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
115 *> SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
116 *> one of a complex conjugate pair of eigenvalues is selected,
117 *> then both complex eigenvalues are selected.
119 *> Note that in the ill-conditioned case, a selected complex
120 *> eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
121 *> BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
128 *> The order of the matrices A, B, VSL, and VSR. N >= 0.
133 *> A is DOUBLE PRECISION array, dimension (LDA, N)
134 *> On entry, the first of the pair of matrices.
135 *> On exit, A has been overwritten by its generalized Schur
142 *> The leading dimension of A. LDA >= max(1,N).
147 *> B is DOUBLE PRECISION array, dimension (LDB, N)
148 *> On entry, the second of the pair of matrices.
149 *> On exit, B has been overwritten by its generalized Schur
156 *> The leading dimension of B. LDB >= max(1,N).
162 *> If SORT = 'N', SDIM = 0.
163 *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
164 *> for which SELCTG is true. (Complex conjugate pairs for which
165 *> SELCTG is true for either eigenvalue count as 2.)
168 *> \param[out] ALPHAR
170 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
173 *> \param[out] ALPHAI
175 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
180 *> BETA is DOUBLE PRECISION array, dimension (N)
181 *> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
182 *> be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i,
183 *> and BETA(j),j=1,...,N are the diagonals of the complex Schur
184 *> form (S,T) that would result if the 2-by-2 diagonal blocks of
185 *> the real Schur form of (A,B) were further reduced to
186 *> triangular form using 2-by-2 complex unitary transformations.
187 *> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
188 *> positive, then the j-th and (j+1)-st eigenvalues are a
189 *> complex conjugate pair, with ALPHAI(j+1) negative.
191 *> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
192 *> may easily over- or underflow, and BETA(j) may even be zero.
193 *> Thus, the user should avoid naively computing the ratio.
194 *> However, ALPHAR and ALPHAI will be always less than and
195 *> usually comparable with norm(A) in magnitude, and BETA always
196 *> less than and usually comparable with norm(B).
201 *> VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
202 *> If JOBVSL = 'V', VSL will contain the left Schur vectors.
203 *> Not referenced if JOBVSL = 'N'.
209 *> The leading dimension of the matrix VSL. LDVSL >=1, and
210 *> if JOBVSL = 'V', LDVSL >= N.
215 *> VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
216 *> If JOBVSR = 'V', VSR will contain the right Schur vectors.
217 *> Not referenced if JOBVSR = 'N'.
223 *> The leading dimension of the matrix VSR. LDVSR >= 1, and
224 *> if JOBVSR = 'V', LDVSR >= N.
229 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
230 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
236 *> The dimension of the array WORK.
237 *> If N = 0, LWORK >= 1, else LWORK >= 8*N+16.
238 *> For good performance , LWORK must generally be larger.
240 *> If LWORK = -1, then a workspace query is assumed; the routine
241 *> only calculates the optimal size of the WORK array, returns
242 *> this value as the first entry of the WORK array, and no error
243 *> message related to LWORK is issued by XERBLA.
248 *> BWORK is LOGICAL array, dimension (N)
249 *> Not referenced if SORT = 'N'.
255 *> = 0: successful exit
256 *> < 0: if INFO = -i, the i-th argument had an illegal value.
258 *> The QZ iteration failed. (A,B) are not in Schur
259 *> form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
260 *> be correct for j=INFO+1,...,N.
261 *> > N: =N+1: other than QZ iteration failed in DHGEQZ.
262 *> =N+2: after reordering, roundoff changed values of
263 *> some complex eigenvalues so that leading
264 *> eigenvalues in the Generalized Schur form no
265 *> longer satisfy SELCTG=.TRUE. This could also
266 *> be caused due to scaling.
267 *> =N+3: reordering failed in DTGSEN.
273 *> \author Univ. of Tennessee
274 *> \author Univ. of California Berkeley
275 *> \author Univ. of Colorado Denver
278 *> \date November 2011
280 *> \ingroup doubleGEeigen
282 * =====================================================================
283 SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
284 $ SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
285 $ LDVSR, WORK, LWORK, BWORK, INFO )
287 * -- LAPACK driver routine (version 3.4.0) --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
292 * .. Scalar Arguments ..
293 CHARACTER JOBVSL, JOBVSR, SORT
294 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
296 * .. Array Arguments ..
298 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
299 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
300 $ VSR( LDVSR, * ), WORK( * )
302 * .. Function Arguments ..
307 * =====================================================================
310 DOUBLE PRECISION ZERO, ONE
311 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
313 * .. Local Scalars ..
314 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
315 $ LQUERY, LST2SL, WANTST
316 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
317 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
319 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
320 $ PVSR, SAFMAX, SAFMIN, SMLNUM
324 DOUBLE PRECISION DIF( 2 )
326 * .. External Subroutines ..
327 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
328 $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
331 * .. External Functions ..
334 DOUBLE PRECISION DLAMCH, DLANGE
335 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
337 * .. Intrinsic Functions ..
338 INTRINSIC ABS, MAX, SQRT
340 * .. Executable Statements ..
342 * Decode the input arguments
344 IF( LSAME( JOBVSL, 'N' ) ) THEN
347 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
355 IF( LSAME( JOBVSR, 'N' ) ) THEN
358 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
366 WANTST = LSAME( SORT, 'S' )
368 * Test the input arguments
371 LQUERY = ( LWORK.EQ.-1 )
372 IF( IJOBVL.LE.0 ) THEN
374 ELSE IF( IJOBVR.LE.0 ) THEN
376 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
378 ELSE IF( N.LT.0 ) THEN
380 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
382 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
384 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
386 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
391 * (Note: Comments in the code beginning "Workspace:" describe the
392 * minimal amount of workspace needed at that point in the code,
393 * as well as the preferred amount for good performance.
394 * NB refers to the optimal block size for the immediately
395 * following subroutine, as returned by ILAENV.)
399 MINWRK = MAX( 8*N, 6*N + 16 )
400 MAXWRK = MINWRK - N +
401 $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
402 MAXWRK = MAX( MAXWRK, MINWRK - N +
403 $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
405 MAXWRK = MAX( MAXWRK, MINWRK - N +
406 $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
414 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
419 CALL XERBLA( 'DGGES ', -INFO )
421 ELSE IF( LQUERY ) THEN
425 * Quick return if possible
432 * Get machine constants
435 SAFMIN = DLAMCH( 'S' )
436 SAFMAX = ONE / SAFMIN
437 CALL DLABAD( SAFMIN, SAFMAX )
438 SMLNUM = SQRT( SAFMIN ) / EPS
439 BIGNUM = ONE / SMLNUM
441 * Scale A if max element outside range [SMLNUM,BIGNUM]
443 ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
445 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
448 ELSE IF( ANRM.GT.BIGNUM ) THEN
453 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
455 * Scale B if max element outside range [SMLNUM,BIGNUM]
457 BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
459 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
462 ELSE IF( BNRM.GT.BIGNUM ) THEN
467 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
469 * Permute the matrix to make it more nearly triangular
470 * (Workspace: need 6*N + 2*N space for storing balancing factors)
475 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
476 $ WORK( IRIGHT ), WORK( IWRK ), IERR )
478 * Reduce B to triangular form (QR decomposition of B)
479 * (Workspace: need N, prefer N*NB)
481 IROWS = IHI + 1 - ILO
485 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
486 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
488 * Apply the orthogonal transformation to matrix A
489 * (Workspace: need N, prefer N*NB)
491 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
492 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
493 $ LWORK+1-IWRK, IERR )
496 * (Workspace: need N, prefer N*NB)
499 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
500 IF( IROWS.GT.1 ) THEN
501 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
502 $ VSL( ILO+1, ILO ), LDVSL )
504 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
505 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
511 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
513 * Reduce to generalized Hessenberg form
514 * (Workspace: none needed)
516 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
517 $ LDVSL, VSR, LDVSR, IERR )
519 * Perform QZ algorithm, computing Schur vectors if desired
520 * (Workspace: need N)
523 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
524 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
525 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
527 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
529 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
537 * Sort eigenvalues ALPHA/BETA if desired
538 * (Workspace: need 4*N+16 )
543 * Undo scaling on eigenvalues before SELCTGing
546 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
548 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
552 $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
557 BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
560 CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
561 $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
562 $ PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
569 * Apply back-permutation to VSL and VSR
570 * (Workspace: none needed)
573 $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
574 $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
577 $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
578 $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
580 * Check if unscaling would cause over/underflow, if so, rescale
581 * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
582 * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
586 IF( ALPHAI( I ).NE.ZERO ) THEN
587 IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
588 $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
589 WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
590 BETA( I ) = BETA( I )*WORK( 1 )
591 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
592 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
593 ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
594 $ ( ANRMTO / ANRM ) .OR.
595 $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
597 WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
598 BETA( I ) = BETA( I )*WORK( 1 )
599 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
600 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
608 IF( ALPHAI( I ).NE.ZERO ) THEN
609 IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
610 $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
611 WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
612 BETA( I ) = BETA( I )*WORK( 1 )
613 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
614 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
623 CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
624 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
625 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
629 CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
630 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
635 * Check if reordering is correct
642 CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
643 IF( ALPHAI( I ).EQ.ZERO ) THEN
647 IF( CURSL .AND. .NOT.LASTSL )
652 * Last eigenvalue of conjugate pair
654 CURSL = CURSL .OR. LASTSL
659 IF( CURSL .AND. .NOT.LST2SL )
663 * First eigenvalue of conjugate pair