1 SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
2 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
3 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
4 $ LIWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2.1) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * -- April 2009 --
10 *
11 * .. Scalar Arguments ..
12 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
13 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
14 $ SDIM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
20 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
21 $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
22 $ WORK( * )
23 * ..
24 * .. Function Arguments ..
25 LOGICAL SELCTG
26 EXTERNAL SELCTG
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * DGGESX computes for a pair of N-by-N real nonsymmetric matrices
33 * (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
34 * optionally, the left and/or right matrices of Schur vectors (VSL and
35 * VSR). This gives the generalized Schur factorization
36 *
37 * (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
38 *
39 * Optionally, it also orders the eigenvalues so that a selected cluster
40 * of eigenvalues appears in the leading diagonal blocks of the upper
41 * quasi-triangular matrix S and the upper triangular matrix T; computes
42 * a reciprocal condition number for the average of the selected
43 * eigenvalues (RCONDE); and computes a reciprocal condition number for
44 * the right and left deflating subspaces corresponding to the selected
45 * eigenvalues (RCONDV). The leading columns of VSL and VSR then form
46 * an orthonormal basis for the corresponding left and right eigenspaces
47 * (deflating subspaces).
48 *
49 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
50 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
51 * usually represented as the pair (alpha,beta), as there is a
52 * reasonable interpretation for beta=0 or for both being zero.
53 *
54 * A pair of matrices (S,T) is in generalized real Schur form if T is
55 * upper triangular with non-negative diagonal and S is block upper
56 * triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
57 * to real generalized eigenvalues, while 2-by-2 blocks of S will be
58 * "standardized" by making the corresponding elements of T have the
59 * form:
60 * [ a 0 ]
61 * [ 0 b ]
62 *
63 * and the pair of corresponding 2-by-2 blocks in S and T will have a
64 * complex conjugate pair of generalized eigenvalues.
65 *
66 *
67 * Arguments
68 * =========
69 *
70 * JOBVSL (input) CHARACTER*1
71 * = 'N': do not compute the left Schur vectors;
72 * = 'V': compute the left Schur vectors.
73 *
74 * JOBVSR (input) CHARACTER*1
75 * = 'N': do not compute the right Schur vectors;
76 * = 'V': compute the right Schur vectors.
77 *
78 * SORT (input) CHARACTER*1
79 * Specifies whether or not to order the eigenvalues on the
80 * diagonal of the generalized Schur form.
81 * = 'N': Eigenvalues are not ordered;
82 * = 'S': Eigenvalues are ordered (see SELCTG).
83 *
84 * SELCTG (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
85 * SELCTG must be declared EXTERNAL in the calling subroutine.
86 * If SORT = 'N', SELCTG is not referenced.
87 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
88 * to the top left of the Schur form.
89 * An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
90 * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
91 * one of a complex conjugate pair of eigenvalues is selected,
92 * then both complex eigenvalues are selected.
93 * Note that a selected complex eigenvalue may no longer satisfy
94 * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
95 * since ordering may change the value of complex eigenvalues
96 * (especially if the eigenvalue is ill-conditioned), in this
97 * case INFO is set to N+3.
98 *
99 * SENSE (input) CHARACTER*1
100 * Determines which reciprocal condition numbers are computed.
101 * = 'N' : None are computed;
102 * = 'E' : Computed for average of selected eigenvalues only;
103 * = 'V' : Computed for selected deflating subspaces only;
104 * = 'B' : Computed for both.
105 * If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
106 *
107 * N (input) INTEGER
108 * The order of the matrices A, B, VSL, and VSR. N >= 0.
109 *
110 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
111 * On entry, the first of the pair of matrices.
112 * On exit, A has been overwritten by its generalized Schur
113 * form S.
114 *
115 * LDA (input) INTEGER
116 * The leading dimension of A. LDA >= max(1,N).
117 *
118 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
119 * On entry, the second of the pair of matrices.
120 * On exit, B has been overwritten by its generalized Schur
121 * form T.
122 *
123 * LDB (input) INTEGER
124 * The leading dimension of B. LDB >= max(1,N).
125 *
126 * SDIM (output) INTEGER
127 * If SORT = 'N', SDIM = 0.
128 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
129 * for which SELCTG is true. (Complex conjugate pairs for which
130 * SELCTG is true for either eigenvalue count as 2.)
131 *
132 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
133 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
134 * BETA (output) DOUBLE PRECISION array, dimension (N)
135 * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
136 * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
137 * and BETA(j),j=1,...,N are the diagonals of the complex Schur
138 * form (S,T) that would result if the 2-by-2 diagonal blocks of
139 * the real Schur form of (A,B) were further reduced to
140 * triangular form using 2-by-2 complex unitary transformations.
141 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
142 * positive, then the j-th and (j+1)-st eigenvalues are a
143 * complex conjugate pair, with ALPHAI(j+1) negative.
144 *
145 * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
146 * may easily over- or underflow, and BETA(j) may even be zero.
147 * Thus, the user should avoid naively computing the ratio.
148 * However, ALPHAR and ALPHAI will be always less than and
149 * usually comparable with norm(A) in magnitude, and BETA always
150 * less than and usually comparable with norm(B).
151 *
152 * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
153 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
154 * Not referenced if JOBVSL = 'N'.
155 *
156 * LDVSL (input) INTEGER
157 * The leading dimension of the matrix VSL. LDVSL >=1, and
158 * if JOBVSL = 'V', LDVSL >= N.
159 *
160 * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
161 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
162 * Not referenced if JOBVSR = 'N'.
163 *
164 * LDVSR (input) INTEGER
165 * The leading dimension of the matrix VSR. LDVSR >= 1, and
166 * if JOBVSR = 'V', LDVSR >= N.
167 *
168 * RCONDE (output) DOUBLE PRECISION array, dimension ( 2 )
169 * If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
170 * reciprocal condition numbers for the average of the selected
171 * eigenvalues.
172 * Not referenced if SENSE = 'N' or 'V'.
173 *
174 * RCONDV (output) DOUBLE PRECISION array, dimension ( 2 )
175 * If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
176 * reciprocal condition numbers for the selected deflating
177 * subspaces.
178 * Not referenced if SENSE = 'N' or 'E'.
179 *
180 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182 *
183 * LWORK (input) INTEGER
184 * The dimension of the array WORK.
185 * If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
186 * LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
187 * LWORK >= max( 8*N, 6*N+16 ).
188 * Note that 2*SDIM*(N-SDIM) <= N*N/2.
189 * Note also that an error is only returned if
190 * LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
191 * this may not be large enough.
192 *
193 * If LWORK = -1, then a workspace query is assumed; the routine
194 * only calculates the bound on the optimal size of the WORK
195 * array and the minimum size of the IWORK array, returns these
196 * values as the first entries of the WORK and IWORK arrays, and
197 * no error message related to LWORK or LIWORK is issued by
198 * XERBLA.
199 *
200 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
201 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
202 *
203 * LIWORK (input) INTEGER
204 * The dimension of the array IWORK.
205 * If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
206 * LIWORK >= N+6.
207 *
208 * If LIWORK = -1, then a workspace query is assumed; the
209 * routine only calculates the bound on the optimal size of the
210 * WORK array and the minimum size of the IWORK array, returns
211 * these values as the first entries of the WORK and IWORK
212 * arrays, and no error message related to LWORK or LIWORK is
213 * issued by XERBLA.
214 *
215 * BWORK (workspace) LOGICAL array, dimension (N)
216 * Not referenced if SORT = 'N'.
217 *
218 * INFO (output) INTEGER
219 * = 0: successful exit
220 * < 0: if INFO = -i, the i-th argument had an illegal value.
221 * = 1,...,N:
222 * The QZ iteration failed. (A,B) are not in Schur
223 * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
224 * be correct for j=INFO+1,...,N.
225 * > N: =N+1: other than QZ iteration failed in DHGEQZ
226 * =N+2: after reordering, roundoff changed values of
227 * some complex eigenvalues so that leading
228 * eigenvalues in the Generalized Schur form no
229 * longer satisfy SELCTG=.TRUE. This could also
230 * be caused due to scaling.
231 * =N+3: reordering failed in DTGSEN.
232 *
233 * Further Details
234 * ===============
235 *
236 * An approximate (asymptotic) bound on the average absolute error of
237 * the selected eigenvalues is
238 *
239 * EPS * norm((A, B)) / RCONDE( 1 ).
240 *
241 * An approximate (asymptotic) bound on the maximum angular error in
242 * the computed deflating subspaces is
243 *
244 * EPS * norm((A, B)) / RCONDV( 2 ).
245 *
246 * See LAPACK User's Guide, section 4.11 for more information.
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251 DOUBLE PRECISION ZERO, ONE
252 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
253 * ..
254 * .. Local Scalars ..
255 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
256 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
257 $ WANTSV
258 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
259 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
260 $ LIWMIN, LWRK, MAXWRK, MINWRK
261 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
262 $ PR, SAFMAX, SAFMIN, SMLNUM
263 * ..
264 * .. Local Arrays ..
265 DOUBLE PRECISION DIF( 2 )
266 * ..
267 * .. External Subroutines ..
268 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
269 $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
270 $ XERBLA
271 * ..
272 * .. External Functions ..
273 LOGICAL LSAME
274 INTEGER ILAENV
275 DOUBLE PRECISION DLAMCH, DLANGE
276 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
277 * ..
278 * .. Intrinsic Functions ..
279 INTRINSIC ABS, MAX, SQRT
280 * ..
281 * .. Executable Statements ..
282 *
283 * Decode the input arguments
284 *
285 IF( LSAME( JOBVSL, 'N' ) ) THEN
286 IJOBVL = 1
287 ILVSL = .FALSE.
288 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
289 IJOBVL = 2
290 ILVSL = .TRUE.
291 ELSE
292 IJOBVL = -1
293 ILVSL = .FALSE.
294 END IF
295 *
296 IF( LSAME( JOBVSR, 'N' ) ) THEN
297 IJOBVR = 1
298 ILVSR = .FALSE.
299 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
300 IJOBVR = 2
301 ILVSR = .TRUE.
302 ELSE
303 IJOBVR = -1
304 ILVSR = .FALSE.
305 END IF
306 *
307 WANTST = LSAME( SORT, 'S' )
308 WANTSN = LSAME( SENSE, 'N' )
309 WANTSE = LSAME( SENSE, 'E' )
310 WANTSV = LSAME( SENSE, 'V' )
311 WANTSB = LSAME( SENSE, 'B' )
312 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
313 IF( WANTSN ) THEN
314 IJOB = 0
315 ELSE IF( WANTSE ) THEN
316 IJOB = 1
317 ELSE IF( WANTSV ) THEN
318 IJOB = 2
319 ELSE IF( WANTSB ) THEN
320 IJOB = 4
321 END IF
322 *
323 * Test the input arguments
324 *
325 INFO = 0
326 IF( IJOBVL.LE.0 ) THEN
327 INFO = -1
328 ELSE IF( IJOBVR.LE.0 ) THEN
329 INFO = -2
330 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
331 INFO = -3
332 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
333 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
334 INFO = -5
335 ELSE IF( N.LT.0 ) THEN
336 INFO = -6
337 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338 INFO = -8
339 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340 INFO = -10
341 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
342 INFO = -16
343 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
344 INFO = -18
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.)
353 *
354 IF( INFO.EQ.0 ) THEN
355 IF( N.GT.0) THEN
356 MINWRK = MAX( 8*N, 6*N + 16 )
357 MAXWRK = MINWRK - N +
358 $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
359 MAXWRK = MAX( MAXWRK, MINWRK - N +
360 $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
361 IF( ILVSL ) THEN
362 MAXWRK = MAX( MAXWRK, MINWRK - N +
363 $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
364 END IF
365 LWRK = MAXWRK
366 IF( IJOB.GE.1 )
367 $ LWRK = MAX( LWRK, N*N/2 )
368 ELSE
369 MINWRK = 1
370 MAXWRK = 1
371 LWRK = 1
372 END IF
373 WORK( 1 ) = LWRK
374 IF( WANTSN .OR. N.EQ.0 ) THEN
375 LIWMIN = 1
376 ELSE
377 LIWMIN = N + 6
378 END IF
379 IWORK( 1 ) = LIWMIN
380 *
381 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
382 INFO = -22
383 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
384 INFO = -24
385 END IF
386 END IF
387 *
388 IF( INFO.NE.0 ) THEN
389 CALL XERBLA( 'DGGESX', -INFO )
390 RETURN
391 ELSE IF (LQUERY) THEN
392 RETURN
393 END IF
394 *
395 * Quick return if possible
396 *
397 IF( N.EQ.0 ) THEN
398 SDIM = 0
399 RETURN
400 END IF
401 *
402 * Get machine constants
403 *
404 EPS = DLAMCH( 'P' )
405 SAFMIN = DLAMCH( 'S' )
406 SAFMAX = ONE / SAFMIN
407 CALL DLABAD( SAFMIN, SAFMAX )
408 SMLNUM = SQRT( SAFMIN ) / EPS
409 BIGNUM = ONE / SMLNUM
410 *
411 * Scale A if max element outside range [SMLNUM,BIGNUM]
412 *
413 ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
414 ILASCL = .FALSE.
415 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
416 ANRMTO = SMLNUM
417 ILASCL = .TRUE.
418 ELSE IF( ANRM.GT.BIGNUM ) THEN
419 ANRMTO = BIGNUM
420 ILASCL = .TRUE.
421 END IF
422 IF( ILASCL )
423 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
424 *
425 * Scale B if max element outside range [SMLNUM,BIGNUM]
426 *
427 BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
428 ILBSCL = .FALSE.
429 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
430 BNRMTO = SMLNUM
431 ILBSCL = .TRUE.
432 ELSE IF( BNRM.GT.BIGNUM ) THEN
433 BNRMTO = BIGNUM
434 ILBSCL = .TRUE.
435 END IF
436 IF( ILBSCL )
437 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
438 *
439 * Permute the matrix to make it more nearly triangular
440 * (Workspace: need 6*N + 2*N for permutation parameters)
441 *
442 ILEFT = 1
443 IRIGHT = N + 1
444 IWRK = IRIGHT + N
445 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
446 $ WORK( IRIGHT ), WORK( IWRK ), IERR )
447 *
448 * Reduce B to triangular form (QR decomposition of B)
449 * (Workspace: need N, prefer N*NB)
450 *
451 IROWS = IHI + 1 - ILO
452 ICOLS = N + 1 - ILO
453 ITAU = IWRK
454 IWRK = ITAU + IROWS
455 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
456 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
457 *
458 * Apply the orthogonal transformation to matrix A
459 * (Workspace: need N, prefer N*NB)
460 *
461 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
462 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
463 $ LWORK+1-IWRK, IERR )
464 *
465 * Initialize VSL
466 * (Workspace: need N, prefer N*NB)
467 *
468 IF( ILVSL ) THEN
469 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
470 IF( IROWS.GT.1 ) THEN
471 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
472 $ VSL( ILO+1, ILO ), LDVSL )
473 END IF
474 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
475 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
476 END IF
477 *
478 * Initialize VSR
479 *
480 IF( ILVSR )
481 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
482 *
483 * Reduce to generalized Hessenberg form
484 * (Workspace: none needed)
485 *
486 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
487 $ LDVSL, VSR, LDVSR, IERR )
488 *
489 SDIM = 0
490 *
491 * Perform QZ algorithm, computing Schur vectors if desired
492 * (Workspace: need N)
493 *
494 IWRK = ITAU
495 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
496 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
497 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
498 IF( IERR.NE.0 ) THEN
499 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
500 INFO = IERR
501 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
502 INFO = IERR - N
503 ELSE
504 INFO = N + 1
505 END IF
506 GO TO 60
507 END IF
508 *
509 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
510 * condition number(s)
511 * (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
512 * otherwise, need 8*(N+1) )
513 *
514 IF( WANTST ) THEN
515 *
516 * Undo scaling on eigenvalues before SELCTGing
517 *
518 IF( ILASCL ) THEN
519 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
520 $ IERR )
521 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
522 $ IERR )
523 END IF
524 IF( ILBSCL )
525 $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
526 *
527 * Select eigenvalues
528 *
529 DO 10 I = 1, N
530 BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
531 10 CONTINUE
532 *
533 * Reorder eigenvalues, transform Generalized Schur vectors, and
534 * compute reciprocal condition numbers
535 *
536 CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
537 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
538 $ SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1,
539 $ IWORK, LIWORK, IERR )
540 *
541 IF( IJOB.GE.1 )
542 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
543 IF( IERR.EQ.-22 ) THEN
544 *
545 * not enough real workspace
546 *
547 INFO = -22
548 ELSE
549 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
550 RCONDE( 1 ) = PL
551 RCONDE( 2 ) = PR
552 END IF
553 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
554 RCONDV( 1 ) = DIF( 1 )
555 RCONDV( 2 ) = DIF( 2 )
556 END IF
557 IF( IERR.EQ.1 )
558 $ INFO = N + 3
559 END IF
560 *
561 END IF
562 *
563 * Apply permutation to VSL and VSR
564 * (Workspace: none needed)
565 *
566 IF( ILVSL )
567 $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
568 $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
569 *
570 IF( ILVSR )
571 $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
572 $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
573 *
574 * Check if unscaling would cause over/underflow, if so, rescale
575 * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
576 * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
577 *
578 IF( ILASCL ) THEN
579 DO 20 I = 1, N
580 IF( ALPHAI( I ).NE.ZERO ) THEN
581 IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
582 $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
583 WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
584 BETA( I ) = BETA( I )*WORK( 1 )
585 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
586 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
587 ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
588 $ ( ANRMTO / ANRM ) .OR.
589 $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
590 $ THEN
591 WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
592 BETA( I ) = BETA( I )*WORK( 1 )
593 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
594 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
595 END IF
596 END IF
597 20 CONTINUE
598 END IF
599 *
600 IF( ILBSCL ) THEN
601 DO 30 I = 1, N
602 IF( ALPHAI( I ).NE.ZERO ) THEN
603 IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
604 $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
605 WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
606 BETA( I ) = BETA( I )*WORK( 1 )
607 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
608 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
609 END IF
610 END IF
611 30 CONTINUE
612 END IF
613 *
614 * Undo scaling
615 *
616 IF( ILASCL ) THEN
617 CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
618 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
619 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
620 END IF
621 *
622 IF( ILBSCL ) THEN
623 CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
624 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
625 END IF
626 *
627 IF( WANTST ) THEN
628 *
629 * Check if reordering is correct
630 *
631 LASTSL = .TRUE.
632 LST2SL = .TRUE.
633 SDIM = 0
634 IP = 0
635 DO 50 I = 1, N
636 CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
637 IF( ALPHAI( I ).EQ.ZERO ) THEN
638 IF( CURSL )
639 $ SDIM = SDIM + 1
640 IP = 0
641 IF( CURSL .AND. .NOT.LASTSL )
642 $ INFO = N + 2
643 ELSE
644 IF( IP.EQ.1 ) THEN
645 *
646 * Last eigenvalue of conjugate pair
647 *
648 CURSL = CURSL .OR. LASTSL
649 LASTSL = CURSL
650 IF( CURSL )
651 $ SDIM = SDIM + 2
652 IP = -1
653 IF( CURSL .AND. .NOT.LST2SL )
654 $ INFO = N + 2
655 ELSE
656 *
657 * First eigenvalue of conjugate pair
658 *
659 IP = 1
660 END IF
661 END IF
662 LST2SL = LASTSL
663 LASTSL = CURSL
664 50 CONTINUE
665 *
666 END IF
667 *
668 60 CONTINUE
669 *
670 WORK( 1 ) = MAXWRK
671 IWORK( 1 ) = LIWMIN
672 *
673 RETURN
674 *
675 * End of DGGESX
676 *
677 END
2 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
3 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
4 $ LIWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2.1) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * -- April 2009 --
10 *
11 * .. Scalar Arguments ..
12 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
13 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
14 $ SDIM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
20 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
21 $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
22 $ WORK( * )
23 * ..
24 * .. Function Arguments ..
25 LOGICAL SELCTG
26 EXTERNAL SELCTG
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * DGGESX computes for a pair of N-by-N real nonsymmetric matrices
33 * (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
34 * optionally, the left and/or right matrices of Schur vectors (VSL and
35 * VSR). This gives the generalized Schur factorization
36 *
37 * (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
38 *
39 * Optionally, it also orders the eigenvalues so that a selected cluster
40 * of eigenvalues appears in the leading diagonal blocks of the upper
41 * quasi-triangular matrix S and the upper triangular matrix T; computes
42 * a reciprocal condition number for the average of the selected
43 * eigenvalues (RCONDE); and computes a reciprocal condition number for
44 * the right and left deflating subspaces corresponding to the selected
45 * eigenvalues (RCONDV). The leading columns of VSL and VSR then form
46 * an orthonormal basis for the corresponding left and right eigenspaces
47 * (deflating subspaces).
48 *
49 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
50 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
51 * usually represented as the pair (alpha,beta), as there is a
52 * reasonable interpretation for beta=0 or for both being zero.
53 *
54 * A pair of matrices (S,T) is in generalized real Schur form if T is
55 * upper triangular with non-negative diagonal and S is block upper
56 * triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
57 * to real generalized eigenvalues, while 2-by-2 blocks of S will be
58 * "standardized" by making the corresponding elements of T have the
59 * form:
60 * [ a 0 ]
61 * [ 0 b ]
62 *
63 * and the pair of corresponding 2-by-2 blocks in S and T will have a
64 * complex conjugate pair of generalized eigenvalues.
65 *
66 *
67 * Arguments
68 * =========
69 *
70 * JOBVSL (input) CHARACTER*1
71 * = 'N': do not compute the left Schur vectors;
72 * = 'V': compute the left Schur vectors.
73 *
74 * JOBVSR (input) CHARACTER*1
75 * = 'N': do not compute the right Schur vectors;
76 * = 'V': compute the right Schur vectors.
77 *
78 * SORT (input) CHARACTER*1
79 * Specifies whether or not to order the eigenvalues on the
80 * diagonal of the generalized Schur form.
81 * = 'N': Eigenvalues are not ordered;
82 * = 'S': Eigenvalues are ordered (see SELCTG).
83 *
84 * SELCTG (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
85 * SELCTG must be declared EXTERNAL in the calling subroutine.
86 * If SORT = 'N', SELCTG is not referenced.
87 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
88 * to the top left of the Schur form.
89 * An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
90 * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
91 * one of a complex conjugate pair of eigenvalues is selected,
92 * then both complex eigenvalues are selected.
93 * Note that a selected complex eigenvalue may no longer satisfy
94 * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
95 * since ordering may change the value of complex eigenvalues
96 * (especially if the eigenvalue is ill-conditioned), in this
97 * case INFO is set to N+3.
98 *
99 * SENSE (input) CHARACTER*1
100 * Determines which reciprocal condition numbers are computed.
101 * = 'N' : None are computed;
102 * = 'E' : Computed for average of selected eigenvalues only;
103 * = 'V' : Computed for selected deflating subspaces only;
104 * = 'B' : Computed for both.
105 * If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
106 *
107 * N (input) INTEGER
108 * The order of the matrices A, B, VSL, and VSR. N >= 0.
109 *
110 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
111 * On entry, the first of the pair of matrices.
112 * On exit, A has been overwritten by its generalized Schur
113 * form S.
114 *
115 * LDA (input) INTEGER
116 * The leading dimension of A. LDA >= max(1,N).
117 *
118 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
119 * On entry, the second of the pair of matrices.
120 * On exit, B has been overwritten by its generalized Schur
121 * form T.
122 *
123 * LDB (input) INTEGER
124 * The leading dimension of B. LDB >= max(1,N).
125 *
126 * SDIM (output) INTEGER
127 * If SORT = 'N', SDIM = 0.
128 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
129 * for which SELCTG is true. (Complex conjugate pairs for which
130 * SELCTG is true for either eigenvalue count as 2.)
131 *
132 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
133 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
134 * BETA (output) DOUBLE PRECISION array, dimension (N)
135 * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
136 * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
137 * and BETA(j),j=1,...,N are the diagonals of the complex Schur
138 * form (S,T) that would result if the 2-by-2 diagonal blocks of
139 * the real Schur form of (A,B) were further reduced to
140 * triangular form using 2-by-2 complex unitary transformations.
141 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
142 * positive, then the j-th and (j+1)-st eigenvalues are a
143 * complex conjugate pair, with ALPHAI(j+1) negative.
144 *
145 * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
146 * may easily over- or underflow, and BETA(j) may even be zero.
147 * Thus, the user should avoid naively computing the ratio.
148 * However, ALPHAR and ALPHAI will be always less than and
149 * usually comparable with norm(A) in magnitude, and BETA always
150 * less than and usually comparable with norm(B).
151 *
152 * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
153 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
154 * Not referenced if JOBVSL = 'N'.
155 *
156 * LDVSL (input) INTEGER
157 * The leading dimension of the matrix VSL. LDVSL >=1, and
158 * if JOBVSL = 'V', LDVSL >= N.
159 *
160 * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
161 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
162 * Not referenced if JOBVSR = 'N'.
163 *
164 * LDVSR (input) INTEGER
165 * The leading dimension of the matrix VSR. LDVSR >= 1, and
166 * if JOBVSR = 'V', LDVSR >= N.
167 *
168 * RCONDE (output) DOUBLE PRECISION array, dimension ( 2 )
169 * If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
170 * reciprocal condition numbers for the average of the selected
171 * eigenvalues.
172 * Not referenced if SENSE = 'N' or 'V'.
173 *
174 * RCONDV (output) DOUBLE PRECISION array, dimension ( 2 )
175 * If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
176 * reciprocal condition numbers for the selected deflating
177 * subspaces.
178 * Not referenced if SENSE = 'N' or 'E'.
179 *
180 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182 *
183 * LWORK (input) INTEGER
184 * The dimension of the array WORK.
185 * If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
186 * LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
187 * LWORK >= max( 8*N, 6*N+16 ).
188 * Note that 2*SDIM*(N-SDIM) <= N*N/2.
189 * Note also that an error is only returned if
190 * LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
191 * this may not be large enough.
192 *
193 * If LWORK = -1, then a workspace query is assumed; the routine
194 * only calculates the bound on the optimal size of the WORK
195 * array and the minimum size of the IWORK array, returns these
196 * values as the first entries of the WORK and IWORK arrays, and
197 * no error message related to LWORK or LIWORK is issued by
198 * XERBLA.
199 *
200 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
201 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
202 *
203 * LIWORK (input) INTEGER
204 * The dimension of the array IWORK.
205 * If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
206 * LIWORK >= N+6.
207 *
208 * If LIWORK = -1, then a workspace query is assumed; the
209 * routine only calculates the bound on the optimal size of the
210 * WORK array and the minimum size of the IWORK array, returns
211 * these values as the first entries of the WORK and IWORK
212 * arrays, and no error message related to LWORK or LIWORK is
213 * issued by XERBLA.
214 *
215 * BWORK (workspace) LOGICAL array, dimension (N)
216 * Not referenced if SORT = 'N'.
217 *
218 * INFO (output) INTEGER
219 * = 0: successful exit
220 * < 0: if INFO = -i, the i-th argument had an illegal value.
221 * = 1,...,N:
222 * The QZ iteration failed. (A,B) are not in Schur
223 * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
224 * be correct for j=INFO+1,...,N.
225 * > N: =N+1: other than QZ iteration failed in DHGEQZ
226 * =N+2: after reordering, roundoff changed values of
227 * some complex eigenvalues so that leading
228 * eigenvalues in the Generalized Schur form no
229 * longer satisfy SELCTG=.TRUE. This could also
230 * be caused due to scaling.
231 * =N+3: reordering failed in DTGSEN.
232 *
233 * Further Details
234 * ===============
235 *
236 * An approximate (asymptotic) bound on the average absolute error of
237 * the selected eigenvalues is
238 *
239 * EPS * norm((A, B)) / RCONDE( 1 ).
240 *
241 * An approximate (asymptotic) bound on the maximum angular error in
242 * the computed deflating subspaces is
243 *
244 * EPS * norm((A, B)) / RCONDV( 2 ).
245 *
246 * See LAPACK User's Guide, section 4.11 for more information.
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251 DOUBLE PRECISION ZERO, ONE
252 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
253 * ..
254 * .. Local Scalars ..
255 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
256 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
257 $ WANTSV
258 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
259 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
260 $ LIWMIN, LWRK, MAXWRK, MINWRK
261 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
262 $ PR, SAFMAX, SAFMIN, SMLNUM
263 * ..
264 * .. Local Arrays ..
265 DOUBLE PRECISION DIF( 2 )
266 * ..
267 * .. External Subroutines ..
268 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
269 $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
270 $ XERBLA
271 * ..
272 * .. External Functions ..
273 LOGICAL LSAME
274 INTEGER ILAENV
275 DOUBLE PRECISION DLAMCH, DLANGE
276 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
277 * ..
278 * .. Intrinsic Functions ..
279 INTRINSIC ABS, MAX, SQRT
280 * ..
281 * .. Executable Statements ..
282 *
283 * Decode the input arguments
284 *
285 IF( LSAME( JOBVSL, 'N' ) ) THEN
286 IJOBVL = 1
287 ILVSL = .FALSE.
288 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
289 IJOBVL = 2
290 ILVSL = .TRUE.
291 ELSE
292 IJOBVL = -1
293 ILVSL = .FALSE.
294 END IF
295 *
296 IF( LSAME( JOBVSR, 'N' ) ) THEN
297 IJOBVR = 1
298 ILVSR = .FALSE.
299 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
300 IJOBVR = 2
301 ILVSR = .TRUE.
302 ELSE
303 IJOBVR = -1
304 ILVSR = .FALSE.
305 END IF
306 *
307 WANTST = LSAME( SORT, 'S' )
308 WANTSN = LSAME( SENSE, 'N' )
309 WANTSE = LSAME( SENSE, 'E' )
310 WANTSV = LSAME( SENSE, 'V' )
311 WANTSB = LSAME( SENSE, 'B' )
312 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
313 IF( WANTSN ) THEN
314 IJOB = 0
315 ELSE IF( WANTSE ) THEN
316 IJOB = 1
317 ELSE IF( WANTSV ) THEN
318 IJOB = 2
319 ELSE IF( WANTSB ) THEN
320 IJOB = 4
321 END IF
322 *
323 * Test the input arguments
324 *
325 INFO = 0
326 IF( IJOBVL.LE.0 ) THEN
327 INFO = -1
328 ELSE IF( IJOBVR.LE.0 ) THEN
329 INFO = -2
330 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
331 INFO = -3
332 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
333 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
334 INFO = -5
335 ELSE IF( N.LT.0 ) THEN
336 INFO = -6
337 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338 INFO = -8
339 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340 INFO = -10
341 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
342 INFO = -16
343 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
344 INFO = -18
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.)
353 *
354 IF( INFO.EQ.0 ) THEN
355 IF( N.GT.0) THEN
356 MINWRK = MAX( 8*N, 6*N + 16 )
357 MAXWRK = MINWRK - N +
358 $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
359 MAXWRK = MAX( MAXWRK, MINWRK - N +
360 $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
361 IF( ILVSL ) THEN
362 MAXWRK = MAX( MAXWRK, MINWRK - N +
363 $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
364 END IF
365 LWRK = MAXWRK
366 IF( IJOB.GE.1 )
367 $ LWRK = MAX( LWRK, N*N/2 )
368 ELSE
369 MINWRK = 1
370 MAXWRK = 1
371 LWRK = 1
372 END IF
373 WORK( 1 ) = LWRK
374 IF( WANTSN .OR. N.EQ.0 ) THEN
375 LIWMIN = 1
376 ELSE
377 LIWMIN = N + 6
378 END IF
379 IWORK( 1 ) = LIWMIN
380 *
381 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
382 INFO = -22
383 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
384 INFO = -24
385 END IF
386 END IF
387 *
388 IF( INFO.NE.0 ) THEN
389 CALL XERBLA( 'DGGESX', -INFO )
390 RETURN
391 ELSE IF (LQUERY) THEN
392 RETURN
393 END IF
394 *
395 * Quick return if possible
396 *
397 IF( N.EQ.0 ) THEN
398 SDIM = 0
399 RETURN
400 END IF
401 *
402 * Get machine constants
403 *
404 EPS = DLAMCH( 'P' )
405 SAFMIN = DLAMCH( 'S' )
406 SAFMAX = ONE / SAFMIN
407 CALL DLABAD( SAFMIN, SAFMAX )
408 SMLNUM = SQRT( SAFMIN ) / EPS
409 BIGNUM = ONE / SMLNUM
410 *
411 * Scale A if max element outside range [SMLNUM,BIGNUM]
412 *
413 ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
414 ILASCL = .FALSE.
415 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
416 ANRMTO = SMLNUM
417 ILASCL = .TRUE.
418 ELSE IF( ANRM.GT.BIGNUM ) THEN
419 ANRMTO = BIGNUM
420 ILASCL = .TRUE.
421 END IF
422 IF( ILASCL )
423 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
424 *
425 * Scale B if max element outside range [SMLNUM,BIGNUM]
426 *
427 BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
428 ILBSCL = .FALSE.
429 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
430 BNRMTO = SMLNUM
431 ILBSCL = .TRUE.
432 ELSE IF( BNRM.GT.BIGNUM ) THEN
433 BNRMTO = BIGNUM
434 ILBSCL = .TRUE.
435 END IF
436 IF( ILBSCL )
437 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
438 *
439 * Permute the matrix to make it more nearly triangular
440 * (Workspace: need 6*N + 2*N for permutation parameters)
441 *
442 ILEFT = 1
443 IRIGHT = N + 1
444 IWRK = IRIGHT + N
445 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
446 $ WORK( IRIGHT ), WORK( IWRK ), IERR )
447 *
448 * Reduce B to triangular form (QR decomposition of B)
449 * (Workspace: need N, prefer N*NB)
450 *
451 IROWS = IHI + 1 - ILO
452 ICOLS = N + 1 - ILO
453 ITAU = IWRK
454 IWRK = ITAU + IROWS
455 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
456 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
457 *
458 * Apply the orthogonal transformation to matrix A
459 * (Workspace: need N, prefer N*NB)
460 *
461 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
462 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
463 $ LWORK+1-IWRK, IERR )
464 *
465 * Initialize VSL
466 * (Workspace: need N, prefer N*NB)
467 *
468 IF( ILVSL ) THEN
469 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
470 IF( IROWS.GT.1 ) THEN
471 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
472 $ VSL( ILO+1, ILO ), LDVSL )
473 END IF
474 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
475 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
476 END IF
477 *
478 * Initialize VSR
479 *
480 IF( ILVSR )
481 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
482 *
483 * Reduce to generalized Hessenberg form
484 * (Workspace: none needed)
485 *
486 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
487 $ LDVSL, VSR, LDVSR, IERR )
488 *
489 SDIM = 0
490 *
491 * Perform QZ algorithm, computing Schur vectors if desired
492 * (Workspace: need N)
493 *
494 IWRK = ITAU
495 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
496 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
497 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
498 IF( IERR.NE.0 ) THEN
499 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
500 INFO = IERR
501 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
502 INFO = IERR - N
503 ELSE
504 INFO = N + 1
505 END IF
506 GO TO 60
507 END IF
508 *
509 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
510 * condition number(s)
511 * (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
512 * otherwise, need 8*(N+1) )
513 *
514 IF( WANTST ) THEN
515 *
516 * Undo scaling on eigenvalues before SELCTGing
517 *
518 IF( ILASCL ) THEN
519 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
520 $ IERR )
521 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
522 $ IERR )
523 END IF
524 IF( ILBSCL )
525 $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
526 *
527 * Select eigenvalues
528 *
529 DO 10 I = 1, N
530 BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
531 10 CONTINUE
532 *
533 * Reorder eigenvalues, transform Generalized Schur vectors, and
534 * compute reciprocal condition numbers
535 *
536 CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
537 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
538 $ SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1,
539 $ IWORK, LIWORK, IERR )
540 *
541 IF( IJOB.GE.1 )
542 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
543 IF( IERR.EQ.-22 ) THEN
544 *
545 * not enough real workspace
546 *
547 INFO = -22
548 ELSE
549 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
550 RCONDE( 1 ) = PL
551 RCONDE( 2 ) = PR
552 END IF
553 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
554 RCONDV( 1 ) = DIF( 1 )
555 RCONDV( 2 ) = DIF( 2 )
556 END IF
557 IF( IERR.EQ.1 )
558 $ INFO = N + 3
559 END IF
560 *
561 END IF
562 *
563 * Apply permutation to VSL and VSR
564 * (Workspace: none needed)
565 *
566 IF( ILVSL )
567 $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
568 $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
569 *
570 IF( ILVSR )
571 $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
572 $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
573 *
574 * Check if unscaling would cause over/underflow, if so, rescale
575 * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
576 * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
577 *
578 IF( ILASCL ) THEN
579 DO 20 I = 1, N
580 IF( ALPHAI( I ).NE.ZERO ) THEN
581 IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
582 $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
583 WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
584 BETA( I ) = BETA( I )*WORK( 1 )
585 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
586 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
587 ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
588 $ ( ANRMTO / ANRM ) .OR.
589 $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
590 $ THEN
591 WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
592 BETA( I ) = BETA( I )*WORK( 1 )
593 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
594 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
595 END IF
596 END IF
597 20 CONTINUE
598 END IF
599 *
600 IF( ILBSCL ) THEN
601 DO 30 I = 1, N
602 IF( ALPHAI( I ).NE.ZERO ) THEN
603 IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
604 $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
605 WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
606 BETA( I ) = BETA( I )*WORK( 1 )
607 ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
608 ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
609 END IF
610 END IF
611 30 CONTINUE
612 END IF
613 *
614 * Undo scaling
615 *
616 IF( ILASCL ) THEN
617 CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
618 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
619 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
620 END IF
621 *
622 IF( ILBSCL ) THEN
623 CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
624 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
625 END IF
626 *
627 IF( WANTST ) THEN
628 *
629 * Check if reordering is correct
630 *
631 LASTSL = .TRUE.
632 LST2SL = .TRUE.
633 SDIM = 0
634 IP = 0
635 DO 50 I = 1, N
636 CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
637 IF( ALPHAI( I ).EQ.ZERO ) THEN
638 IF( CURSL )
639 $ SDIM = SDIM + 1
640 IP = 0
641 IF( CURSL .AND. .NOT.LASTSL )
642 $ INFO = N + 2
643 ELSE
644 IF( IP.EQ.1 ) THEN
645 *
646 * Last eigenvalue of conjugate pair
647 *
648 CURSL = CURSL .OR. LASTSL
649 LASTSL = CURSL
650 IF( CURSL )
651 $ SDIM = SDIM + 2
652 IP = -1
653 IF( CURSL .AND. .NOT.LST2SL )
654 $ INFO = N + 2
655 ELSE
656 *
657 * First eigenvalue of conjugate pair
658 *
659 IP = 1
660 END IF
661 END IF
662 LST2SL = LASTSL
663 LASTSL = CURSL
664 50 CONTINUE
665 *
666 END IF
667 *
668 60 CONTINUE
669 *
670 WORK( 1 ) = MAXWRK
671 IWORK( 1 ) = LIWMIN
672 *
673 RETURN
674 *
675 * End of DGGESX
676 *
677 END