Annotation of rpl/lapack/lapack/zggesx.f, revision 1.16
1.8 bertrand 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: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.14 bertrand 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">
1.8 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.8 bertrand 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 )
1.14 bertrand 25: *
1.8 bertrand 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: * ..
1.14 bertrand 43: *
1.8 bertrand 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
1.16 ! bertrand 107: *> SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments
1.8 bertrand 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: *
1.14 bertrand 318: *> \author Univ. of Tennessee
319: *> \author Univ. of California Berkeley
320: *> \author Univ. of Colorado Denver
321: *> \author NAG Ltd.
1.8 bertrand 322: *
1.16 ! bertrand 323: *> \date June 2017
1.8 bertrand 324: *
325: *> \ingroup complex16GEeigen
326: *
327: * =====================================================================
1.1 bertrand 328: SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
329: $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
330: $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
331: $ IWORK, LIWORK, BWORK, INFO )
332: *
1.16 ! bertrand 333: * -- LAPACK driver routine (version 3.7.1) --
1.1 bertrand 334: * -- LAPACK is a software package provided by Univ. of Tennessee, --
335: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 ! bertrand 336: * June 2017
1.1 bertrand 337: *
338: * .. Scalar Arguments ..
339: CHARACTER JOBVSL, JOBVSR, SENSE, SORT
340: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
341: $ SDIM
342: * ..
343: * .. Array Arguments ..
344: LOGICAL BWORK( * )
345: INTEGER IWORK( * )
346: DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
347: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
348: $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
349: $ WORK( * )
350: * ..
351: * .. Function Arguments ..
352: LOGICAL SELCTG
353: EXTERNAL SELCTG
354: * ..
355: *
356: * =====================================================================
357: *
358: * .. Parameters ..
359: DOUBLE PRECISION ZERO, ONE
360: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
361: COMPLEX*16 CZERO, CONE
362: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
363: $ CONE = ( 1.0D+0, 0.0D+0 ) )
364: * ..
365: * .. Local Scalars ..
366: LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
367: $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV
368: INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
369: $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK,
370: $ LIWMIN, LWRK, MAXWRK, MINWRK
371: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
372: $ PR, SMLNUM
373: * ..
374: * .. Local Arrays ..
375: DOUBLE PRECISION DIF( 2 )
376: * ..
377: * .. External Subroutines ..
378: EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
379: $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
380: $ ZUNMQR
381: * ..
382: * .. External Functions ..
383: LOGICAL LSAME
384: INTEGER ILAENV
385: DOUBLE PRECISION DLAMCH, ZLANGE
386: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
387: * ..
388: * .. Intrinsic Functions ..
389: INTRINSIC MAX, SQRT
390: * ..
391: * .. Executable Statements ..
392: *
393: * Decode the input arguments
394: *
395: IF( LSAME( JOBVSL, 'N' ) ) THEN
396: IJOBVL = 1
397: ILVSL = .FALSE.
398: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
399: IJOBVL = 2
400: ILVSL = .TRUE.
401: ELSE
402: IJOBVL = -1
403: ILVSL = .FALSE.
404: END IF
405: *
406: IF( LSAME( JOBVSR, 'N' ) ) THEN
407: IJOBVR = 1
408: ILVSR = .FALSE.
409: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
410: IJOBVR = 2
411: ILVSR = .TRUE.
412: ELSE
413: IJOBVR = -1
414: ILVSR = .FALSE.
415: END IF
416: *
417: WANTST = LSAME( SORT, 'S' )
418: WANTSN = LSAME( SENSE, 'N' )
419: WANTSE = LSAME( SENSE, 'E' )
420: WANTSV = LSAME( SENSE, 'V' )
421: WANTSB = LSAME( SENSE, 'B' )
422: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
423: IF( WANTSN ) THEN
424: IJOB = 0
425: ELSE IF( WANTSE ) THEN
426: IJOB = 1
427: ELSE IF( WANTSV ) THEN
428: IJOB = 2
429: ELSE IF( WANTSB ) THEN
430: IJOB = 4
431: END IF
432: *
433: * Test the input arguments
434: *
435: INFO = 0
436: IF( IJOBVL.LE.0 ) THEN
437: INFO = -1
438: ELSE IF( IJOBVR.LE.0 ) THEN
439: INFO = -2
440: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
441: INFO = -3
442: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
443: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
444: INFO = -5
445: ELSE IF( N.LT.0 ) THEN
446: INFO = -6
447: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
448: INFO = -8
449: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
450: INFO = -10
451: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
452: INFO = -15
453: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
454: INFO = -17
455: END IF
456: *
457: * Compute workspace
458: * (Note: Comments in the code beginning "Workspace:" describe the
459: * minimal amount of workspace needed at that point in the code,
460: * as well as the preferred amount for good performance.
461: * NB refers to the optimal block size for the immediately
462: * following subroutine, as returned by ILAENV.)
463: *
464: IF( INFO.EQ.0 ) THEN
465: IF( N.GT.0) THEN
466: MINWRK = 2*N
467: MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
468: MAXWRK = MAX( MAXWRK, N*( 1 +
469: $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) )
470: IF( ILVSL ) THEN
471: MAXWRK = MAX( MAXWRK, N*( 1 +
472: $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) )
473: END IF
474: LWRK = MAXWRK
475: IF( IJOB.GE.1 )
476: $ LWRK = MAX( LWRK, N*N/2 )
477: ELSE
478: MINWRK = 1
479: MAXWRK = 1
480: LWRK = 1
481: END IF
482: WORK( 1 ) = LWRK
483: IF( WANTSN .OR. N.EQ.0 ) THEN
484: LIWMIN = 1
485: ELSE
486: LIWMIN = N + 2
487: END IF
488: IWORK( 1 ) = LIWMIN
489: *
490: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
491: INFO = -21
492: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN
493: INFO = -24
494: END IF
495: END IF
496: *
497: IF( INFO.NE.0 ) THEN
498: CALL XERBLA( 'ZGGESX', -INFO )
499: RETURN
500: ELSE IF (LQUERY) THEN
501: RETURN
502: END IF
503: *
504: * Quick return if possible
505: *
506: IF( N.EQ.0 ) THEN
507: SDIM = 0
508: RETURN
509: END IF
510: *
511: * Get machine constants
512: *
513: EPS = DLAMCH( 'P' )
514: SMLNUM = DLAMCH( 'S' )
515: BIGNUM = ONE / SMLNUM
516: CALL DLABAD( SMLNUM, BIGNUM )
517: SMLNUM = SQRT( SMLNUM ) / EPS
518: BIGNUM = ONE / SMLNUM
519: *
520: * Scale A if max element outside range [SMLNUM,BIGNUM]
521: *
522: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
523: ILASCL = .FALSE.
524: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
525: ANRMTO = SMLNUM
526: ILASCL = .TRUE.
527: ELSE IF( ANRM.GT.BIGNUM ) THEN
528: ANRMTO = BIGNUM
529: ILASCL = .TRUE.
530: END IF
531: IF( ILASCL )
532: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
533: *
534: * Scale B if max element outside range [SMLNUM,BIGNUM]
535: *
536: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
537: ILBSCL = .FALSE.
538: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
539: BNRMTO = SMLNUM
540: ILBSCL = .TRUE.
541: ELSE IF( BNRM.GT.BIGNUM ) THEN
542: BNRMTO = BIGNUM
543: ILBSCL = .TRUE.
544: END IF
545: IF( ILBSCL )
546: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
547: *
548: * Permute the matrix to make it more nearly triangular
549: * (Real Workspace: need 6*N)
550: *
551: ILEFT = 1
552: IRIGHT = N + 1
553: IRWRK = IRIGHT + N
554: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
555: $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
556: *
557: * Reduce B to triangular form (QR decomposition of B)
558: * (Complex Workspace: need N, prefer N*NB)
559: *
560: IROWS = IHI + 1 - ILO
561: ICOLS = N + 1 - ILO
562: ITAU = 1
563: IWRK = ITAU + IROWS
564: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
565: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
566: *
567: * Apply the unitary transformation to matrix A
568: * (Complex Workspace: need N, prefer N*NB)
569: *
570: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
571: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
572: $ LWORK+1-IWRK, IERR )
573: *
574: * Initialize VSL
575: * (Complex Workspace: need N, prefer N*NB)
576: *
577: IF( ILVSL ) THEN
578: CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
579: IF( IROWS.GT.1 ) THEN
580: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
581: $ VSL( ILO+1, ILO ), LDVSL )
582: END IF
583: CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
584: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
585: END IF
586: *
587: * Initialize VSR
588: *
589: IF( ILVSR )
590: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
591: *
592: * Reduce to generalized Hessenberg form
593: * (Workspace: none needed)
594: *
595: CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
596: $ LDVSL, VSR, LDVSR, IERR )
597: *
598: SDIM = 0
599: *
600: * Perform QZ algorithm, computing Schur vectors if desired
601: * (Complex Workspace: need N)
602: * (Real Workspace: need N)
603: *
604: IWRK = ITAU
605: CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
606: $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
607: $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
608: IF( IERR.NE.0 ) THEN
609: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
610: INFO = IERR
611: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
612: INFO = IERR - N
613: ELSE
614: INFO = N + 1
615: END IF
616: GO TO 40
617: END IF
618: *
619: * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
620: * condition number(s)
621: *
622: IF( WANTST ) THEN
623: *
624: * Undo scaling on eigenvalues before SELCTGing
625: *
626: IF( ILASCL )
627: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
628: IF( ILBSCL )
629: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
630: *
631: * Select eigenvalues
632: *
633: DO 10 I = 1, N
634: BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
635: 10 CONTINUE
636: *
637: * Reorder eigenvalues, transform Generalized Schur vectors, and
638: * compute reciprocal condition numbers
639: * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
640: * otherwise, need 1 )
641: *
642: CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
643: $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
644: $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
645: $ IERR )
646: *
647: IF( IJOB.GE.1 )
648: $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
649: IF( IERR.EQ.-21 ) THEN
650: *
651: * not enough complex workspace
652: *
653: INFO = -21
654: ELSE
655: IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
656: RCONDE( 1 ) = PL
657: RCONDE( 2 ) = PR
658: END IF
659: IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
660: RCONDV( 1 ) = DIF( 1 )
661: RCONDV( 2 ) = DIF( 2 )
662: END IF
663: IF( IERR.EQ.1 )
664: $ INFO = N + 3
665: END IF
666: *
667: END IF
668: *
669: * Apply permutation to VSL and VSR
670: * (Workspace: none needed)
671: *
672: IF( ILVSL )
673: $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
674: $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
675: *
676: IF( ILVSR )
677: $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
678: $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
679: *
680: * Undo scaling
681: *
682: IF( ILASCL ) THEN
683: CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
684: CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
685: END IF
686: *
687: IF( ILBSCL ) THEN
688: CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
689: CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
690: END IF
691: *
692: IF( WANTST ) THEN
693: *
694: * Check if reordering is correct
695: *
696: LASTSL = .TRUE.
697: SDIM = 0
698: DO 30 I = 1, N
699: CURSL = SELCTG( ALPHA( I ), BETA( I ) )
700: IF( CURSL )
701: $ SDIM = SDIM + 1
702: IF( CURSL .AND. .NOT.LASTSL )
703: $ INFO = N + 2
704: LASTSL = CURSL
705: 30 CONTINUE
706: *
707: END IF
708: *
709: 40 CONTINUE
710: *
711: WORK( 1 ) = MAXWRK
712: IWORK( 1 ) = LIWMIN
713: *
714: RETURN
715: *
716: * End of ZGGESX
717: *
718: 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