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