1: *> \brief <b> ZGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors 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 ZGGEVX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggevx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggevx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggevx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
22: * ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
23: * LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
24: * WORK, LWORK, RWORK, IWORK, BWORK, INFO )
25: *
26: * .. Scalar Arguments ..
27: * CHARACTER BALANC, JOBVL, JOBVR, SENSE
28: * INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
29: * DOUBLE PRECISION ABNRM, BBNRM
30: * ..
31: * .. Array Arguments ..
32: * LOGICAL BWORK( * )
33: * INTEGER IWORK( * )
34: * DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
35: * $ RSCALE( * ), RWORK( * )
36: * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
37: * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
38: * $ WORK( * )
39: * ..
40: *
41: *
42: *> \par Purpose:
43: * =============
44: *>
45: *> \verbatim
46: *>
47: *> ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
48: *> (A,B) the generalized eigenvalues, and optionally, the left and/or
49: *> right generalized eigenvectors.
50: *>
51: *> Optionally, it also computes a balancing transformation to improve
52: *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
53: *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
54: *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
55: *> right eigenvectors (RCONDV).
56: *>
57: *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
58: *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
59: *> singular. It is usually represented as the pair (alpha,beta), as
60: *> there is a reasonable interpretation for beta=0, and even for both
61: *> being zero.
62: *>
63: *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
64: *> of (A,B) satisfies
65: *> A * v(j) = lambda(j) * B * v(j) .
66: *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
67: *> of (A,B) satisfies
68: *> u(j)**H * A = lambda(j) * u(j)**H * B.
69: *> where u(j)**H is the conjugate-transpose of u(j).
70: *>
71: *> \endverbatim
72: *
73: * Arguments:
74: * ==========
75: *
76: *> \param[in] BALANC
77: *> \verbatim
78: *> BALANC is CHARACTER*1
79: *> Specifies the balance option to be performed:
80: *> = 'N': do not diagonally scale or permute;
81: *> = 'P': permute only;
82: *> = 'S': scale only;
83: *> = 'B': both permute and scale.
84: *> Computed reciprocal condition numbers will be for the
85: *> matrices after permuting and/or balancing. Permuting does
86: *> not change condition numbers (in exact arithmetic), but
87: *> balancing does.
88: *> \endverbatim
89: *>
90: *> \param[in] JOBVL
91: *> \verbatim
92: *> JOBVL is CHARACTER*1
93: *> = 'N': do not compute the left generalized eigenvectors;
94: *> = 'V': compute the left generalized eigenvectors.
95: *> \endverbatim
96: *>
97: *> \param[in] JOBVR
98: *> \verbatim
99: *> JOBVR is CHARACTER*1
100: *> = 'N': do not compute the right generalized eigenvectors;
101: *> = 'V': compute the right generalized eigenvectors.
102: *> \endverbatim
103: *>
104: *> \param[in] SENSE
105: *> \verbatim
106: *> SENSE is CHARACTER*1
107: *> Determines which reciprocal condition numbers are computed.
108: *> = 'N': none are computed;
109: *> = 'E': computed for eigenvalues only;
110: *> = 'V': computed for eigenvectors only;
111: *> = 'B': computed for eigenvalues and eigenvectors.
112: *> \endverbatim
113: *>
114: *> \param[in] N
115: *> \verbatim
116: *> N is INTEGER
117: *> The order of the matrices A, B, VL, and VR. N >= 0.
118: *> \endverbatim
119: *>
120: *> \param[in,out] A
121: *> \verbatim
122: *> A is COMPLEX*16 array, dimension (LDA, N)
123: *> On entry, the matrix A in the pair (A,B).
124: *> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
125: *> or both, then A contains the first part of the complex Schur
126: *> form of the "balanced" versions of the input A and B.
127: *> \endverbatim
128: *>
129: *> \param[in] LDA
130: *> \verbatim
131: *> LDA is INTEGER
132: *> The leading dimension of A. LDA >= max(1,N).
133: *> \endverbatim
134: *>
135: *> \param[in,out] B
136: *> \verbatim
137: *> B is COMPLEX*16 array, dimension (LDB, N)
138: *> On entry, the matrix B in the pair (A,B).
139: *> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
140: *> or both, then B contains the second part of the complex
141: *> Schur form of the "balanced" versions of the input A and B.
142: *> \endverbatim
143: *>
144: *> \param[in] LDB
145: *> \verbatim
146: *> LDB is INTEGER
147: *> The leading dimension of B. LDB >= max(1,N).
148: *> \endverbatim
149: *>
150: *> \param[out] ALPHA
151: *> \verbatim
152: *> ALPHA is COMPLEX*16 array, dimension (N)
153: *> \endverbatim
154: *>
155: *> \param[out] BETA
156: *> \verbatim
157: *> BETA is COMPLEX*16 array, dimension (N)
158: *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
159: *> eigenvalues.
160: *>
161: *> Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
162: *> underflow, and BETA(j) may even be zero. Thus, the user
163: *> should avoid naively computing the ratio ALPHA/BETA.
164: *> However, ALPHA will be always less than and usually
165: *> comparable with norm(A) in magnitude, and BETA always less
166: *> than and usually comparable with norm(B).
167: *> \endverbatim
168: *>
169: *> \param[out] VL
170: *> \verbatim
171: *> VL is COMPLEX*16 array, dimension (LDVL,N)
172: *> If JOBVL = 'V', the left generalized eigenvectors u(j) are
173: *> stored one after another in the columns of VL, in the same
174: *> order as their eigenvalues.
175: *> Each eigenvector will be scaled so the largest component
176: *> will have abs(real part) + abs(imag. part) = 1.
177: *> Not referenced if JOBVL = 'N'.
178: *> \endverbatim
179: *>
180: *> \param[in] LDVL
181: *> \verbatim
182: *> LDVL is INTEGER
183: *> The leading dimension of the matrix VL. LDVL >= 1, and
184: *> if JOBVL = 'V', LDVL >= N.
185: *> \endverbatim
186: *>
187: *> \param[out] VR
188: *> \verbatim
189: *> VR is COMPLEX*16 array, dimension (LDVR,N)
190: *> If JOBVR = 'V', the right generalized eigenvectors v(j) are
191: *> stored one after another in the columns of VR, in the same
192: *> order as their eigenvalues.
193: *> Each eigenvector will be scaled so the largest component
194: *> will have abs(real part) + abs(imag. part) = 1.
195: *> Not referenced if JOBVR = 'N'.
196: *> \endverbatim
197: *>
198: *> \param[in] LDVR
199: *> \verbatim
200: *> LDVR is INTEGER
201: *> The leading dimension of the matrix VR. LDVR >= 1, and
202: *> if JOBVR = 'V', LDVR >= N.
203: *> \endverbatim
204: *>
205: *> \param[out] ILO
206: *> \verbatim
207: *> ILO is INTEGER
208: *> \endverbatim
209: *>
210: *> \param[out] IHI
211: *> \verbatim
212: *> IHI is INTEGER
213: *> ILO and IHI are integer values such that on exit
214: *> A(i,j) = 0 and B(i,j) = 0 if i > j and
215: *> j = 1,...,ILO-1 or i = IHI+1,...,N.
216: *> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
217: *> \endverbatim
218: *>
219: *> \param[out] LSCALE
220: *> \verbatim
221: *> LSCALE is DOUBLE PRECISION array, dimension (N)
222: *> Details of the permutations and scaling factors applied
223: *> to the left side of A and B. If PL(j) is the index of the
224: *> row interchanged with row j, and DL(j) is the scaling
225: *> factor applied to row j, then
226: *> LSCALE(j) = PL(j) for j = 1,...,ILO-1
227: *> = DL(j) for j = ILO,...,IHI
228: *> = PL(j) for j = IHI+1,...,N.
229: *> The order in which the interchanges are made is N to IHI+1,
230: *> then 1 to ILO-1.
231: *> \endverbatim
232: *>
233: *> \param[out] RSCALE
234: *> \verbatim
235: *> RSCALE is DOUBLE PRECISION array, dimension (N)
236: *> Details of the permutations and scaling factors applied
237: *> to the right side of A and B. If PR(j) is the index of the
238: *> column interchanged with column j, and DR(j) is the scaling
239: *> factor applied to column j, then
240: *> RSCALE(j) = PR(j) for j = 1,...,ILO-1
241: *> = DR(j) for j = ILO,...,IHI
242: *> = PR(j) for j = IHI+1,...,N
243: *> The order in which the interchanges are made is N to IHI+1,
244: *> then 1 to ILO-1.
245: *> \endverbatim
246: *>
247: *> \param[out] ABNRM
248: *> \verbatim
249: *> ABNRM is DOUBLE PRECISION
250: *> The one-norm of the balanced matrix A.
251: *> \endverbatim
252: *>
253: *> \param[out] BBNRM
254: *> \verbatim
255: *> BBNRM is DOUBLE PRECISION
256: *> The one-norm of the balanced matrix B.
257: *> \endverbatim
258: *>
259: *> \param[out] RCONDE
260: *> \verbatim
261: *> RCONDE is DOUBLE PRECISION array, dimension (N)
262: *> If SENSE = 'E' or 'B', the reciprocal condition numbers of
263: *> the eigenvalues, stored in consecutive elements of the array.
264: *> If SENSE = 'N' or 'V', RCONDE is not referenced.
265: *> \endverbatim
266: *>
267: *> \param[out] RCONDV
268: *> \verbatim
269: *> RCONDV is DOUBLE PRECISION array, dimension (N)
270: *> If JOB = 'V' or 'B', the estimated reciprocal condition
271: *> numbers of the eigenvectors, stored in consecutive elements
272: *> of the array. If the eigenvalues cannot be reordered to
273: *> compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
274: *> when the true value would be very small anyway.
275: *> If SENSE = 'N' or 'E', RCONDV is not referenced.
276: *> \endverbatim
277: *>
278: *> \param[out] WORK
279: *> \verbatim
280: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
281: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
282: *> \endverbatim
283: *>
284: *> \param[in] LWORK
285: *> \verbatim
286: *> LWORK is INTEGER
287: *> The dimension of the array WORK. LWORK >= max(1,2*N).
288: *> If SENSE = 'E', LWORK >= max(1,4*N).
289: *> If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
290: *>
291: *> If LWORK = -1, then a workspace query is assumed; the routine
292: *> only calculates the optimal size of the WORK array, returns
293: *> this value as the first entry of the WORK array, and no error
294: *> message related to LWORK is issued by XERBLA.
295: *> \endverbatim
296: *>
297: *> \param[out] RWORK
298: *> \verbatim
299: *> RWORK is DOUBLE PRECISION array, dimension (lrwork)
300: *> lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
301: *> and at least max(1,2*N) otherwise.
302: *> Real workspace.
303: *> \endverbatim
304: *>
305: *> \param[out] IWORK
306: *> \verbatim
307: *> IWORK is INTEGER array, dimension (N+2)
308: *> If SENSE = 'E', IWORK is not referenced.
309: *> \endverbatim
310: *>
311: *> \param[out] BWORK
312: *> \verbatim
313: *> BWORK is LOGICAL array, dimension (N)
314: *> If SENSE = 'N', BWORK is not referenced.
315: *> \endverbatim
316: *>
317: *> \param[out] INFO
318: *> \verbatim
319: *> INFO is INTEGER
320: *> = 0: successful exit
321: *> < 0: if INFO = -i, the i-th argument had an illegal value.
322: *> = 1,...,N:
323: *> The QZ iteration failed. No eigenvectors have been
324: *> calculated, but ALPHA(j) and BETA(j) should be correct
325: *> for j=INFO+1,...,N.
326: *> > N: =N+1: other than QZ iteration failed in ZHGEQZ.
327: *> =N+2: error return from ZTGEVC.
328: *> \endverbatim
329: *
330: * Authors:
331: * ========
332: *
333: *> \author Univ. of Tennessee
334: *> \author Univ. of California Berkeley
335: *> \author Univ. of Colorado Denver
336: *> \author NAG Ltd.
337: *
338: *> \ingroup complex16GEeigen
339: *
340: *> \par Further Details:
341: * =====================
342: *>
343: *> \verbatim
344: *>
345: *> Balancing a matrix pair (A,B) includes, first, permuting rows and
346: *> columns to isolate eigenvalues, second, applying diagonal similarity
347: *> transformation to the rows and columns to make the rows and columns
348: *> as close in norm as possible. The computed reciprocal condition
349: *> numbers correspond to the balanced matrix. Permuting rows and columns
350: *> will not change the condition numbers (in exact arithmetic) but
351: *> diagonal scaling will. For further explanation of balancing, see
352: *> section 4.11.1.2 of LAPACK Users' Guide.
353: *>
354: *> An approximate error bound on the chordal distance between the i-th
355: *> computed generalized eigenvalue w and the corresponding exact
356: *> eigenvalue lambda is
357: *>
358: *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
359: *>
360: *> An approximate error bound for the angle between the i-th computed
361: *> eigenvector VL(i) or VR(i) is given by
362: *>
363: *> EPS * norm(ABNRM, BBNRM) / DIF(i).
364: *>
365: *> For further explanation of the reciprocal condition numbers RCONDE
366: *> and RCONDV, see section 4.11 of LAPACK User's Guide.
367: *> \endverbatim
368: *>
369: * =====================================================================
370: SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
371: $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
372: $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
373: $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
374: *
375: * -- LAPACK driver routine --
376: * -- LAPACK is a software package provided by Univ. of Tennessee, --
377: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
378: *
379: * .. Scalar Arguments ..
380: CHARACTER BALANC, JOBVL, JOBVR, SENSE
381: INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
382: DOUBLE PRECISION ABNRM, BBNRM
383: * ..
384: * .. Array Arguments ..
385: LOGICAL BWORK( * )
386: INTEGER IWORK( * )
387: DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
388: $ RSCALE( * ), RWORK( * )
389: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
390: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
391: $ WORK( * )
392: * ..
393: *
394: * =====================================================================
395: *
396: * .. Parameters ..
397: DOUBLE PRECISION ZERO, ONE
398: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
399: COMPLEX*16 CZERO, CONE
400: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
401: $ CONE = ( 1.0D+0, 0.0D+0 ) )
402: * ..
403: * .. Local Scalars ..
404: LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
405: $ WANTSB, WANTSE, WANTSN, WANTSV
406: CHARACTER CHTEMP
407: INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
408: $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
409: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
410: $ SMLNUM, TEMP
411: COMPLEX*16 X
412: * ..
413: * .. Local Arrays ..
414: LOGICAL LDUMMA( 1 )
415: * ..
416: * .. External Subroutines ..
417: EXTERNAL DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
418: $ ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
419: $ ZTGSNA, ZUNGQR, ZUNMQR
420: * ..
421: * .. External Functions ..
422: LOGICAL LSAME
423: INTEGER ILAENV
424: DOUBLE PRECISION DLAMCH, ZLANGE
425: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
426: * ..
427: * .. Intrinsic Functions ..
428: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
429: * ..
430: * .. Statement Functions ..
431: DOUBLE PRECISION ABS1
432: * ..
433: * .. Statement Function definitions ..
434: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
435: * ..
436: * .. Executable Statements ..
437: *
438: * Decode the input arguments
439: *
440: IF( LSAME( JOBVL, 'N' ) ) THEN
441: IJOBVL = 1
442: ILVL = .FALSE.
443: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
444: IJOBVL = 2
445: ILVL = .TRUE.
446: ELSE
447: IJOBVL = -1
448: ILVL = .FALSE.
449: END IF
450: *
451: IF( LSAME( JOBVR, 'N' ) ) THEN
452: IJOBVR = 1
453: ILVR = .FALSE.
454: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
455: IJOBVR = 2
456: ILVR = .TRUE.
457: ELSE
458: IJOBVR = -1
459: ILVR = .FALSE.
460: END IF
461: ILV = ILVL .OR. ILVR
462: *
463: NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
464: WANTSN = LSAME( SENSE, 'N' )
465: WANTSE = LSAME( SENSE, 'E' )
466: WANTSV = LSAME( SENSE, 'V' )
467: WANTSB = LSAME( SENSE, 'B' )
468: *
469: * Test the input arguments
470: *
471: INFO = 0
472: LQUERY = ( LWORK.EQ.-1 )
473: IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
474: $ LSAME( BALANC, 'B' ) ) ) THEN
475: INFO = -1
476: ELSE IF( IJOBVL.LE.0 ) THEN
477: INFO = -2
478: ELSE IF( IJOBVR.LE.0 ) THEN
479: INFO = -3
480: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
481: $ THEN
482: INFO = -4
483: ELSE IF( N.LT.0 ) THEN
484: INFO = -5
485: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
486: INFO = -7
487: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
488: INFO = -9
489: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
490: INFO = -13
491: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
492: INFO = -15
493: END IF
494: *
495: * Compute workspace
496: * (Note: Comments in the code beginning "Workspace:" describe the
497: * minimal amount of workspace needed at that point in the code,
498: * as well as the preferred amount for good performance.
499: * NB refers to the optimal block size for the immediately
500: * following subroutine, as returned by ILAENV. The workspace is
501: * computed assuming ILO = 1 and IHI = N, the worst case.)
502: *
503: IF( INFO.EQ.0 ) THEN
504: IF( N.EQ.0 ) THEN
505: MINWRK = 1
506: MAXWRK = 1
507: ELSE
508: MINWRK = 2*N
509: IF( WANTSE ) THEN
510: MINWRK = 4*N
511: ELSE IF( WANTSV .OR. WANTSB ) THEN
512: MINWRK = 2*N*( N + 1)
513: END IF
514: MAXWRK = MINWRK
515: MAXWRK = MAX( MAXWRK,
516: $ N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
517: MAXWRK = MAX( MAXWRK,
518: $ N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
519: IF( ILVL ) THEN
520: MAXWRK = MAX( MAXWRK, N +
521: $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
522: END IF
523: END IF
524: WORK( 1 ) = MAXWRK
525: *
526: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
527: INFO = -25
528: END IF
529: END IF
530: *
531: IF( INFO.NE.0 ) THEN
532: CALL XERBLA( 'ZGGEVX', -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 )
541: $ RETURN
542: *
543: * Get machine constants
544: *
545: EPS = DLAMCH( 'P' )
546: SMLNUM = DLAMCH( 'S' )
547: BIGNUM = ONE / SMLNUM
548: CALL DLABAD( SMLNUM, BIGNUM )
549: SMLNUM = SQRT( SMLNUM ) / EPS
550: BIGNUM = ONE / SMLNUM
551: *
552: * Scale A if max element outside range [SMLNUM,BIGNUM]
553: *
554: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
555: ILASCL = .FALSE.
556: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
557: ANRMTO = SMLNUM
558: ILASCL = .TRUE.
559: ELSE IF( ANRM.GT.BIGNUM ) THEN
560: ANRMTO = BIGNUM
561: ILASCL = .TRUE.
562: END IF
563: IF( ILASCL )
564: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
565: *
566: * Scale B if max element outside range [SMLNUM,BIGNUM]
567: *
568: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
569: ILBSCL = .FALSE.
570: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
571: BNRMTO = SMLNUM
572: ILBSCL = .TRUE.
573: ELSE IF( BNRM.GT.BIGNUM ) THEN
574: BNRMTO = BIGNUM
575: ILBSCL = .TRUE.
576: END IF
577: IF( ILBSCL )
578: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
579: *
580: * Permute and/or balance the matrix pair (A,B)
581: * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
582: *
583: CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
584: $ RWORK, IERR )
585: *
586: * Compute ABNRM and BBNRM
587: *
588: ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
589: IF( ILASCL ) THEN
590: RWORK( 1 ) = ABNRM
591: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
592: $ IERR )
593: ABNRM = RWORK( 1 )
594: END IF
595: *
596: BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
597: IF( ILBSCL ) THEN
598: RWORK( 1 ) = BBNRM
599: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
600: $ IERR )
601: BBNRM = RWORK( 1 )
602: END IF
603: *
604: * Reduce B to triangular form (QR decomposition of B)
605: * (Complex Workspace: need N, prefer N*NB )
606: *
607: IROWS = IHI + 1 - ILO
608: IF( ILV .OR. .NOT.WANTSN ) THEN
609: ICOLS = N + 1 - ILO
610: ELSE
611: ICOLS = IROWS
612: END IF
613: ITAU = 1
614: IWRK = ITAU + IROWS
615: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
616: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
617: *
618: * Apply the unitary transformation to A
619: * (Complex Workspace: need N, prefer N*NB)
620: *
621: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
622: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
623: $ LWORK+1-IWRK, IERR )
624: *
625: * Initialize VL and/or VR
626: * (Workspace: need N, prefer N*NB)
627: *
628: IF( ILVL ) THEN
629: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
630: IF( IROWS.GT.1 ) THEN
631: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
632: $ VL( ILO+1, ILO ), LDVL )
633: END IF
634: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
635: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
636: END IF
637: *
638: IF( ILVR )
639: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
640: *
641: * Reduce to generalized Hessenberg form
642: * (Workspace: none needed)
643: *
644: IF( ILV .OR. .NOT.WANTSN ) THEN
645: *
646: * Eigenvectors requested -- work on whole matrix.
647: *
648: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
649: $ LDVL, VR, LDVR, IERR )
650: ELSE
651: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
652: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
653: END IF
654: *
655: * Perform QZ algorithm (Compute eigenvalues, and optionally, the
656: * Schur forms and Schur vectors)
657: * (Complex Workspace: need N)
658: * (Real Workspace: need N)
659: *
660: IWRK = ITAU
661: IF( ILV .OR. .NOT.WANTSN ) THEN
662: CHTEMP = 'S'
663: ELSE
664: CHTEMP = 'E'
665: END IF
666: *
667: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
668: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
669: $ LWORK+1-IWRK, RWORK, IERR )
670: IF( IERR.NE.0 ) THEN
671: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
672: INFO = IERR
673: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
674: INFO = IERR - N
675: ELSE
676: INFO = N + 1
677: END IF
678: GO TO 90
679: END IF
680: *
681: * Compute Eigenvectors and estimate condition numbers if desired
682: * ZTGEVC: (Complex Workspace: need 2*N )
683: * (Real Workspace: need 2*N )
684: * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
685: * (Integer Workspace: need N+2 )
686: *
687: IF( ILV .OR. .NOT.WANTSN ) THEN
688: IF( ILV ) THEN
689: IF( ILVL ) THEN
690: IF( ILVR ) THEN
691: CHTEMP = 'B'
692: ELSE
693: CHTEMP = 'L'
694: END IF
695: ELSE
696: CHTEMP = 'R'
697: END IF
698: *
699: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
700: $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
701: $ IERR )
702: IF( IERR.NE.0 ) THEN
703: INFO = N + 2
704: GO TO 90
705: END IF
706: END IF
707: *
708: IF( .NOT.WANTSN ) THEN
709: *
710: * compute eigenvectors (ZTGEVC) and estimate condition
711: * numbers (ZTGSNA). Note that the definition of the condition
712: * number is not invariant under transformation (u,v) to
713: * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
714: * Schur form (S,T), Q and Z are orthogonal matrices. In order
715: * to avoid using extra 2*N*N workspace, we have to
716: * re-calculate eigenvectors and estimate the condition numbers
717: * one at a time.
718: *
719: DO 20 I = 1, N
720: *
721: DO 10 J = 1, N
722: BWORK( J ) = .FALSE.
723: 10 CONTINUE
724: BWORK( I ) = .TRUE.
725: *
726: IWRK = N + 1
727: IWRK1 = IWRK + N
728: *
729: IF( WANTSE .OR. WANTSB ) THEN
730: CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
731: $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
732: $ WORK( IWRK1 ), RWORK, IERR )
733: IF( IERR.NE.0 ) THEN
734: INFO = N + 2
735: GO TO 90
736: END IF
737: END IF
738: *
739: CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
740: $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
741: $ RCONDV( I ), 1, M, WORK( IWRK1 ),
742: $ LWORK-IWRK1+1, IWORK, IERR )
743: *
744: 20 CONTINUE
745: END IF
746: END IF
747: *
748: * Undo balancing on VL and VR and normalization
749: * (Workspace: none needed)
750: *
751: IF( ILVL ) THEN
752: CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
753: $ LDVL, IERR )
754: *
755: DO 50 JC = 1, N
756: TEMP = ZERO
757: DO 30 JR = 1, N
758: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
759: 30 CONTINUE
760: IF( TEMP.LT.SMLNUM )
761: $ GO TO 50
762: TEMP = ONE / TEMP
763: DO 40 JR = 1, N
764: VL( JR, JC ) = VL( JR, JC )*TEMP
765: 40 CONTINUE
766: 50 CONTINUE
767: END IF
768: *
769: IF( ILVR ) THEN
770: CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
771: $ LDVR, IERR )
772: DO 80 JC = 1, N
773: TEMP = ZERO
774: DO 60 JR = 1, N
775: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
776: 60 CONTINUE
777: IF( TEMP.LT.SMLNUM )
778: $ GO TO 80
779: TEMP = ONE / TEMP
780: DO 70 JR = 1, N
781: VR( JR, JC ) = VR( JR, JC )*TEMP
782: 70 CONTINUE
783: 80 CONTINUE
784: END IF
785: *
786: * Undo scaling if necessary
787: *
788: 90 CONTINUE
789: *
790: IF( ILASCL )
791: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
792: *
793: IF( ILBSCL )
794: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
795: *
796: WORK( 1 ) = MAXWRK
797: RETURN
798: *
799: * End of ZGGEVX
800: *
801: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggevx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack