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