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: *> \date June 2017
341: *
342: *> \ingroup doubleGEeigen
343: *
344: *> \par Further Details:
345: * =====================
346: *>
347: *> \verbatim
348: *>
349: *> An approximate (asymptotic) bound on the average absolute error of
350: *> the selected eigenvalues is
351: *>
352: *> EPS * norm((A, B)) / RCONDE( 1 ).
353: *>
354: *> An approximate (asymptotic) bound on the maximum angular error in
355: *> the computed deflating subspaces is
356: *>
357: *> EPS * norm((A, B)) / RCONDV( 2 ).
358: *>
359: *> See LAPACK User's Guide, section 4.11 for more information.
360: *> \endverbatim
361: *>
362: * =====================================================================
363: SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
364: $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
365: $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
366: $ LIWORK, BWORK, INFO )
367: *
368: * -- LAPACK driver routine (version 3.7.1) --
369: * -- LAPACK is a software package provided by Univ. of Tennessee, --
370: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
371: * June 2017
372: *
373: * .. Scalar Arguments ..
374: CHARACTER JOBVSL, JOBVSR, SENSE, SORT
375: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
376: $ SDIM
377: * ..
378: * .. Array Arguments ..
379: LOGICAL BWORK( * )
380: INTEGER IWORK( * )
381: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
382: $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
383: $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
384: $ WORK( * )
385: * ..
386: * .. Function Arguments ..
387: LOGICAL SELCTG
388: EXTERNAL SELCTG
389: * ..
390: *
391: * =====================================================================
392: *
393: * .. Parameters ..
394: DOUBLE PRECISION ZERO, ONE
395: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
396: * ..
397: * .. Local Scalars ..
398: LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
399: $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
400: $ WANTSV
401: INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
402: $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
403: $ LIWMIN, LWRK, MAXWRK, MINWRK
404: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
405: $ PR, SAFMAX, SAFMIN, SMLNUM
406: * ..
407: * .. Local Arrays ..
408: DOUBLE PRECISION DIF( 2 )
409: * ..
410: * .. External Subroutines ..
411: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
412: $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
413: $ XERBLA
414: * ..
415: * .. External Functions ..
416: LOGICAL LSAME
417: INTEGER ILAENV
418: DOUBLE PRECISION DLAMCH, DLANGE
419: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
420: * ..
421: * .. Intrinsic Functions ..
422: INTRINSIC ABS, MAX, SQRT
423: * ..
424: * .. Executable Statements ..
425: *
426: * Decode the input arguments
427: *
428: IF( LSAME( JOBVSL, 'N' ) ) THEN
429: IJOBVL = 1
430: ILVSL = .FALSE.
431: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
432: IJOBVL = 2
433: ILVSL = .TRUE.
434: ELSE
435: IJOBVL = -1
436: ILVSL = .FALSE.
437: END IF
438: *
439: IF( LSAME( JOBVSR, 'N' ) ) THEN
440: IJOBVR = 1
441: ILVSR = .FALSE.
442: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
443: IJOBVR = 2
444: ILVSR = .TRUE.
445: ELSE
446: IJOBVR = -1
447: ILVSR = .FALSE.
448: END IF
449: *
450: WANTST = LSAME( SORT, 'S' )
451: WANTSN = LSAME( SENSE, 'N' )
452: WANTSE = LSAME( SENSE, 'E' )
453: WANTSV = LSAME( SENSE, 'V' )
454: WANTSB = LSAME( SENSE, 'B' )
455: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
456: IF( WANTSN ) THEN
457: IJOB = 0
458: ELSE IF( WANTSE ) THEN
459: IJOB = 1
460: ELSE IF( WANTSV ) THEN
461: IJOB = 2
462: ELSE IF( WANTSB ) THEN
463: IJOB = 4
464: END IF
465: *
466: * Test the input arguments
467: *
468: INFO = 0
469: IF( IJOBVL.LE.0 ) THEN
470: INFO = -1
471: ELSE IF( IJOBVR.LE.0 ) THEN
472: INFO = -2
473: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
474: INFO = -3
475: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
476: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
477: INFO = -5
478: ELSE IF( N.LT.0 ) THEN
479: INFO = -6
480: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
481: INFO = -8
482: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
483: INFO = -10
484: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
485: INFO = -16
486: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
487: INFO = -18
488: END IF
489: *
490: * Compute workspace
491: * (Note: Comments in the code beginning "Workspace:" describe the
492: * minimal amount of workspace needed at that point in the code,
493: * as well as the preferred amount for good performance.
494: * NB refers to the optimal block size for the immediately
495: * following subroutine, as returned by ILAENV.)
496: *
497: IF( INFO.EQ.0 ) THEN
498: IF( N.GT.0) THEN
499: MINWRK = MAX( 8*N, 6*N + 16 )
500: MAXWRK = MINWRK - N +
501: $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
502: MAXWRK = MAX( MAXWRK, MINWRK - N +
503: $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
504: IF( ILVSL ) THEN
505: MAXWRK = MAX( MAXWRK, MINWRK - N +
506: $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
507: END IF
508: LWRK = MAXWRK
509: IF( IJOB.GE.1 )
510: $ LWRK = MAX( LWRK, N*N/2 )
511: ELSE
512: MINWRK = 1
513: MAXWRK = 1
514: LWRK = 1
515: END IF
516: WORK( 1 ) = LWRK
517: IF( WANTSN .OR. N.EQ.0 ) THEN
518: LIWMIN = 1
519: ELSE
520: LIWMIN = N + 6
521: END IF
522: IWORK( 1 ) = LIWMIN
523: *
524: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
525: INFO = -22
526: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
527: INFO = -24
528: END IF
529: END IF
530: *
531: IF( INFO.NE.0 ) THEN
532: CALL XERBLA( 'DGGESX', -INFO )
533: RETURN
534: ELSE IF (LQUERY) THEN
535: RETURN
536: END IF
537: *
538: * Quick return if possible
539: *
540: IF( N.EQ.0 ) THEN
541: SDIM = 0
542: RETURN
543: END IF
544: *
545: * Get machine constants
546: *
547: EPS = DLAMCH( 'P' )
548: SAFMIN = DLAMCH( 'S' )
549: SAFMAX = ONE / SAFMIN
550: CALL DLABAD( SAFMIN, SAFMAX )
551: SMLNUM = SQRT( SAFMIN ) / EPS
552: BIGNUM = ONE / SMLNUM
553: *
554: * Scale A if max element outside range [SMLNUM,BIGNUM]
555: *
556: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
557: ILASCL = .FALSE.
558: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
559: ANRMTO = SMLNUM
560: ILASCL = .TRUE.
561: ELSE IF( ANRM.GT.BIGNUM ) THEN
562: ANRMTO = BIGNUM
563: ILASCL = .TRUE.
564: END IF
565: IF( ILASCL )
566: $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
567: *
568: * Scale B if max element outside range [SMLNUM,BIGNUM]
569: *
570: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
571: ILBSCL = .FALSE.
572: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
573: BNRMTO = SMLNUM
574: ILBSCL = .TRUE.
575: ELSE IF( BNRM.GT.BIGNUM ) THEN
576: BNRMTO = BIGNUM
577: ILBSCL = .TRUE.
578: END IF
579: IF( ILBSCL )
580: $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
581: *
582: * Permute the matrix to make it more nearly triangular
583: * (Workspace: need 6*N + 2*N for permutation parameters)
584: *
585: ILEFT = 1
586: IRIGHT = N + 1
587: IWRK = IRIGHT + N
588: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
589: $ WORK( IRIGHT ), WORK( IWRK ), IERR )
590: *
591: * Reduce B to triangular form (QR decomposition of B)
592: * (Workspace: need N, prefer N*NB)
593: *
594: IROWS = IHI + 1 - ILO
595: ICOLS = N + 1 - ILO
596: ITAU = IWRK
597: IWRK = ITAU + IROWS
598: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
599: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
600: *
601: * Apply the orthogonal transformation to matrix A
602: * (Workspace: need N, prefer N*NB)
603: *
604: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
605: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
606: $ LWORK+1-IWRK, IERR )
607: *
608: * Initialize VSL
609: * (Workspace: need N, prefer N*NB)
610: *
611: IF( ILVSL ) THEN
612: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
613: IF( IROWS.GT.1 ) THEN
614: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
615: $ VSL( ILO+1, ILO ), LDVSL )
616: END IF
617: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
618: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
619: END IF
620: *
621: * Initialize VSR
622: *
623: IF( ILVSR )
624: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
625: *
626: * Reduce to generalized Hessenberg form
627: * (Workspace: none needed)
628: *
629: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
630: $ LDVSL, VSR, LDVSR, IERR )
631: *
632: SDIM = 0
633: *
634: * Perform QZ algorithm, computing Schur vectors if desired
635: * (Workspace: need N)
636: *
637: IWRK = ITAU
638: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
639: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
640: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
641: IF( IERR.NE.0 ) THEN
642: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
643: INFO = IERR
644: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
645: INFO = IERR - N
646: ELSE
647: INFO = N + 1
648: END IF
649: GO TO 60
650: END IF
651: *
652: * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
653: * condition number(s)
654: * (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
655: * otherwise, need 8*(N+1) )
656: *
657: IF( WANTST ) THEN
658: *
659: * Undo scaling on eigenvalues before SELCTGing
660: *
661: IF( ILASCL ) THEN
662: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
663: $ IERR )
664: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
665: $ IERR )
666: END IF
667: IF( ILBSCL )
668: $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
669: *
670: * Select eigenvalues
671: *
672: DO 10 I = 1, N
673: BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
674: 10 CONTINUE
675: *
676: * Reorder eigenvalues, transform Generalized Schur vectors, and
677: * compute reciprocal condition numbers
678: *
679: CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
680: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
681: $ SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1,
682: $ IWORK, LIWORK, IERR )
683: *
684: IF( IJOB.GE.1 )
685: $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
686: IF( IERR.EQ.-22 ) THEN
687: *
688: * not enough real workspace
689: *
690: INFO = -22
691: ELSE
692: IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
693: RCONDE( 1 ) = PL
694: RCONDE( 2 ) = PR
695: END IF
696: IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
697: RCONDV( 1 ) = DIF( 1 )
698: RCONDV( 2 ) = DIF( 2 )
699: END IF
700: IF( IERR.EQ.1 )
701: $ INFO = N + 3
702: END IF
703: *
704: END IF
705: *
706: * Apply permutation to VSL and VSR
707: * (Workspace: none needed)
708: *
709: IF( ILVSL )
710: $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
711: $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
712: *
713: IF( ILVSR )
714: $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
715: $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
716: *
717: * Check if unscaling would cause over/underflow, if so, rescale
718: * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
719: * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
720: *
721: IF( ILASCL ) THEN
722: DO 20 I = 1, N
723: IF( ALPHAI( I ).NE.ZERO ) THEN
724: IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
725: $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
726: WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
727: BETA( I ) = BETA( I )*WORK( 1 )
728: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
729: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
730: ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
731: $ ( ANRMTO / ANRM ) .OR.
732: $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
733: $ THEN
734: WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
735: BETA( I ) = BETA( I )*WORK( 1 )
736: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
737: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
738: END IF
739: END IF
740: 20 CONTINUE
741: END IF
742: *
743: IF( ILBSCL ) THEN
744: DO 30 I = 1, N
745: IF( ALPHAI( I ).NE.ZERO ) THEN
746: IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
747: $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
748: WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
749: BETA( I ) = BETA( I )*WORK( 1 )
750: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
751: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
752: END IF
753: END IF
754: 30 CONTINUE
755: END IF
756: *
757: * Undo scaling
758: *
759: IF( ILASCL ) THEN
760: CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
761: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
762: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
763: END IF
764: *
765: IF( ILBSCL ) THEN
766: CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
767: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
768: END IF
769: *
770: IF( WANTST ) THEN
771: *
772: * Check if reordering is correct
773: *
774: LASTSL = .TRUE.
775: LST2SL = .TRUE.
776: SDIM = 0
777: IP = 0
778: DO 50 I = 1, N
779: CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
780: IF( ALPHAI( I ).EQ.ZERO ) THEN
781: IF( CURSL )
782: $ SDIM = SDIM + 1
783: IP = 0
784: IF( CURSL .AND. .NOT.LASTSL )
785: $ INFO = N + 2
786: ELSE
787: IF( IP.EQ.1 ) THEN
788: *
789: * Last eigenvalue of conjugate pair
790: *
791: CURSL = CURSL .OR. LASTSL
792: LASTSL = CURSL
793: IF( CURSL )
794: $ SDIM = SDIM + 2
795: IP = -1
796: IF( CURSL .AND. .NOT.LST2SL )
797: $ INFO = N + 2
798: ELSE
799: *
800: * First eigenvalue of conjugate pair
801: *
802: IP = 1
803: END IF
804: END IF
805: LST2SL = LASTSL
806: LASTSL = CURSL
807: 50 CONTINUE
808: *
809: END IF
810: *
811: 60 CONTINUE
812: *
813: WORK( 1 ) = MAXWRK
814: IWORK( 1 ) = LIWMIN
815: *
816: RETURN
817: *
818: * End of DGGESX
819: *
820: 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