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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack