1 SUBROUTINE ZGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
2 $ SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
3 $ LWORK, RWORK, BWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBVSL, JOBVSR, SORT
12 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 DOUBLE PRECISION RWORK( * )
17 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
18 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
19 $ WORK( * )
20 * ..
21 * .. Function Arguments ..
22 LOGICAL SELCTG
23 EXTERNAL SELCTG
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
30 * (A,B), the generalized eigenvalues, the generalized complex Schur
31 * form (S, T), and optionally left and/or right Schur vectors (VSL
32 * and VSR). This gives the generalized Schur factorization
33 *
34 * (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
35 *
36 * where (VSR)**H is the conjugate-transpose of VSR.
37 *
38 * Optionally, it also orders the eigenvalues so that a selected cluster
39 * of eigenvalues appears in the leading diagonal blocks of the upper
40 * triangular matrix S and the upper triangular matrix T. The leading
41 * columns of VSL and VSR then form an unitary basis for the
42 * corresponding left and right eigenspaces (deflating subspaces).
43 *
44 * (If only the generalized eigenvalues are needed, use the driver
45 * ZGGEV instead, which is faster.)
46 *
47 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
48 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
49 * usually represented as the pair (alpha,beta), as there is a
50 * reasonable interpretation for beta=0, and even for both being zero.
51 *
52 * A pair of matrices (S,T) is in generalized complex Schur form if S
53 * and T are upper triangular and, in addition, the diagonal elements
54 * of T are non-negative real numbers.
55 *
56 * Arguments
57 * =========
58 *
59 * JOBVSL (input) CHARACTER*1
60 * = 'N': do not compute the left Schur vectors;
61 * = 'V': compute the left Schur vectors.
62 *
63 * JOBVSR (input) CHARACTER*1
64 * = 'N': do not compute the right Schur vectors;
65 * = 'V': compute the right Schur vectors.
66 *
67 * SORT (input) CHARACTER*1
68 * Specifies whether or not to order the eigenvalues on the
69 * diagonal of the generalized Schur form.
70 * = 'N': Eigenvalues are not ordered;
71 * = 'S': Eigenvalues are ordered (see SELCTG).
72 *
73 * SELCTG (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
74 * SELCTG must be declared EXTERNAL in the calling subroutine.
75 * If SORT = 'N', SELCTG is not referenced.
76 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
77 * to the top left of the Schur form.
78 * An eigenvalue ALPHA(j)/BETA(j) is selected if
79 * SELCTG(ALPHA(j),BETA(j)) is true.
80 *
81 * Note that a selected complex eigenvalue may no longer satisfy
82 * SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
83 * ordering may change the value of complex eigenvalues
84 * (especially if the eigenvalue is ill-conditioned), in this
85 * case INFO is set to N+2 (See INFO below).
86 *
87 * N (input) INTEGER
88 * The order of the matrices A, B, VSL, and VSR. N >= 0.
89 *
90 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
91 * On entry, the first of the pair of matrices.
92 * On exit, A has been overwritten by its generalized Schur
93 * form S.
94 *
95 * LDA (input) INTEGER
96 * The leading dimension of A. LDA >= max(1,N).
97 *
98 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
99 * On entry, the second of the pair of matrices.
100 * On exit, B has been overwritten by its generalized Schur
101 * form T.
102 *
103 * LDB (input) INTEGER
104 * The leading dimension of B. LDB >= max(1,N).
105 *
106 * SDIM (output) INTEGER
107 * If SORT = 'N', SDIM = 0.
108 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
109 * for which SELCTG is true.
110 *
111 * ALPHA (output) COMPLEX*16 array, dimension (N)
112 * BETA (output) COMPLEX*16 array, dimension (N)
113 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
114 * generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
115 * j=1,...,N are the diagonals of the complex Schur form (A,B)
116 * output by ZGGES. The BETA(j) will be non-negative real.
117 *
118 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
119 * underflow, and BETA(j) may even be zero. Thus, the user
120 * should avoid naively computing the ratio alpha/beta.
121 * However, ALPHA will be always less than and usually
122 * comparable with norm(A) in magnitude, and BETA always less
123 * than and usually comparable with norm(B).
124 *
125 * VSL (output) COMPLEX*16 array, dimension (LDVSL,N)
126 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
127 * Not referenced if JOBVSL = 'N'.
128 *
129 * LDVSL (input) INTEGER
130 * The leading dimension of the matrix VSL. LDVSL >= 1, and
131 * if JOBVSL = 'V', LDVSL >= N.
132 *
133 * VSR (output) COMPLEX*16 array, dimension (LDVSR,N)
134 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
135 * Not referenced if JOBVSR = 'N'.
136 *
137 * LDVSR (input) INTEGER
138 * The leading dimension of the matrix VSR. LDVSR >= 1, and
139 * if JOBVSR = 'V', LDVSR >= N.
140 *
141 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
142 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143 *
144 * LWORK (input) INTEGER
145 * The dimension of the array WORK. LWORK >= max(1,2*N).
146 * For good performance, LWORK must generally be larger.
147 *
148 * If LWORK = -1, then a workspace query is assumed; the routine
149 * only calculates the optimal size of the WORK array, returns
150 * this value as the first entry of the WORK array, and no error
151 * message related to LWORK is issued by XERBLA.
152 *
153 * RWORK (workspace) DOUBLE PRECISION array, dimension (8*N)
154 *
155 * BWORK (workspace) LOGICAL array, dimension (N)
156 * Not referenced if SORT = 'N'.
157 *
158 * INFO (output) INTEGER
159 * = 0: successful exit
160 * < 0: if INFO = -i, the i-th argument had an illegal value.
161 * =1,...,N:
162 * The QZ iteration failed. (A,B) are not in Schur
163 * form, but ALPHA(j) and BETA(j) should be correct for
164 * j=INFO+1,...,N.
165 * > N: =N+1: other than QZ iteration failed in ZHGEQZ
166 * =N+2: after reordering, roundoff changed values of
167 * some complex eigenvalues so that leading
168 * eigenvalues in the Generalized Schur form no
169 * longer satisfy SELCTG=.TRUE. This could also
170 * be caused due to scaling.
171 * =N+3: reordering falied in ZTGSEN.
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176 DOUBLE PRECISION ZERO, ONE
177 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
178 COMPLEX*16 CZERO, CONE
179 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
180 $ CONE = ( 1.0D0, 0.0D0 ) )
181 * ..
182 * .. Local Scalars ..
183 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
184 $ LQUERY, WANTST
185 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
186 $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
187 $ LWKOPT
188 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
189 $ PVSR, SMLNUM
190 * ..
191 * .. Local Arrays ..
192 INTEGER IDUM( 1 )
193 DOUBLE PRECISION DIF( 2 )
194 * ..
195 * .. External Subroutines ..
196 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
197 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
198 $ ZUNMQR
199 * ..
200 * .. External Functions ..
201 LOGICAL LSAME
202 INTEGER ILAENV
203 DOUBLE PRECISION DLAMCH, ZLANGE
204 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
205 * ..
206 * .. Intrinsic Functions ..
207 INTRINSIC MAX, SQRT
208 * ..
209 * .. Executable Statements ..
210 *
211 * Decode the input arguments
212 *
213 IF( LSAME( JOBVSL, 'N' ) ) THEN
214 IJOBVL = 1
215 ILVSL = .FALSE.
216 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
217 IJOBVL = 2
218 ILVSL = .TRUE.
219 ELSE
220 IJOBVL = -1
221 ILVSL = .FALSE.
222 END IF
223 *
224 IF( LSAME( JOBVSR, 'N' ) ) THEN
225 IJOBVR = 1
226 ILVSR = .FALSE.
227 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
228 IJOBVR = 2
229 ILVSR = .TRUE.
230 ELSE
231 IJOBVR = -1
232 ILVSR = .FALSE.
233 END IF
234 *
235 WANTST = LSAME( SORT, 'S' )
236 *
237 * Test the input arguments
238 *
239 INFO = 0
240 LQUERY = ( LWORK.EQ.-1 )
241 IF( IJOBVL.LE.0 ) THEN
242 INFO = -1
243 ELSE IF( IJOBVR.LE.0 ) THEN
244 INFO = -2
245 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
246 INFO = -3
247 ELSE IF( N.LT.0 ) THEN
248 INFO = -5
249 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
250 INFO = -7
251 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
252 INFO = -9
253 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
254 INFO = -14
255 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
256 INFO = -16
257 END IF
258 *
259 * Compute workspace
260 * (Note: Comments in the code beginning "Workspace:" describe the
261 * minimal amount of workspace needed at that point in the code,
262 * as well as the preferred amount for good performance.
263 * NB refers to the optimal block size for the immediately
264 * following subroutine, as returned by ILAENV.)
265 *
266 IF( INFO.EQ.0 ) THEN
267 LWKMIN = MAX( 1, 2*N )
268 LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
269 LWKOPT = MAX( LWKOPT, N +
270 $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) )
271 IF( ILVSL ) THEN
272 LWKOPT = MAX( LWKOPT, N +
273 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
274 END IF
275 WORK( 1 ) = LWKOPT
276 *
277 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
278 $ INFO = -18
279 END IF
280 *
281 IF( INFO.NE.0 ) THEN
282 CALL XERBLA( 'ZGGES ', -INFO )
283 RETURN
284 ELSE IF( LQUERY ) THEN
285 RETURN
286 END IF
287 *
288 * Quick return if possible
289 *
290 IF( N.EQ.0 ) THEN
291 SDIM = 0
292 RETURN
293 END IF
294 *
295 * Get machine constants
296 *
297 EPS = DLAMCH( 'P' )
298 SMLNUM = DLAMCH( 'S' )
299 BIGNUM = ONE / SMLNUM
300 CALL DLABAD( SMLNUM, BIGNUM )
301 SMLNUM = SQRT( SMLNUM ) / EPS
302 BIGNUM = ONE / SMLNUM
303 *
304 * Scale A if max element outside range [SMLNUM,BIGNUM]
305 *
306 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
307 ILASCL = .FALSE.
308 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
309 ANRMTO = SMLNUM
310 ILASCL = .TRUE.
311 ELSE IF( ANRM.GT.BIGNUM ) THEN
312 ANRMTO = BIGNUM
313 ILASCL = .TRUE.
314 END IF
315 *
316 IF( ILASCL )
317 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
318 *
319 * Scale B if max element outside range [SMLNUM,BIGNUM]
320 *
321 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
322 ILBSCL = .FALSE.
323 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
324 BNRMTO = SMLNUM
325 ILBSCL = .TRUE.
326 ELSE IF( BNRM.GT.BIGNUM ) THEN
327 BNRMTO = BIGNUM
328 ILBSCL = .TRUE.
329 END IF
330 *
331 IF( ILBSCL )
332 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
333 *
334 * Permute the matrix to make it more nearly triangular
335 * (Real Workspace: need 6*N)
336 *
337 ILEFT = 1
338 IRIGHT = N + 1
339 IRWRK = IRIGHT + N
340 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
341 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
342 *
343 * Reduce B to triangular form (QR decomposition of B)
344 * (Complex Workspace: need N, prefer N*NB)
345 *
346 IROWS = IHI + 1 - ILO
347 ICOLS = N + 1 - ILO
348 ITAU = 1
349 IWRK = ITAU + IROWS
350 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
351 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
352 *
353 * Apply the orthogonal transformation to matrix A
354 * (Complex Workspace: need N, prefer N*NB)
355 *
356 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
357 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
358 $ LWORK+1-IWRK, IERR )
359 *
360 * Initialize VSL
361 * (Complex Workspace: need N, prefer N*NB)
362 *
363 IF( ILVSL ) THEN
364 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
365 IF( IROWS.GT.1 ) THEN
366 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
367 $ VSL( ILO+1, ILO ), LDVSL )
368 END IF
369 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
370 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
371 END IF
372 *
373 * Initialize VSR
374 *
375 IF( ILVSR )
376 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
377 *
378 * Reduce to generalized Hessenberg form
379 * (Workspace: none needed)
380 *
381 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
382 $ LDVSL, VSR, LDVSR, IERR )
383 *
384 SDIM = 0
385 *
386 * Perform QZ algorithm, computing Schur vectors if desired
387 * (Complex Workspace: need N)
388 * (Real Workspace: need N)
389 *
390 IWRK = ITAU
391 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
392 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
393 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
394 IF( IERR.NE.0 ) THEN
395 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
396 INFO = IERR
397 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
398 INFO = IERR - N
399 ELSE
400 INFO = N + 1
401 END IF
402 GO TO 30
403 END IF
404 *
405 * Sort eigenvalues ALPHA/BETA if desired
406 * (Workspace: none needed)
407 *
408 IF( WANTST ) THEN
409 *
410 * Undo scaling on eigenvalues before selecting
411 *
412 IF( ILASCL )
413 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
414 IF( ILBSCL )
415 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
416 *
417 * Select eigenvalues
418 *
419 DO 10 I = 1, N
420 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
421 10 CONTINUE
422 *
423 CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
424 $ BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
425 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
426 IF( IERR.EQ.1 )
427 $ INFO = N + 3
428 *
429 END IF
430 *
431 * Apply back-permutation to VSL and VSR
432 * (Workspace: none needed)
433 *
434 IF( ILVSL )
435 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
436 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
437 IF( ILVSR )
438 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
439 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
440 *
441 * Undo scaling
442 *
443 IF( ILASCL ) THEN
444 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
445 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
446 END IF
447 *
448 IF( ILBSCL ) THEN
449 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
450 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
451 END IF
452 *
453 IF( WANTST ) THEN
454 *
455 * Check if reordering is correct
456 *
457 LASTSL = .TRUE.
458 SDIM = 0
459 DO 20 I = 1, N
460 CURSL = SELCTG( ALPHA( I ), BETA( I ) )
461 IF( CURSL )
462 $ SDIM = SDIM + 1
463 IF( CURSL .AND. .NOT.LASTSL )
464 $ INFO = N + 2
465 LASTSL = CURSL
466 20 CONTINUE
467 *
468 END IF
469 *
470 30 CONTINUE
471 *
472 WORK( 1 ) = LWKOPT
473 *
474 RETURN
475 *
476 * End of ZGGES
477 *
478 END
2 $ SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
3 $ LWORK, RWORK, BWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBVSL, JOBVSR, SORT
12 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 DOUBLE PRECISION RWORK( * )
17 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
18 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
19 $ WORK( * )
20 * ..
21 * .. Function Arguments ..
22 LOGICAL SELCTG
23 EXTERNAL SELCTG
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
30 * (A,B), the generalized eigenvalues, the generalized complex Schur
31 * form (S, T), and optionally left and/or right Schur vectors (VSL
32 * and VSR). This gives the generalized Schur factorization
33 *
34 * (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
35 *
36 * where (VSR)**H is the conjugate-transpose of VSR.
37 *
38 * Optionally, it also orders the eigenvalues so that a selected cluster
39 * of eigenvalues appears in the leading diagonal blocks of the upper
40 * triangular matrix S and the upper triangular matrix T. The leading
41 * columns of VSL and VSR then form an unitary basis for the
42 * corresponding left and right eigenspaces (deflating subspaces).
43 *
44 * (If only the generalized eigenvalues are needed, use the driver
45 * ZGGEV instead, which is faster.)
46 *
47 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
48 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
49 * usually represented as the pair (alpha,beta), as there is a
50 * reasonable interpretation for beta=0, and even for both being zero.
51 *
52 * A pair of matrices (S,T) is in generalized complex Schur form if S
53 * and T are upper triangular and, in addition, the diagonal elements
54 * of T are non-negative real numbers.
55 *
56 * Arguments
57 * =========
58 *
59 * JOBVSL (input) CHARACTER*1
60 * = 'N': do not compute the left Schur vectors;
61 * = 'V': compute the left Schur vectors.
62 *
63 * JOBVSR (input) CHARACTER*1
64 * = 'N': do not compute the right Schur vectors;
65 * = 'V': compute the right Schur vectors.
66 *
67 * SORT (input) CHARACTER*1
68 * Specifies whether or not to order the eigenvalues on the
69 * diagonal of the generalized Schur form.
70 * = 'N': Eigenvalues are not ordered;
71 * = 'S': Eigenvalues are ordered (see SELCTG).
72 *
73 * SELCTG (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
74 * SELCTG must be declared EXTERNAL in the calling subroutine.
75 * If SORT = 'N', SELCTG is not referenced.
76 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
77 * to the top left of the Schur form.
78 * An eigenvalue ALPHA(j)/BETA(j) is selected if
79 * SELCTG(ALPHA(j),BETA(j)) is true.
80 *
81 * Note that a selected complex eigenvalue may no longer satisfy
82 * SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
83 * ordering may change the value of complex eigenvalues
84 * (especially if the eigenvalue is ill-conditioned), in this
85 * case INFO is set to N+2 (See INFO below).
86 *
87 * N (input) INTEGER
88 * The order of the matrices A, B, VSL, and VSR. N >= 0.
89 *
90 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
91 * On entry, the first of the pair of matrices.
92 * On exit, A has been overwritten by its generalized Schur
93 * form S.
94 *
95 * LDA (input) INTEGER
96 * The leading dimension of A. LDA >= max(1,N).
97 *
98 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
99 * On entry, the second of the pair of matrices.
100 * On exit, B has been overwritten by its generalized Schur
101 * form T.
102 *
103 * LDB (input) INTEGER
104 * The leading dimension of B. LDB >= max(1,N).
105 *
106 * SDIM (output) INTEGER
107 * If SORT = 'N', SDIM = 0.
108 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
109 * for which SELCTG is true.
110 *
111 * ALPHA (output) COMPLEX*16 array, dimension (N)
112 * BETA (output) COMPLEX*16 array, dimension (N)
113 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
114 * generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
115 * j=1,...,N are the diagonals of the complex Schur form (A,B)
116 * output by ZGGES. The BETA(j) will be non-negative real.
117 *
118 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
119 * underflow, and BETA(j) may even be zero. Thus, the user
120 * should avoid naively computing the ratio alpha/beta.
121 * However, ALPHA will be always less than and usually
122 * comparable with norm(A) in magnitude, and BETA always less
123 * than and usually comparable with norm(B).
124 *
125 * VSL (output) COMPLEX*16 array, dimension (LDVSL,N)
126 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
127 * Not referenced if JOBVSL = 'N'.
128 *
129 * LDVSL (input) INTEGER
130 * The leading dimension of the matrix VSL. LDVSL >= 1, and
131 * if JOBVSL = 'V', LDVSL >= N.
132 *
133 * VSR (output) COMPLEX*16 array, dimension (LDVSR,N)
134 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
135 * Not referenced if JOBVSR = 'N'.
136 *
137 * LDVSR (input) INTEGER
138 * The leading dimension of the matrix VSR. LDVSR >= 1, and
139 * if JOBVSR = 'V', LDVSR >= N.
140 *
141 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
142 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143 *
144 * LWORK (input) INTEGER
145 * The dimension of the array WORK. LWORK >= max(1,2*N).
146 * For good performance, LWORK must generally be larger.
147 *
148 * If LWORK = -1, then a workspace query is assumed; the routine
149 * only calculates the optimal size of the WORK array, returns
150 * this value as the first entry of the WORK array, and no error
151 * message related to LWORK is issued by XERBLA.
152 *
153 * RWORK (workspace) DOUBLE PRECISION array, dimension (8*N)
154 *
155 * BWORK (workspace) LOGICAL array, dimension (N)
156 * Not referenced if SORT = 'N'.
157 *
158 * INFO (output) INTEGER
159 * = 0: successful exit
160 * < 0: if INFO = -i, the i-th argument had an illegal value.
161 * =1,...,N:
162 * The QZ iteration failed. (A,B) are not in Schur
163 * form, but ALPHA(j) and BETA(j) should be correct for
164 * j=INFO+1,...,N.
165 * > N: =N+1: other than QZ iteration failed in ZHGEQZ
166 * =N+2: after reordering, roundoff changed values of
167 * some complex eigenvalues so that leading
168 * eigenvalues in the Generalized Schur form no
169 * longer satisfy SELCTG=.TRUE. This could also
170 * be caused due to scaling.
171 * =N+3: reordering falied in ZTGSEN.
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176 DOUBLE PRECISION ZERO, ONE
177 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
178 COMPLEX*16 CZERO, CONE
179 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
180 $ CONE = ( 1.0D0, 0.0D0 ) )
181 * ..
182 * .. Local Scalars ..
183 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
184 $ LQUERY, WANTST
185 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
186 $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
187 $ LWKOPT
188 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
189 $ PVSR, SMLNUM
190 * ..
191 * .. Local Arrays ..
192 INTEGER IDUM( 1 )
193 DOUBLE PRECISION DIF( 2 )
194 * ..
195 * .. External Subroutines ..
196 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
197 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
198 $ ZUNMQR
199 * ..
200 * .. External Functions ..
201 LOGICAL LSAME
202 INTEGER ILAENV
203 DOUBLE PRECISION DLAMCH, ZLANGE
204 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
205 * ..
206 * .. Intrinsic Functions ..
207 INTRINSIC MAX, SQRT
208 * ..
209 * .. Executable Statements ..
210 *
211 * Decode the input arguments
212 *
213 IF( LSAME( JOBVSL, 'N' ) ) THEN
214 IJOBVL = 1
215 ILVSL = .FALSE.
216 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
217 IJOBVL = 2
218 ILVSL = .TRUE.
219 ELSE
220 IJOBVL = -1
221 ILVSL = .FALSE.
222 END IF
223 *
224 IF( LSAME( JOBVSR, 'N' ) ) THEN
225 IJOBVR = 1
226 ILVSR = .FALSE.
227 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
228 IJOBVR = 2
229 ILVSR = .TRUE.
230 ELSE
231 IJOBVR = -1
232 ILVSR = .FALSE.
233 END IF
234 *
235 WANTST = LSAME( SORT, 'S' )
236 *
237 * Test the input arguments
238 *
239 INFO = 0
240 LQUERY = ( LWORK.EQ.-1 )
241 IF( IJOBVL.LE.0 ) THEN
242 INFO = -1
243 ELSE IF( IJOBVR.LE.0 ) THEN
244 INFO = -2
245 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
246 INFO = -3
247 ELSE IF( N.LT.0 ) THEN
248 INFO = -5
249 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
250 INFO = -7
251 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
252 INFO = -9
253 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
254 INFO = -14
255 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
256 INFO = -16
257 END IF
258 *
259 * Compute workspace
260 * (Note: Comments in the code beginning "Workspace:" describe the
261 * minimal amount of workspace needed at that point in the code,
262 * as well as the preferred amount for good performance.
263 * NB refers to the optimal block size for the immediately
264 * following subroutine, as returned by ILAENV.)
265 *
266 IF( INFO.EQ.0 ) THEN
267 LWKMIN = MAX( 1, 2*N )
268 LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
269 LWKOPT = MAX( LWKOPT, N +
270 $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) )
271 IF( ILVSL ) THEN
272 LWKOPT = MAX( LWKOPT, N +
273 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
274 END IF
275 WORK( 1 ) = LWKOPT
276 *
277 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
278 $ INFO = -18
279 END IF
280 *
281 IF( INFO.NE.0 ) THEN
282 CALL XERBLA( 'ZGGES ', -INFO )
283 RETURN
284 ELSE IF( LQUERY ) THEN
285 RETURN
286 END IF
287 *
288 * Quick return if possible
289 *
290 IF( N.EQ.0 ) THEN
291 SDIM = 0
292 RETURN
293 END IF
294 *
295 * Get machine constants
296 *
297 EPS = DLAMCH( 'P' )
298 SMLNUM = DLAMCH( 'S' )
299 BIGNUM = ONE / SMLNUM
300 CALL DLABAD( SMLNUM, BIGNUM )
301 SMLNUM = SQRT( SMLNUM ) / EPS
302 BIGNUM = ONE / SMLNUM
303 *
304 * Scale A if max element outside range [SMLNUM,BIGNUM]
305 *
306 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
307 ILASCL = .FALSE.
308 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
309 ANRMTO = SMLNUM
310 ILASCL = .TRUE.
311 ELSE IF( ANRM.GT.BIGNUM ) THEN
312 ANRMTO = BIGNUM
313 ILASCL = .TRUE.
314 END IF
315 *
316 IF( ILASCL )
317 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
318 *
319 * Scale B if max element outside range [SMLNUM,BIGNUM]
320 *
321 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
322 ILBSCL = .FALSE.
323 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
324 BNRMTO = SMLNUM
325 ILBSCL = .TRUE.
326 ELSE IF( BNRM.GT.BIGNUM ) THEN
327 BNRMTO = BIGNUM
328 ILBSCL = .TRUE.
329 END IF
330 *
331 IF( ILBSCL )
332 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
333 *
334 * Permute the matrix to make it more nearly triangular
335 * (Real Workspace: need 6*N)
336 *
337 ILEFT = 1
338 IRIGHT = N + 1
339 IRWRK = IRIGHT + N
340 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
341 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
342 *
343 * Reduce B to triangular form (QR decomposition of B)
344 * (Complex Workspace: need N, prefer N*NB)
345 *
346 IROWS = IHI + 1 - ILO
347 ICOLS = N + 1 - ILO
348 ITAU = 1
349 IWRK = ITAU + IROWS
350 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
351 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
352 *
353 * Apply the orthogonal transformation to matrix A
354 * (Complex Workspace: need N, prefer N*NB)
355 *
356 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
357 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
358 $ LWORK+1-IWRK, IERR )
359 *
360 * Initialize VSL
361 * (Complex Workspace: need N, prefer N*NB)
362 *
363 IF( ILVSL ) THEN
364 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
365 IF( IROWS.GT.1 ) THEN
366 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
367 $ VSL( ILO+1, ILO ), LDVSL )
368 END IF
369 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
370 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
371 END IF
372 *
373 * Initialize VSR
374 *
375 IF( ILVSR )
376 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
377 *
378 * Reduce to generalized Hessenberg form
379 * (Workspace: none needed)
380 *
381 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
382 $ LDVSL, VSR, LDVSR, IERR )
383 *
384 SDIM = 0
385 *
386 * Perform QZ algorithm, computing Schur vectors if desired
387 * (Complex Workspace: need N)
388 * (Real Workspace: need N)
389 *
390 IWRK = ITAU
391 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
392 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
393 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
394 IF( IERR.NE.0 ) THEN
395 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
396 INFO = IERR
397 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
398 INFO = IERR - N
399 ELSE
400 INFO = N + 1
401 END IF
402 GO TO 30
403 END IF
404 *
405 * Sort eigenvalues ALPHA/BETA if desired
406 * (Workspace: none needed)
407 *
408 IF( WANTST ) THEN
409 *
410 * Undo scaling on eigenvalues before selecting
411 *
412 IF( ILASCL )
413 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
414 IF( ILBSCL )
415 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
416 *
417 * Select eigenvalues
418 *
419 DO 10 I = 1, N
420 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
421 10 CONTINUE
422 *
423 CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
424 $ BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
425 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
426 IF( IERR.EQ.1 )
427 $ INFO = N + 3
428 *
429 END IF
430 *
431 * Apply back-permutation to VSL and VSR
432 * (Workspace: none needed)
433 *
434 IF( ILVSL )
435 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
436 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
437 IF( ILVSR )
438 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
439 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
440 *
441 * Undo scaling
442 *
443 IF( ILASCL ) THEN
444 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
445 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
446 END IF
447 *
448 IF( ILBSCL ) THEN
449 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
450 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
451 END IF
452 *
453 IF( WANTST ) THEN
454 *
455 * Check if reordering is correct
456 *
457 LASTSL = .TRUE.
458 SDIM = 0
459 DO 20 I = 1, N
460 CURSL = SELCTG( ALPHA( I ), BETA( I ) )
461 IF( CURSL )
462 $ SDIM = SDIM + 1
463 IF( CURSL .AND. .NOT.LASTSL )
464 $ INFO = N + 2
465 LASTSL = CURSL
466 20 CONTINUE
467 *
468 END IF
469 *
470 30 CONTINUE
471 *
472 WORK( 1 ) = LWKOPT
473 *
474 RETURN
475 *
476 * End of ZGGES
477 *
478 END