1 *> \brief <b> CGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CGGEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggevx.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggevx.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggevx.f">
21 * SUBROUTINE CGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
22 * ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
23 * LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
24 * WORK, LWORK, RWORK, IWORK, BWORK, INFO )
26 * .. Scalar Arguments ..
27 * CHARACTER BALANC, JOBVL, JOBVR, SENSE
28 * INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
31 * .. Array Arguments ..
34 * REAL LSCALE( * ), RCONDE( * ), RCONDV( * ),
35 * $ RSCALE( * ), RWORK( * )
36 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
37 * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
47 *> CGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
48 *> (A,B) the generalized eigenvalues, and optionally, the left and/or
49 *> right generalized eigenvectors.
51 *> Optionally, it also computes a balancing transformation to improve
52 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
53 *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
54 *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
55 *> right eigenvectors (RCONDV).
57 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
58 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
59 *> singular. It is usually represented as the pair (alpha,beta), as
60 *> there is a reasonable interpretation for beta=0, and even for both
63 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
65 *> A * v(j) = lambda(j) * B * v(j) .
66 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
68 *> u(j)**H * A = lambda(j) * u(j)**H * B.
69 *> where u(j)**H is the conjugate-transpose of u(j).
78 *> BALANC is CHARACTER*1
79 *> Specifies the balance option to be performed:
80 *> = 'N': do not diagonally scale or permute;
81 *> = 'P': permute only;
83 *> = 'B': both permute and scale.
84 *> Computed reciprocal condition numbers will be for the
85 *> matrices after permuting and/or balancing. Permuting does
86 *> not change condition numbers (in exact arithmetic), but
92 *> JOBVL is CHARACTER*1
93 *> = 'N': do not compute the left generalized eigenvectors;
94 *> = 'V': compute the left generalized eigenvectors.
99 *> JOBVR is CHARACTER*1
100 *> = 'N': do not compute the right generalized eigenvectors;
101 *> = 'V': compute the right generalized eigenvectors.
106 *> SENSE is CHARACTER*1
107 *> Determines which reciprocal condition numbers are computed.
108 *> = 'N': none are computed;
109 *> = 'E': computed for eigenvalues only;
110 *> = 'V': computed for eigenvectors only;
111 *> = 'B': computed for eigenvalues and eigenvectors.
117 *> The order of the matrices A, B, VL, and VR. N >= 0.
122 *> A is COMPLEX array, dimension (LDA, N)
123 *> On entry, the matrix A in the pair (A,B).
124 *> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
125 *> or both, then A contains the first part of the complex Schur
126 *> form of the "balanced" versions of the input A and B.
132 *> The leading dimension of A. LDA >= max(1,N).
137 *> B is COMPLEX array, dimension (LDB, N)
138 *> On entry, the matrix B in the pair (A,B).
139 *> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
140 *> or both, then B contains the second part of the complex
141 *> Schur form of the "balanced" versions of the input A and B.
147 *> The leading dimension of B. LDB >= max(1,N).
152 *> ALPHA is COMPLEX array, dimension (N)
157 *> BETA is COMPLEX array, dimension (N)
158 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
161 *> Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
162 *> underflow, and BETA(j) may even be zero. Thus, the user
163 *> should avoid naively computing the ratio ALPHA/BETA.
164 *> However, ALPHA will be always less than and usually
165 *> comparable with norm(A) in magnitude, and BETA always less
166 *> than and usually comparable with norm(B).
171 *> VL is COMPLEX array, dimension (LDVL,N)
172 *> If JOBVL = 'V', the left generalized eigenvectors u(j) are
173 *> stored one after another in the columns of VL, in the same
174 *> order as their eigenvalues.
175 *> Each eigenvector will be scaled so the largest component
176 *> will have abs(real part) + abs(imag. part) = 1.
177 *> Not referenced if JOBVL = 'N'.
183 *> The leading dimension of the matrix VL. LDVL >= 1, and
184 *> if JOBVL = 'V', LDVL >= N.
189 *> VR is COMPLEX array, dimension (LDVR,N)
190 *> If JOBVR = 'V', the right generalized eigenvectors v(j) are
191 *> stored one after another in the columns of VR, in the same
192 *> order as their eigenvalues.
193 *> Each eigenvector will be scaled so the largest component
194 *> will have abs(real part) + abs(imag. part) = 1.
195 *> Not referenced if JOBVR = 'N'.
201 *> The leading dimension of the matrix VR. LDVR >= 1, and
202 *> if JOBVR = 'V', LDVR >= N.
213 *> ILO and IHI are integer values such that on exit
214 *> A(i,j) = 0 and B(i,j) = 0 if i > j and
215 *> j = 1,...,ILO-1 or i = IHI+1,...,N.
216 *> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
219 *> \param[out] LSCALE
221 *> LSCALE is REAL array, dimension (N)
222 *> Details of the permutations and scaling factors applied
223 *> to the left side of A and B. If PL(j) is the index of the
224 *> row interchanged with row j, and DL(j) is the scaling
225 *> factor applied to row j, then
226 *> LSCALE(j) = PL(j) for j = 1,...,ILO-1
227 *> = DL(j) for j = ILO,...,IHI
228 *> = PL(j) for j = IHI+1,...,N.
229 *> The order in which the interchanges are made is N to IHI+1,
233 *> \param[out] RSCALE
235 *> RSCALE is REAL array, dimension (N)
236 *> Details of the permutations and scaling factors applied
237 *> to the right side of A and B. If PR(j) is the index of the
238 *> column interchanged with column j, and DR(j) is the scaling
239 *> factor applied to column j, then
240 *> RSCALE(j) = PR(j) for j = 1,...,ILO-1
241 *> = DR(j) for j = ILO,...,IHI
242 *> = PR(j) for j = IHI+1,...,N
243 *> The order in which the interchanges are made is N to IHI+1,
250 *> The one-norm of the balanced matrix A.
256 *> The one-norm of the balanced matrix B.
259 *> \param[out] RCONDE
261 *> RCONDE is REAL array, dimension (N)
262 *> If SENSE = 'E' or 'B', the reciprocal condition numbers of
263 *> the eigenvalues, stored in consecutive elements of the array.
264 *> If SENSE = 'N' or 'V', RCONDE is not referenced.
267 *> \param[out] RCONDV
269 *> RCONDV is REAL array, dimension (N)
270 *> If SENSE = 'V' or 'B', the estimated reciprocal condition
271 *> numbers of the eigenvectors, stored in consecutive elements
272 *> of the array. If the eigenvalues cannot be reordered to
273 *> compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
274 *> when the true value would be very small anyway.
275 *> If SENSE = 'N' or 'E', RCONDV is not referenced.
280 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
281 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
287 *> The dimension of the array WORK. LWORK >= max(1,2*N).
288 *> If SENSE = 'E', LWORK >= max(1,4*N).
289 *> If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
291 *> If LWORK = -1, then a workspace query is assumed; the routine
292 *> only calculates the optimal size of the WORK array, returns
293 *> this value as the first entry of the WORK array, and no error
294 *> message related to LWORK is issued by XERBLA.
299 *> RWORK is REAL array, dimension (lrwork)
300 *> lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
301 *> and at least max(1,2*N) otherwise.
307 *> IWORK is INTEGER array, dimension (N+2)
308 *> If SENSE = 'E', IWORK is not referenced.
313 *> BWORK is LOGICAL array, dimension (N)
314 *> If SENSE = 'N', BWORK is not referenced.
320 *> = 0: successful exit
321 *> < 0: if INFO = -i, the i-th argument had an illegal value.
323 *> The QZ iteration failed. No eigenvectors have been
324 *> calculated, but ALPHA(j) and BETA(j) should be correct
325 *> for j=INFO+1,...,N.
326 *> > N: =N+1: other than QZ iteration failed in CHGEQZ.
327 *> =N+2: error return from CTGEVC.
333 *> \author Univ. of Tennessee
334 *> \author Univ. of California Berkeley
335 *> \author Univ. of Colorado Denver
340 *> \ingroup complexGEeigen
342 *> \par Further Details:
343 * =====================
347 *> Balancing a matrix pair (A,B) includes, first, permuting rows and
348 *> columns to isolate eigenvalues, second, applying diagonal similarity
349 *> transformation to the rows and columns to make the rows and columns
350 *> as close in norm as possible. The computed reciprocal condition
351 *> numbers correspond to the balanced matrix. Permuting rows and columns
352 *> will not change the condition numbers (in exact arithmetic) but
353 *> diagonal scaling will. For further explanation of balancing, see
354 *> section 4.11.1.2 of LAPACK Users' Guide.
356 *> An approximate error bound on the chordal distance between the i-th
357 *> computed generalized eigenvalue w and the corresponding exact
358 *> eigenvalue lambda is
360 *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
362 *> An approximate error bound for the angle between the i-th computed
363 *> eigenvector VL(i) or VR(i) is given by
365 *> EPS * norm(ABNRM, BBNRM) / DIF(i).
367 *> For further explanation of the reciprocal condition numbers RCONDE
368 *> and RCONDV, see section 4.11 of LAPACK User's Guide.
371 * =====================================================================
372 SUBROUTINE CGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
373 $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
374 $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
375 $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
377 * -- LAPACK driver routine (version 3.4.1) --
378 * -- LAPACK is a software package provided by Univ. of Tennessee, --
379 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
382 * .. Scalar Arguments ..
383 CHARACTER BALANC, JOBVL, JOBVR, SENSE
384 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
387 * .. Array Arguments ..
390 REAL LSCALE( * ), RCONDE( * ), RCONDV( * ),
391 $ RSCALE( * ), RWORK( * )
392 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
393 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
397 * =====================================================================
401 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
403 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
404 $ CONE = ( 1.0E+0, 0.0E+0 ) )
406 * .. Local Scalars ..
407 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
408 $ WANTSB, WANTSE, WANTSN, WANTSV
410 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
411 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
412 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
419 * .. External Subroutines ..
420 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
421 $ CLASCL, CLASET, CTGEVC, CTGSNA, CUNGQR, CUNMQR,
422 $ SLABAD, SLASCL, XERBLA
424 * .. External Functions ..
428 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH
430 * .. Intrinsic Functions ..
431 INTRINSIC ABS, AIMAG, MAX, REAL, SQRT
433 * .. Statement Functions ..
436 * .. Statement Function definitions ..
437 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
439 * .. Executable Statements ..
441 * Decode the input arguments
443 IF( LSAME( JOBVL, 'N' ) ) THEN
446 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
454 IF( LSAME( JOBVR, 'N' ) ) THEN
457 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
466 NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
467 WANTSN = LSAME( SENSE, 'N' )
468 WANTSE = LSAME( SENSE, 'E' )
469 WANTSV = LSAME( SENSE, 'V' )
470 WANTSB = LSAME( SENSE, 'B' )
472 * Test the input arguments
475 LQUERY = ( LWORK.EQ.-1 )
476 IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
477 $ LSAME( BALANC, 'B' ) ) ) THEN
479 ELSE IF( IJOBVL.LE.0 ) THEN
481 ELSE IF( IJOBVR.LE.0 ) THEN
483 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
486 ELSE IF( N.LT.0 ) THEN
488 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
490 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
492 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
494 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
499 * (Note: Comments in the code beginning "Workspace:" describe the
500 * minimal amount of workspace needed at that point in the code,
501 * as well as the preferred amount for good performance.
502 * NB refers to the optimal block size for the immediately
503 * following subroutine, as returned by ILAENV. The workspace is
504 * computed assuming ILO = 1 and IHI = N, the worst case.)
514 ELSE IF( WANTSV .OR. WANTSB ) THEN
515 MINWRK = 2*N*( N + 1)
518 MAXWRK = MAX( MAXWRK,
519 $ N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) )
520 MAXWRK = MAX( MAXWRK,
521 $ N + N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) )
523 MAXWRK = MAX( MAXWRK, N +
524 $ N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, 0 ) )
529 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
535 CALL XERBLA( 'CGGEVX', -INFO )
537 ELSE IF( LQUERY ) THEN
541 * Quick return if possible
546 * Get machine constants
549 SMLNUM = SLAMCH( 'S' )
550 BIGNUM = ONE / SMLNUM
551 CALL SLABAD( SMLNUM, BIGNUM )
552 SMLNUM = SQRT( SMLNUM ) / EPS
553 BIGNUM = ONE / SMLNUM
555 * Scale A if max element outside range [SMLNUM,BIGNUM]
557 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
559 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
562 ELSE IF( ANRM.GT.BIGNUM ) THEN
567 $ CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
569 * Scale B if max element outside range [SMLNUM,BIGNUM]
571 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
573 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
576 ELSE IF( BNRM.GT.BIGNUM ) THEN
581 $ CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
583 * Permute and/or balance the matrix pair (A,B)
584 * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
586 CALL CGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
589 * Compute ABNRM and BBNRM
591 ABNRM = CLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
594 CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
599 BBNRM = CLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
602 CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
607 * Reduce B to triangular form (QR decomposition of B)
608 * (Complex Workspace: need N, prefer N*NB )
610 IROWS = IHI + 1 - ILO
611 IF( ILV .OR. .NOT.WANTSN ) THEN
618 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
619 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
621 * Apply the unitary transformation to A
622 * (Complex Workspace: need N, prefer N*NB)
624 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
625 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
626 $ LWORK+1-IWRK, IERR )
628 * Initialize VL and/or VR
629 * (Workspace: need N, prefer N*NB)
632 CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
633 IF( IROWS.GT.1 ) THEN
634 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
635 $ VL( ILO+1, ILO ), LDVL )
637 CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
638 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
642 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
644 * Reduce to generalized Hessenberg form
645 * (Workspace: none needed)
647 IF( ILV .OR. .NOT.WANTSN ) THEN
649 * Eigenvectors requested -- work on whole matrix.
651 CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
652 $ LDVL, VR, LDVR, IERR )
654 CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
655 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
658 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
659 * Schur forms and Schur vectors)
660 * (Complex Workspace: need N)
661 * (Real Workspace: need N)
664 IF( ILV .OR. .NOT.WANTSN ) THEN
670 CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
671 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
672 $ LWORK+1-IWRK, RWORK, IERR )
674 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
676 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
684 * Compute Eigenvectors and estimate condition numbers if desired
685 * CTGEVC: (Complex Workspace: need 2*N )
686 * (Real Workspace: need 2*N )
687 * CTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
688 * (Integer Workspace: need N+2 )
690 IF( ILV .OR. .NOT.WANTSN ) THEN
702 CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
703 $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
711 IF( .NOT.WANTSN ) THEN
713 * compute eigenvectors (STGEVC) and estimate condition
714 * numbers (STGSNA). Note that the definition of the condition
715 * number is not invariant under transformation (u,v) to
716 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
717 * Schur form (S,T), Q and Z are orthogonal matrices. In order
718 * to avoid using extra 2*N*N workspace, we have to
719 * re-calculate eigenvectors and estimate the condition numbers
732 IF( WANTSE .OR. WANTSB ) THEN
733 CALL CTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
734 $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
735 $ WORK( IWRK1 ), RWORK, IERR )
742 CALL CTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
743 $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
744 $ RCONDV( I ), 1, M, WORK( IWRK1 ),
745 $ LWORK-IWRK1+1, IWORK, IERR )
751 * Undo balancing on VL and VR and normalization
752 * (Workspace: none needed)
755 CALL CGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
761 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
767 VL( JR, JC ) = VL( JR, JC )*TEMP
773 CALL CGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
778 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
784 VR( JR, JC ) = VR( JR, JC )*TEMP
789 * Undo scaling if necessary
794 $ CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
797 $ CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )