1 SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
2 $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
3 $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
4 $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 CHARACTER BALANC, JOBVL, JOBVR, SENSE
13 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
14 DOUBLE PRECISION ABNRM, BBNRM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
20 $ RSCALE( * ), RWORK( * )
21 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
22 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
23 $ WORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
30 * (A,B) the generalized eigenvalues, and optionally, the left and/or
31 * right generalized eigenvectors.
32 *
33 * Optionally, it also computes a balancing transformation to improve
34 * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
35 * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
36 * the eigenvalues (RCONDE), and reciprocal condition numbers for the
37 * right eigenvectors (RCONDV).
38 *
39 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
40 * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
41 * singular. It is usually represented as the pair (alpha,beta), as
42 * there is a reasonable interpretation for beta=0, and even for both
43 * being zero.
44 *
45 * The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
46 * of (A,B) satisfies
47 * A * v(j) = lambda(j) * B * v(j) .
48 * The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
49 * of (A,B) satisfies
50 * u(j)**H * A = lambda(j) * u(j)**H * B.
51 * where u(j)**H is the conjugate-transpose of u(j).
52 *
53 *
54 * Arguments
55 * =========
56 *
57 * BALANC (input) CHARACTER*1
58 * Specifies the balance option to be performed:
59 * = 'N': do not diagonally scale or permute;
60 * = 'P': permute only;
61 * = 'S': scale only;
62 * = 'B': both permute and scale.
63 * Computed reciprocal condition numbers will be for the
64 * matrices after permuting and/or balancing. Permuting does
65 * not change condition numbers (in exact arithmetic), but
66 * balancing does.
67 *
68 * JOBVL (input) CHARACTER*1
69 * = 'N': do not compute the left generalized eigenvectors;
70 * = 'V': compute the left generalized eigenvectors.
71 *
72 * JOBVR (input) CHARACTER*1
73 * = 'N': do not compute the right generalized eigenvectors;
74 * = 'V': compute the right generalized eigenvectors.
75 *
76 * SENSE (input) CHARACTER*1
77 * Determines which reciprocal condition numbers are computed.
78 * = 'N': none are computed;
79 * = 'E': computed for eigenvalues only;
80 * = 'V': computed for eigenvectors only;
81 * = 'B': computed for eigenvalues and eigenvectors.
82 *
83 * N (input) INTEGER
84 * The order of the matrices A, B, VL, and VR. N >= 0.
85 *
86 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
87 * On entry, the matrix A in the pair (A,B).
88 * On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
89 * or both, then A contains the first part of the complex Schur
90 * form of the "balanced" versions of the input A and B.
91 *
92 * LDA (input) INTEGER
93 * The leading dimension of A. LDA >= max(1,N).
94 *
95 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
96 * On entry, the matrix B in the pair (A,B).
97 * On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
98 * or both, then B contains the second part of the complex
99 * Schur form of the "balanced" versions of the input A and B.
100 *
101 * LDB (input) INTEGER
102 * The leading dimension of B. LDB >= max(1,N).
103 *
104 * ALPHA (output) COMPLEX*16 array, dimension (N)
105 * BETA (output) COMPLEX*16 array, dimension (N)
106 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
107 * eigenvalues.
108 *
109 * Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
110 * underflow, and BETA(j) may even be zero. Thus, the user
111 * should avoid naively computing the ratio ALPHA/BETA.
112 * However, ALPHA will be always less than and usually
113 * comparable with norm(A) in magnitude, and BETA always less
114 * than and usually comparable with norm(B).
115 *
116 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
117 * If JOBVL = 'V', the left generalized eigenvectors u(j) are
118 * stored one after another in the columns of VL, in the same
119 * order as their eigenvalues.
120 * Each eigenvector will be scaled so the largest component
121 * will have abs(real part) + abs(imag. part) = 1.
122 * Not referenced if JOBVL = 'N'.
123 *
124 * LDVL (input) INTEGER
125 * The leading dimension of the matrix VL. LDVL >= 1, and
126 * if JOBVL = 'V', LDVL >= N.
127 *
128 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
129 * If JOBVR = 'V', the right generalized eigenvectors v(j) are
130 * stored one after another in the columns of VR, in the same
131 * order as their eigenvalues.
132 * Each eigenvector will be scaled so the largest component
133 * will have abs(real part) + abs(imag. part) = 1.
134 * Not referenced if JOBVR = 'N'.
135 *
136 * LDVR (input) INTEGER
137 * The leading dimension of the matrix VR. LDVR >= 1, and
138 * if JOBVR = 'V', LDVR >= N.
139 *
140 * ILO (output) INTEGER
141 * IHI (output) INTEGER
142 * ILO and IHI are integer values such that on exit
143 * A(i,j) = 0 and B(i,j) = 0 if i > j and
144 * j = 1,...,ILO-1 or i = IHI+1,...,N.
145 * If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
146 *
147 * LSCALE (output) DOUBLE PRECISION array, dimension (N)
148 * Details of the permutations and scaling factors applied
149 * to the left side of A and B. If PL(j) is the index of the
150 * row interchanged with row j, and DL(j) is the scaling
151 * factor applied to row j, then
152 * LSCALE(j) = PL(j) for j = 1,...,ILO-1
153 * = DL(j) for j = ILO,...,IHI
154 * = PL(j) for j = IHI+1,...,N.
155 * The order in which the interchanges are made is N to IHI+1,
156 * then 1 to ILO-1.
157 *
158 * RSCALE (output) DOUBLE PRECISION array, dimension (N)
159 * Details of the permutations and scaling factors applied
160 * to the right side of A and B. If PR(j) is the index of the
161 * column interchanged with column j, and DR(j) is the scaling
162 * factor applied to column j, then
163 * RSCALE(j) = PR(j) for j = 1,...,ILO-1
164 * = DR(j) for j = ILO,...,IHI
165 * = PR(j) for j = IHI+1,...,N
166 * The order in which the interchanges are made is N to IHI+1,
167 * then 1 to ILO-1.
168 *
169 * ABNRM (output) DOUBLE PRECISION
170 * The one-norm of the balanced matrix A.
171 *
172 * BBNRM (output) DOUBLE PRECISION
173 * The one-norm of the balanced matrix B.
174 *
175 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
176 * If SENSE = 'E' or 'B', the reciprocal condition numbers of
177 * the eigenvalues, stored in consecutive elements of the array.
178 * If SENSE = 'N' or 'V', RCONDE is not referenced.
179 *
180 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
181 * If JOB = 'V' or 'B', the estimated reciprocal condition
182 * numbers of the eigenvectors, stored in consecutive elements
183 * of the array. If the eigenvalues cannot be reordered to
184 * compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
185 * when the true value would be very small anyway.
186 * If SENSE = 'N' or 'E', RCONDV is not referenced.
187 *
188 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
189 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
190 *
191 * LWORK (input) INTEGER
192 * The dimension of the array WORK. LWORK >= max(1,2*N).
193 * If SENSE = 'E', LWORK >= max(1,4*N).
194 * If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
195 *
196 * If LWORK = -1, then a workspace query is assumed; the routine
197 * only calculates the optimal size of the WORK array, returns
198 * this value as the first entry of the WORK array, and no error
199 * message related to LWORK is issued by XERBLA.
200 *
201 * RWORK (workspace) REAL array, dimension (lrwork)
202 * lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
203 * and at least max(1,2*N) otherwise.
204 * Real workspace.
205 *
206 * IWORK (workspace) INTEGER array, dimension (N+2)
207 * If SENSE = 'E', IWORK is not referenced.
208 *
209 * BWORK (workspace) LOGICAL array, dimension (N)
210 * If SENSE = 'N', BWORK is not referenced.
211 *
212 * INFO (output) INTEGER
213 * = 0: successful exit
214 * < 0: if INFO = -i, the i-th argument had an illegal value.
215 * = 1,...,N:
216 * The QZ iteration failed. No eigenvectors have been
217 * calculated, but ALPHA(j) and BETA(j) should be correct
218 * for j=INFO+1,...,N.
219 * > N: =N+1: other than QZ iteration failed in ZHGEQZ.
220 * =N+2: error return from ZTGEVC.
221 *
222 * Further Details
223 * ===============
224 *
225 * Balancing a matrix pair (A,B) includes, first, permuting rows and
226 * columns to isolate eigenvalues, second, applying diagonal similarity
227 * transformation to the rows and columns to make the rows and columns
228 * as close in norm as possible. The computed reciprocal condition
229 * numbers correspond to the balanced matrix. Permuting rows and columns
230 * will not change the condition numbers (in exact arithmetic) but
231 * diagonal scaling will. For further explanation of balancing, see
232 * section 4.11.1.2 of LAPACK Users' Guide.
233 *
234 * An approximate error bound on the chordal distance between the i-th
235 * computed generalized eigenvalue w and the corresponding exact
236 * eigenvalue lambda is
237 *
238 * chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
239 *
240 * An approximate error bound for the angle between the i-th computed
241 * eigenvector VL(i) or VR(i) is given by
242 *
243 * EPS * norm(ABNRM, BBNRM) / DIF(i).
244 *
245 * For further explanation of the reciprocal condition numbers RCONDE
246 * and RCONDV, see section 4.11 of LAPACK User's Guide.
247 *
248 * .. Parameters ..
249 DOUBLE PRECISION ZERO, ONE
250 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
251 COMPLEX*16 CZERO, CONE
252 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
253 $ CONE = ( 1.0D+0, 0.0D+0 ) )
254 * ..
255 * .. Local Scalars ..
256 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
257 $ WANTSB, WANTSE, WANTSN, WANTSV
258 CHARACTER CHTEMP
259 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
260 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
261 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
262 $ SMLNUM, TEMP
263 COMPLEX*16 X
264 * ..
265 * .. Local Arrays ..
266 LOGICAL LDUMMA( 1 )
267 * ..
268 * .. External Subroutines ..
269 EXTERNAL DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
270 $ ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
271 $ ZTGSNA, ZUNGQR, ZUNMQR
272 * ..
273 * .. External Functions ..
274 LOGICAL LSAME
275 INTEGER ILAENV
276 DOUBLE PRECISION DLAMCH, ZLANGE
277 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
278 * ..
279 * .. Intrinsic Functions ..
280 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
281 * ..
282 * .. Statement Functions ..
283 DOUBLE PRECISION ABS1
284 * ..
285 * .. Statement Function definitions ..
286 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
287 * ..
288 * .. Executable Statements ..
289 *
290 * Decode the input arguments
291 *
292 IF( LSAME( JOBVL, 'N' ) ) THEN
293 IJOBVL = 1
294 ILVL = .FALSE.
295 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
296 IJOBVL = 2
297 ILVL = .TRUE.
298 ELSE
299 IJOBVL = -1
300 ILVL = .FALSE.
301 END IF
302 *
303 IF( LSAME( JOBVR, 'N' ) ) THEN
304 IJOBVR = 1
305 ILVR = .FALSE.
306 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
307 IJOBVR = 2
308 ILVR = .TRUE.
309 ELSE
310 IJOBVR = -1
311 ILVR = .FALSE.
312 END IF
313 ILV = ILVL .OR. ILVR
314 *
315 NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
316 WANTSN = LSAME( SENSE, 'N' )
317 WANTSE = LSAME( SENSE, 'E' )
318 WANTSV = LSAME( SENSE, 'V' )
319 WANTSB = LSAME( SENSE, 'B' )
320 *
321 * Test the input arguments
322 *
323 INFO = 0
324 LQUERY = ( LWORK.EQ.-1 )
325 IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
326 $ LSAME( BALANC, 'B' ) ) ) THEN
327 INFO = -1
328 ELSE IF( IJOBVL.LE.0 ) THEN
329 INFO = -2
330 ELSE IF( IJOBVR.LE.0 ) THEN
331 INFO = -3
332 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
333 $ THEN
334 INFO = -4
335 ELSE IF( N.LT.0 ) THEN
336 INFO = -5
337 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338 INFO = -7
339 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340 INFO = -9
341 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
342 INFO = -13
343 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
344 INFO = -15
345 END IF
346 *
347 * Compute workspace
348 * (Note: Comments in the code beginning "Workspace:" describe the
349 * minimal amount of workspace needed at that point in the code,
350 * as well as the preferred amount for good performance.
351 * NB refers to the optimal block size for the immediately
352 * following subroutine, as returned by ILAENV. The workspace is
353 * computed assuming ILO = 1 and IHI = N, the worst case.)
354 *
355 IF( INFO.EQ.0 ) THEN
356 IF( N.EQ.0 ) THEN
357 MINWRK = 1
358 MAXWRK = 1
359 ELSE
360 MINWRK = 2*N
361 IF( WANTSE ) THEN
362 MINWRK = 4*N
363 ELSE IF( WANTSV .OR. WANTSB ) THEN
364 MINWRK = 2*N*( N + 1)
365 END IF
366 MAXWRK = MINWRK
367 MAXWRK = MAX( MAXWRK,
368 $ N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
369 MAXWRK = MAX( MAXWRK,
370 $ N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
371 IF( ILVL ) THEN
372 MAXWRK = MAX( MAXWRK, N +
373 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
374 END IF
375 END IF
376 WORK( 1 ) = MAXWRK
377 *
378 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
379 INFO = -25
380 END IF
381 END IF
382 *
383 IF( INFO.NE.0 ) THEN
384 CALL XERBLA( 'ZGGEVX', -INFO )
385 RETURN
386 ELSE IF( LQUERY ) THEN
387 RETURN
388 END IF
389 *
390 * Quick return if possible
391 *
392 IF( N.EQ.0 )
393 $ RETURN
394 *
395 * Get machine constants
396 *
397 EPS = DLAMCH( 'P' )
398 SMLNUM = DLAMCH( 'S' )
399 BIGNUM = ONE / SMLNUM
400 CALL DLABAD( SMLNUM, BIGNUM )
401 SMLNUM = SQRT( SMLNUM ) / EPS
402 BIGNUM = ONE / SMLNUM
403 *
404 * Scale A if max element outside range [SMLNUM,BIGNUM]
405 *
406 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
407 ILASCL = .FALSE.
408 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
409 ANRMTO = SMLNUM
410 ILASCL = .TRUE.
411 ELSE IF( ANRM.GT.BIGNUM ) THEN
412 ANRMTO = BIGNUM
413 ILASCL = .TRUE.
414 END IF
415 IF( ILASCL )
416 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
417 *
418 * Scale B if max element outside range [SMLNUM,BIGNUM]
419 *
420 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
421 ILBSCL = .FALSE.
422 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
423 BNRMTO = SMLNUM
424 ILBSCL = .TRUE.
425 ELSE IF( BNRM.GT.BIGNUM ) THEN
426 BNRMTO = BIGNUM
427 ILBSCL = .TRUE.
428 END IF
429 IF( ILBSCL )
430 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
431 *
432 * Permute and/or balance the matrix pair (A,B)
433 * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
434 *
435 CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
436 $ RWORK, IERR )
437 *
438 * Compute ABNRM and BBNRM
439 *
440 ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
441 IF( ILASCL ) THEN
442 RWORK( 1 ) = ABNRM
443 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
444 $ IERR )
445 ABNRM = RWORK( 1 )
446 END IF
447 *
448 BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
449 IF( ILBSCL ) THEN
450 RWORK( 1 ) = BBNRM
451 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
452 $ IERR )
453 BBNRM = RWORK( 1 )
454 END IF
455 *
456 * Reduce B to triangular form (QR decomposition of B)
457 * (Complex Workspace: need N, prefer N*NB )
458 *
459 IROWS = IHI + 1 - ILO
460 IF( ILV .OR. .NOT.WANTSN ) THEN
461 ICOLS = N + 1 - ILO
462 ELSE
463 ICOLS = IROWS
464 END IF
465 ITAU = 1
466 IWRK = ITAU + IROWS
467 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
468 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
469 *
470 * Apply the unitary transformation to A
471 * (Complex Workspace: need N, prefer N*NB)
472 *
473 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
474 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
475 $ LWORK+1-IWRK, IERR )
476 *
477 * Initialize VL and/or VR
478 * (Workspace: need N, prefer N*NB)
479 *
480 IF( ILVL ) THEN
481 CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
482 IF( IROWS.GT.1 ) THEN
483 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
484 $ VL( ILO+1, ILO ), LDVL )
485 END IF
486 CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
487 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
488 END IF
489 *
490 IF( ILVR )
491 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
492 *
493 * Reduce to generalized Hessenberg form
494 * (Workspace: none needed)
495 *
496 IF( ILV .OR. .NOT.WANTSN ) THEN
497 *
498 * Eigenvectors requested -- work on whole matrix.
499 *
500 CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
501 $ LDVL, VR, LDVR, IERR )
502 ELSE
503 CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
504 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
505 END IF
506 *
507 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
508 * Schur forms and Schur vectors)
509 * (Complex Workspace: need N)
510 * (Real Workspace: need N)
511 *
512 IWRK = ITAU
513 IF( ILV .OR. .NOT.WANTSN ) THEN
514 CHTEMP = 'S'
515 ELSE
516 CHTEMP = 'E'
517 END IF
518 *
519 CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
520 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
521 $ LWORK+1-IWRK, RWORK, IERR )
522 IF( IERR.NE.0 ) THEN
523 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
524 INFO = IERR
525 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
526 INFO = IERR - N
527 ELSE
528 INFO = N + 1
529 END IF
530 GO TO 90
531 END IF
532 *
533 * Compute Eigenvectors and estimate condition numbers if desired
534 * ZTGEVC: (Complex Workspace: need 2*N )
535 * (Real Workspace: need 2*N )
536 * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
537 * (Integer Workspace: need N+2 )
538 *
539 IF( ILV .OR. .NOT.WANTSN ) THEN
540 IF( ILV ) THEN
541 IF( ILVL ) THEN
542 IF( ILVR ) THEN
543 CHTEMP = 'B'
544 ELSE
545 CHTEMP = 'L'
546 END IF
547 ELSE
548 CHTEMP = 'R'
549 END IF
550 *
551 CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
552 $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
553 $ IERR )
554 IF( IERR.NE.0 ) THEN
555 INFO = N + 2
556 GO TO 90
557 END IF
558 END IF
559 *
560 IF( .NOT.WANTSN ) THEN
561 *
562 * compute eigenvectors (DTGEVC) and estimate condition
563 * numbers (DTGSNA). Note that the definition of the condition
564 * number is not invariant under transformation (u,v) to
565 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
566 * Schur form (S,T), Q and Z are orthogonal matrices. In order
567 * to avoid using extra 2*N*N workspace, we have to
568 * re-calculate eigenvectors and estimate the condition numbers
569 * one at a time.
570 *
571 DO 20 I = 1, N
572 *
573 DO 10 J = 1, N
574 BWORK( J ) = .FALSE.
575 10 CONTINUE
576 BWORK( I ) = .TRUE.
577 *
578 IWRK = N + 1
579 IWRK1 = IWRK + N
580 *
581 IF( WANTSE .OR. WANTSB ) THEN
582 CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
583 $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
584 $ WORK( IWRK1 ), RWORK, IERR )
585 IF( IERR.NE.0 ) THEN
586 INFO = N + 2
587 GO TO 90
588 END IF
589 END IF
590 *
591 CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
592 $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
593 $ RCONDV( I ), 1, M, WORK( IWRK1 ),
594 $ LWORK-IWRK1+1, IWORK, IERR )
595 *
596 20 CONTINUE
597 END IF
598 END IF
599 *
600 * Undo balancing on VL and VR and normalization
601 * (Workspace: none needed)
602 *
603 IF( ILVL ) THEN
604 CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
605 $ LDVL, IERR )
606 *
607 DO 50 JC = 1, N
608 TEMP = ZERO
609 DO 30 JR = 1, N
610 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
611 30 CONTINUE
612 IF( TEMP.LT.SMLNUM )
613 $ GO TO 50
614 TEMP = ONE / TEMP
615 DO 40 JR = 1, N
616 VL( JR, JC ) = VL( JR, JC )*TEMP
617 40 CONTINUE
618 50 CONTINUE
619 END IF
620 *
621 IF( ILVR ) THEN
622 CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
623 $ LDVR, IERR )
624 DO 80 JC = 1, N
625 TEMP = ZERO
626 DO 60 JR = 1, N
627 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
628 60 CONTINUE
629 IF( TEMP.LT.SMLNUM )
630 $ GO TO 80
631 TEMP = ONE / TEMP
632 DO 70 JR = 1, N
633 VR( JR, JC ) = VR( JR, JC )*TEMP
634 70 CONTINUE
635 80 CONTINUE
636 END IF
637 *
638 * Undo scaling if necessary
639 *
640 IF( ILASCL )
641 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
642 *
643 IF( ILBSCL )
644 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
645 *
646 90 CONTINUE
647 WORK( 1 ) = MAXWRK
648 *
649 RETURN
650 *
651 * End of ZGGEVX
652 *
653 END
2 $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
3 $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
4 $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 CHARACTER BALANC, JOBVL, JOBVR, SENSE
13 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
14 DOUBLE PRECISION ABNRM, BBNRM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
20 $ RSCALE( * ), RWORK( * )
21 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
22 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
23 $ WORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
30 * (A,B) the generalized eigenvalues, and optionally, the left and/or
31 * right generalized eigenvectors.
32 *
33 * Optionally, it also computes a balancing transformation to improve
34 * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
35 * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
36 * the eigenvalues (RCONDE), and reciprocal condition numbers for the
37 * right eigenvectors (RCONDV).
38 *
39 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
40 * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
41 * singular. It is usually represented as the pair (alpha,beta), as
42 * there is a reasonable interpretation for beta=0, and even for both
43 * being zero.
44 *
45 * The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
46 * of (A,B) satisfies
47 * A * v(j) = lambda(j) * B * v(j) .
48 * The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
49 * of (A,B) satisfies
50 * u(j)**H * A = lambda(j) * u(j)**H * B.
51 * where u(j)**H is the conjugate-transpose of u(j).
52 *
53 *
54 * Arguments
55 * =========
56 *
57 * BALANC (input) CHARACTER*1
58 * Specifies the balance option to be performed:
59 * = 'N': do not diagonally scale or permute;
60 * = 'P': permute only;
61 * = 'S': scale only;
62 * = 'B': both permute and scale.
63 * Computed reciprocal condition numbers will be for the
64 * matrices after permuting and/or balancing. Permuting does
65 * not change condition numbers (in exact arithmetic), but
66 * balancing does.
67 *
68 * JOBVL (input) CHARACTER*1
69 * = 'N': do not compute the left generalized eigenvectors;
70 * = 'V': compute the left generalized eigenvectors.
71 *
72 * JOBVR (input) CHARACTER*1
73 * = 'N': do not compute the right generalized eigenvectors;
74 * = 'V': compute the right generalized eigenvectors.
75 *
76 * SENSE (input) CHARACTER*1
77 * Determines which reciprocal condition numbers are computed.
78 * = 'N': none are computed;
79 * = 'E': computed for eigenvalues only;
80 * = 'V': computed for eigenvectors only;
81 * = 'B': computed for eigenvalues and eigenvectors.
82 *
83 * N (input) INTEGER
84 * The order of the matrices A, B, VL, and VR. N >= 0.
85 *
86 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
87 * On entry, the matrix A in the pair (A,B).
88 * On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
89 * or both, then A contains the first part of the complex Schur
90 * form of the "balanced" versions of the input A and B.
91 *
92 * LDA (input) INTEGER
93 * The leading dimension of A. LDA >= max(1,N).
94 *
95 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
96 * On entry, the matrix B in the pair (A,B).
97 * On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
98 * or both, then B contains the second part of the complex
99 * Schur form of the "balanced" versions of the input A and B.
100 *
101 * LDB (input) INTEGER
102 * The leading dimension of B. LDB >= max(1,N).
103 *
104 * ALPHA (output) COMPLEX*16 array, dimension (N)
105 * BETA (output) COMPLEX*16 array, dimension (N)
106 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
107 * eigenvalues.
108 *
109 * Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
110 * underflow, and BETA(j) may even be zero. Thus, the user
111 * should avoid naively computing the ratio ALPHA/BETA.
112 * However, ALPHA will be always less than and usually
113 * comparable with norm(A) in magnitude, and BETA always less
114 * than and usually comparable with norm(B).
115 *
116 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
117 * If JOBVL = 'V', the left generalized eigenvectors u(j) are
118 * stored one after another in the columns of VL, in the same
119 * order as their eigenvalues.
120 * Each eigenvector will be scaled so the largest component
121 * will have abs(real part) + abs(imag. part) = 1.
122 * Not referenced if JOBVL = 'N'.
123 *
124 * LDVL (input) INTEGER
125 * The leading dimension of the matrix VL. LDVL >= 1, and
126 * if JOBVL = 'V', LDVL >= N.
127 *
128 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
129 * If JOBVR = 'V', the right generalized eigenvectors v(j) are
130 * stored one after another in the columns of VR, in the same
131 * order as their eigenvalues.
132 * Each eigenvector will be scaled so the largest component
133 * will have abs(real part) + abs(imag. part) = 1.
134 * Not referenced if JOBVR = 'N'.
135 *
136 * LDVR (input) INTEGER
137 * The leading dimension of the matrix VR. LDVR >= 1, and
138 * if JOBVR = 'V', LDVR >= N.
139 *
140 * ILO (output) INTEGER
141 * IHI (output) INTEGER
142 * ILO and IHI are integer values such that on exit
143 * A(i,j) = 0 and B(i,j) = 0 if i > j and
144 * j = 1,...,ILO-1 or i = IHI+1,...,N.
145 * If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
146 *
147 * LSCALE (output) DOUBLE PRECISION array, dimension (N)
148 * Details of the permutations and scaling factors applied
149 * to the left side of A and B. If PL(j) is the index of the
150 * row interchanged with row j, and DL(j) is the scaling
151 * factor applied to row j, then
152 * LSCALE(j) = PL(j) for j = 1,...,ILO-1
153 * = DL(j) for j = ILO,...,IHI
154 * = PL(j) for j = IHI+1,...,N.
155 * The order in which the interchanges are made is N to IHI+1,
156 * then 1 to ILO-1.
157 *
158 * RSCALE (output) DOUBLE PRECISION array, dimension (N)
159 * Details of the permutations and scaling factors applied
160 * to the right side of A and B. If PR(j) is the index of the
161 * column interchanged with column j, and DR(j) is the scaling
162 * factor applied to column j, then
163 * RSCALE(j) = PR(j) for j = 1,...,ILO-1
164 * = DR(j) for j = ILO,...,IHI
165 * = PR(j) for j = IHI+1,...,N
166 * The order in which the interchanges are made is N to IHI+1,
167 * then 1 to ILO-1.
168 *
169 * ABNRM (output) DOUBLE PRECISION
170 * The one-norm of the balanced matrix A.
171 *
172 * BBNRM (output) DOUBLE PRECISION
173 * The one-norm of the balanced matrix B.
174 *
175 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
176 * If SENSE = 'E' or 'B', the reciprocal condition numbers of
177 * the eigenvalues, stored in consecutive elements of the array.
178 * If SENSE = 'N' or 'V', RCONDE is not referenced.
179 *
180 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
181 * If JOB = 'V' or 'B', the estimated reciprocal condition
182 * numbers of the eigenvectors, stored in consecutive elements
183 * of the array. If the eigenvalues cannot be reordered to
184 * compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
185 * when the true value would be very small anyway.
186 * If SENSE = 'N' or 'E', RCONDV is not referenced.
187 *
188 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
189 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
190 *
191 * LWORK (input) INTEGER
192 * The dimension of the array WORK. LWORK >= max(1,2*N).
193 * If SENSE = 'E', LWORK >= max(1,4*N).
194 * If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
195 *
196 * If LWORK = -1, then a workspace query is assumed; the routine
197 * only calculates the optimal size of the WORK array, returns
198 * this value as the first entry of the WORK array, and no error
199 * message related to LWORK is issued by XERBLA.
200 *
201 * RWORK (workspace) REAL array, dimension (lrwork)
202 * lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
203 * and at least max(1,2*N) otherwise.
204 * Real workspace.
205 *
206 * IWORK (workspace) INTEGER array, dimension (N+2)
207 * If SENSE = 'E', IWORK is not referenced.
208 *
209 * BWORK (workspace) LOGICAL array, dimension (N)
210 * If SENSE = 'N', BWORK is not referenced.
211 *
212 * INFO (output) INTEGER
213 * = 0: successful exit
214 * < 0: if INFO = -i, the i-th argument had an illegal value.
215 * = 1,...,N:
216 * The QZ iteration failed. No eigenvectors have been
217 * calculated, but ALPHA(j) and BETA(j) should be correct
218 * for j=INFO+1,...,N.
219 * > N: =N+1: other than QZ iteration failed in ZHGEQZ.
220 * =N+2: error return from ZTGEVC.
221 *
222 * Further Details
223 * ===============
224 *
225 * Balancing a matrix pair (A,B) includes, first, permuting rows and
226 * columns to isolate eigenvalues, second, applying diagonal similarity
227 * transformation to the rows and columns to make the rows and columns
228 * as close in norm as possible. The computed reciprocal condition
229 * numbers correspond to the balanced matrix. Permuting rows and columns
230 * will not change the condition numbers (in exact arithmetic) but
231 * diagonal scaling will. For further explanation of balancing, see
232 * section 4.11.1.2 of LAPACK Users' Guide.
233 *
234 * An approximate error bound on the chordal distance between the i-th
235 * computed generalized eigenvalue w and the corresponding exact
236 * eigenvalue lambda is
237 *
238 * chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
239 *
240 * An approximate error bound for the angle between the i-th computed
241 * eigenvector VL(i) or VR(i) is given by
242 *
243 * EPS * norm(ABNRM, BBNRM) / DIF(i).
244 *
245 * For further explanation of the reciprocal condition numbers RCONDE
246 * and RCONDV, see section 4.11 of LAPACK User's Guide.
247 *
248 * .. Parameters ..
249 DOUBLE PRECISION ZERO, ONE
250 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
251 COMPLEX*16 CZERO, CONE
252 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
253 $ CONE = ( 1.0D+0, 0.0D+0 ) )
254 * ..
255 * .. Local Scalars ..
256 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
257 $ WANTSB, WANTSE, WANTSN, WANTSV
258 CHARACTER CHTEMP
259 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
260 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
261 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
262 $ SMLNUM, TEMP
263 COMPLEX*16 X
264 * ..
265 * .. Local Arrays ..
266 LOGICAL LDUMMA( 1 )
267 * ..
268 * .. External Subroutines ..
269 EXTERNAL DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
270 $ ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
271 $ ZTGSNA, ZUNGQR, ZUNMQR
272 * ..
273 * .. External Functions ..
274 LOGICAL LSAME
275 INTEGER ILAENV
276 DOUBLE PRECISION DLAMCH, ZLANGE
277 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
278 * ..
279 * .. Intrinsic Functions ..
280 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
281 * ..
282 * .. Statement Functions ..
283 DOUBLE PRECISION ABS1
284 * ..
285 * .. Statement Function definitions ..
286 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
287 * ..
288 * .. Executable Statements ..
289 *
290 * Decode the input arguments
291 *
292 IF( LSAME( JOBVL, 'N' ) ) THEN
293 IJOBVL = 1
294 ILVL = .FALSE.
295 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
296 IJOBVL = 2
297 ILVL = .TRUE.
298 ELSE
299 IJOBVL = -1
300 ILVL = .FALSE.
301 END IF
302 *
303 IF( LSAME( JOBVR, 'N' ) ) THEN
304 IJOBVR = 1
305 ILVR = .FALSE.
306 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
307 IJOBVR = 2
308 ILVR = .TRUE.
309 ELSE
310 IJOBVR = -1
311 ILVR = .FALSE.
312 END IF
313 ILV = ILVL .OR. ILVR
314 *
315 NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
316 WANTSN = LSAME( SENSE, 'N' )
317 WANTSE = LSAME( SENSE, 'E' )
318 WANTSV = LSAME( SENSE, 'V' )
319 WANTSB = LSAME( SENSE, 'B' )
320 *
321 * Test the input arguments
322 *
323 INFO = 0
324 LQUERY = ( LWORK.EQ.-1 )
325 IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
326 $ LSAME( BALANC, 'B' ) ) ) THEN
327 INFO = -1
328 ELSE IF( IJOBVL.LE.0 ) THEN
329 INFO = -2
330 ELSE IF( IJOBVR.LE.0 ) THEN
331 INFO = -3
332 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
333 $ THEN
334 INFO = -4
335 ELSE IF( N.LT.0 ) THEN
336 INFO = -5
337 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338 INFO = -7
339 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340 INFO = -9
341 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
342 INFO = -13
343 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
344 INFO = -15
345 END IF
346 *
347 * Compute workspace
348 * (Note: Comments in the code beginning "Workspace:" describe the
349 * minimal amount of workspace needed at that point in the code,
350 * as well as the preferred amount for good performance.
351 * NB refers to the optimal block size for the immediately
352 * following subroutine, as returned by ILAENV. The workspace is
353 * computed assuming ILO = 1 and IHI = N, the worst case.)
354 *
355 IF( INFO.EQ.0 ) THEN
356 IF( N.EQ.0 ) THEN
357 MINWRK = 1
358 MAXWRK = 1
359 ELSE
360 MINWRK = 2*N
361 IF( WANTSE ) THEN
362 MINWRK = 4*N
363 ELSE IF( WANTSV .OR. WANTSB ) THEN
364 MINWRK = 2*N*( N + 1)
365 END IF
366 MAXWRK = MINWRK
367 MAXWRK = MAX( MAXWRK,
368 $ N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
369 MAXWRK = MAX( MAXWRK,
370 $ N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
371 IF( ILVL ) THEN
372 MAXWRK = MAX( MAXWRK, N +
373 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
374 END IF
375 END IF
376 WORK( 1 ) = MAXWRK
377 *
378 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
379 INFO = -25
380 END IF
381 END IF
382 *
383 IF( INFO.NE.0 ) THEN
384 CALL XERBLA( 'ZGGEVX', -INFO )
385 RETURN
386 ELSE IF( LQUERY ) THEN
387 RETURN
388 END IF
389 *
390 * Quick return if possible
391 *
392 IF( N.EQ.0 )
393 $ RETURN
394 *
395 * Get machine constants
396 *
397 EPS = DLAMCH( 'P' )
398 SMLNUM = DLAMCH( 'S' )
399 BIGNUM = ONE / SMLNUM
400 CALL DLABAD( SMLNUM, BIGNUM )
401 SMLNUM = SQRT( SMLNUM ) / EPS
402 BIGNUM = ONE / SMLNUM
403 *
404 * Scale A if max element outside range [SMLNUM,BIGNUM]
405 *
406 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
407 ILASCL = .FALSE.
408 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
409 ANRMTO = SMLNUM
410 ILASCL = .TRUE.
411 ELSE IF( ANRM.GT.BIGNUM ) THEN
412 ANRMTO = BIGNUM
413 ILASCL = .TRUE.
414 END IF
415 IF( ILASCL )
416 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
417 *
418 * Scale B if max element outside range [SMLNUM,BIGNUM]
419 *
420 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
421 ILBSCL = .FALSE.
422 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
423 BNRMTO = SMLNUM
424 ILBSCL = .TRUE.
425 ELSE IF( BNRM.GT.BIGNUM ) THEN
426 BNRMTO = BIGNUM
427 ILBSCL = .TRUE.
428 END IF
429 IF( ILBSCL )
430 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
431 *
432 * Permute and/or balance the matrix pair (A,B)
433 * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
434 *
435 CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
436 $ RWORK, IERR )
437 *
438 * Compute ABNRM and BBNRM
439 *
440 ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
441 IF( ILASCL ) THEN
442 RWORK( 1 ) = ABNRM
443 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
444 $ IERR )
445 ABNRM = RWORK( 1 )
446 END IF
447 *
448 BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
449 IF( ILBSCL ) THEN
450 RWORK( 1 ) = BBNRM
451 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
452 $ IERR )
453 BBNRM = RWORK( 1 )
454 END IF
455 *
456 * Reduce B to triangular form (QR decomposition of B)
457 * (Complex Workspace: need N, prefer N*NB )
458 *
459 IROWS = IHI + 1 - ILO
460 IF( ILV .OR. .NOT.WANTSN ) THEN
461 ICOLS = N + 1 - ILO
462 ELSE
463 ICOLS = IROWS
464 END IF
465 ITAU = 1
466 IWRK = ITAU + IROWS
467 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
468 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
469 *
470 * Apply the unitary transformation to A
471 * (Complex Workspace: need N, prefer N*NB)
472 *
473 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
474 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
475 $ LWORK+1-IWRK, IERR )
476 *
477 * Initialize VL and/or VR
478 * (Workspace: need N, prefer N*NB)
479 *
480 IF( ILVL ) THEN
481 CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
482 IF( IROWS.GT.1 ) THEN
483 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
484 $ VL( ILO+1, ILO ), LDVL )
485 END IF
486 CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
487 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
488 END IF
489 *
490 IF( ILVR )
491 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
492 *
493 * Reduce to generalized Hessenberg form
494 * (Workspace: none needed)
495 *
496 IF( ILV .OR. .NOT.WANTSN ) THEN
497 *
498 * Eigenvectors requested -- work on whole matrix.
499 *
500 CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
501 $ LDVL, VR, LDVR, IERR )
502 ELSE
503 CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
504 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
505 END IF
506 *
507 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
508 * Schur forms and Schur vectors)
509 * (Complex Workspace: need N)
510 * (Real Workspace: need N)
511 *
512 IWRK = ITAU
513 IF( ILV .OR. .NOT.WANTSN ) THEN
514 CHTEMP = 'S'
515 ELSE
516 CHTEMP = 'E'
517 END IF
518 *
519 CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
520 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
521 $ LWORK+1-IWRK, RWORK, IERR )
522 IF( IERR.NE.0 ) THEN
523 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
524 INFO = IERR
525 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
526 INFO = IERR - N
527 ELSE
528 INFO = N + 1
529 END IF
530 GO TO 90
531 END IF
532 *
533 * Compute Eigenvectors and estimate condition numbers if desired
534 * ZTGEVC: (Complex Workspace: need 2*N )
535 * (Real Workspace: need 2*N )
536 * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
537 * (Integer Workspace: need N+2 )
538 *
539 IF( ILV .OR. .NOT.WANTSN ) THEN
540 IF( ILV ) THEN
541 IF( ILVL ) THEN
542 IF( ILVR ) THEN
543 CHTEMP = 'B'
544 ELSE
545 CHTEMP = 'L'
546 END IF
547 ELSE
548 CHTEMP = 'R'
549 END IF
550 *
551 CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
552 $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
553 $ IERR )
554 IF( IERR.NE.0 ) THEN
555 INFO = N + 2
556 GO TO 90
557 END IF
558 END IF
559 *
560 IF( .NOT.WANTSN ) THEN
561 *
562 * compute eigenvectors (DTGEVC) and estimate condition
563 * numbers (DTGSNA). Note that the definition of the condition
564 * number is not invariant under transformation (u,v) to
565 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
566 * Schur form (S,T), Q and Z are orthogonal matrices. In order
567 * to avoid using extra 2*N*N workspace, we have to
568 * re-calculate eigenvectors and estimate the condition numbers
569 * one at a time.
570 *
571 DO 20 I = 1, N
572 *
573 DO 10 J = 1, N
574 BWORK( J ) = .FALSE.
575 10 CONTINUE
576 BWORK( I ) = .TRUE.
577 *
578 IWRK = N + 1
579 IWRK1 = IWRK + N
580 *
581 IF( WANTSE .OR. WANTSB ) THEN
582 CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
583 $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
584 $ WORK( IWRK1 ), RWORK, IERR )
585 IF( IERR.NE.0 ) THEN
586 INFO = N + 2
587 GO TO 90
588 END IF
589 END IF
590 *
591 CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
592 $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
593 $ RCONDV( I ), 1, M, WORK( IWRK1 ),
594 $ LWORK-IWRK1+1, IWORK, IERR )
595 *
596 20 CONTINUE
597 END IF
598 END IF
599 *
600 * Undo balancing on VL and VR and normalization
601 * (Workspace: none needed)
602 *
603 IF( ILVL ) THEN
604 CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
605 $ LDVL, IERR )
606 *
607 DO 50 JC = 1, N
608 TEMP = ZERO
609 DO 30 JR = 1, N
610 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
611 30 CONTINUE
612 IF( TEMP.LT.SMLNUM )
613 $ GO TO 50
614 TEMP = ONE / TEMP
615 DO 40 JR = 1, N
616 VL( JR, JC ) = VL( JR, JC )*TEMP
617 40 CONTINUE
618 50 CONTINUE
619 END IF
620 *
621 IF( ILVR ) THEN
622 CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
623 $ LDVR, IERR )
624 DO 80 JC = 1, N
625 TEMP = ZERO
626 DO 60 JR = 1, N
627 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
628 60 CONTINUE
629 IF( TEMP.LT.SMLNUM )
630 $ GO TO 80
631 TEMP = ONE / TEMP
632 DO 70 JR = 1, N
633 VR( JR, JC ) = VR( JR, JC )*TEMP
634 70 CONTINUE
635 80 CONTINUE
636 END IF
637 *
638 * Undo scaling if necessary
639 *
640 IF( ILASCL )
641 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
642 *
643 IF( ILBSCL )
644 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
645 *
646 90 CONTINUE
647 WORK( 1 ) = MAXWRK
648 *
649 RETURN
650 *
651 * End of ZGGEVX
652 *
653 END