1 SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
2 $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBVL, JOBVR
11 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
16 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
17 $ WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
24 * (A,B), the generalized eigenvalues, and optionally, the left and/or
25 * right generalized eigenvectors.
26 *
27 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
28 * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
29 * singular. It is usually represented as the pair (alpha,beta), as
30 * there is a reasonable interpretation for beta=0, and even for both
31 * being zero.
32 *
33 * The right generalized eigenvector v(j) corresponding to the
34 * generalized eigenvalue lambda(j) of (A,B) satisfies
35 *
36 * A * v(j) = lambda(j) * B * v(j).
37 *
38 * The left generalized eigenvector u(j) corresponding to the
39 * generalized eigenvalues lambda(j) of (A,B) satisfies
40 *
41 * u(j)**H * A = lambda(j) * u(j)**H * B
42 *
43 * where u(j)**H is the conjugate-transpose of u(j).
44 *
45 * Arguments
46 * =========
47 *
48 * JOBVL (input) CHARACTER*1
49 * = 'N': do not compute the left generalized eigenvectors;
50 * = 'V': compute the left generalized eigenvectors.
51 *
52 * JOBVR (input) CHARACTER*1
53 * = 'N': do not compute the right generalized eigenvectors;
54 * = 'V': compute the right generalized eigenvectors.
55 *
56 * N (input) INTEGER
57 * The order of the matrices A, B, VL, and VR. N >= 0.
58 *
59 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
60 * On entry, the matrix A in the pair (A,B).
61 * On exit, A has been overwritten.
62 *
63 * LDA (input) INTEGER
64 * The leading dimension of A. LDA >= max(1,N).
65 *
66 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
67 * On entry, the matrix B in the pair (A,B).
68 * On exit, B has been overwritten.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of B. LDB >= max(1,N).
72 *
73 * ALPHA (output) COMPLEX*16 array, dimension (N)
74 * BETA (output) COMPLEX*16 array, dimension (N)
75 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
76 * generalized eigenvalues.
77 *
78 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
79 * underflow, and BETA(j) may even be zero. Thus, the user
80 * should avoid naively computing the ratio alpha/beta.
81 * However, ALPHA will be always less than and usually
82 * comparable with norm(A) in magnitude, and BETA always less
83 * than and usually comparable with norm(B).
84 *
85 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
86 * If JOBVL = 'V', the left generalized eigenvectors u(j) are
87 * stored one after another in the columns of VL, in the same
88 * order as their eigenvalues.
89 * Each eigenvector is scaled so the largest component has
90 * abs(real part) + abs(imag. part) = 1.
91 * Not referenced if JOBVL = 'N'.
92 *
93 * LDVL (input) INTEGER
94 * The leading dimension of the matrix VL. LDVL >= 1, and
95 * if JOBVL = 'V', LDVL >= N.
96 *
97 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
98 * If JOBVR = 'V', the right generalized eigenvectors v(j) are
99 * stored one after another in the columns of VR, in the same
100 * order as their eigenvalues.
101 * Each eigenvector is scaled so the largest component has
102 * abs(real part) + abs(imag. part) = 1.
103 * Not referenced if JOBVR = 'N'.
104 *
105 * LDVR (input) INTEGER
106 * The leading dimension of the matrix VR. LDVR >= 1, and
107 * if JOBVR = 'V', LDVR >= N.
108 *
109 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
110 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
111 *
112 * LWORK (input) INTEGER
113 * The dimension of the array WORK. LWORK >= max(1,2*N).
114 * For good performance, LWORK must generally be larger.
115 *
116 * If LWORK = -1, then a workspace query is assumed; the routine
117 * only calculates the optimal size of the WORK array, returns
118 * this value as the first entry of the WORK array, and no error
119 * message related to LWORK is issued by XERBLA.
120 *
121 * RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N)
122 *
123 * INFO (output) INTEGER
124 * = 0: successful exit
125 * < 0: if INFO = -i, the i-th argument had an illegal value.
126 * =1,...,N:
127 * The QZ iteration failed. No eigenvectors have been
128 * calculated, but ALPHA(j) and BETA(j) should be
129 * correct for j=INFO+1,...,N.
130 * > N: =N+1: other then QZ iteration failed in DHGEQZ,
131 * =N+2: error return from DTGEVC.
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136 DOUBLE PRECISION ZERO, ONE
137 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
138 COMPLEX*16 CZERO, CONE
139 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
140 $ CONE = ( 1.0D0, 0.0D0 ) )
141 * ..
142 * .. Local Scalars ..
143 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
144 CHARACTER CHTEMP
145 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
146 $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
147 $ LWKMIN, LWKOPT
148 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
149 $ SMLNUM, TEMP
150 COMPLEX*16 X
151 * ..
152 * .. Local Arrays ..
153 LOGICAL LDUMMA( 1 )
154 * ..
155 * .. External Subroutines ..
156 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
157 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
158 $ ZUNMQR
159 * ..
160 * .. External Functions ..
161 LOGICAL LSAME
162 INTEGER ILAENV
163 DOUBLE PRECISION DLAMCH, ZLANGE
164 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
165 * ..
166 * .. Intrinsic Functions ..
167 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
168 * ..
169 * .. Statement Functions ..
170 DOUBLE PRECISION ABS1
171 * ..
172 * .. Statement Function definitions ..
173 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
174 * ..
175 * .. Executable Statements ..
176 *
177 * Decode the input arguments
178 *
179 IF( LSAME( JOBVL, 'N' ) ) THEN
180 IJOBVL = 1
181 ILVL = .FALSE.
182 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
183 IJOBVL = 2
184 ILVL = .TRUE.
185 ELSE
186 IJOBVL = -1
187 ILVL = .FALSE.
188 END IF
189 *
190 IF( LSAME( JOBVR, 'N' ) ) THEN
191 IJOBVR = 1
192 ILVR = .FALSE.
193 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
194 IJOBVR = 2
195 ILVR = .TRUE.
196 ELSE
197 IJOBVR = -1
198 ILVR = .FALSE.
199 END IF
200 ILV = ILVL .OR. ILVR
201 *
202 * Test the input arguments
203 *
204 INFO = 0
205 LQUERY = ( LWORK.EQ.-1 )
206 IF( IJOBVL.LE.0 ) THEN
207 INFO = -1
208 ELSE IF( IJOBVR.LE.0 ) THEN
209 INFO = -2
210 ELSE IF( N.LT.0 ) THEN
211 INFO = -3
212 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
213 INFO = -5
214 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
215 INFO = -7
216 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
217 INFO = -11
218 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
219 INFO = -13
220 END IF
221 *
222 * Compute workspace
223 * (Note: Comments in the code beginning "Workspace:" describe the
224 * minimal amount of workspace needed at that point in the code,
225 * as well as the preferred amount for good performance.
226 * NB refers to the optimal block size for the immediately
227 * following subroutine, as returned by ILAENV. The workspace is
228 * computed assuming ILO = 1 and IHI = N, the worst case.)
229 *
230 IF( INFO.EQ.0 ) THEN
231 LWKMIN = MAX( 1, 2*N )
232 LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
233 LWKOPT = MAX( LWKOPT, N +
234 $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
235 IF( ILVL ) THEN
236 LWKOPT = MAX( LWKOPT, N +
237 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
238 END IF
239 WORK( 1 ) = LWKOPT
240 *
241 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
242 $ INFO = -15
243 END IF
244 *
245 IF( INFO.NE.0 ) THEN
246 CALL XERBLA( 'ZGGEV ', -INFO )
247 RETURN
248 ELSE IF( LQUERY ) THEN
249 RETURN
250 END IF
251 *
252 * Quick return if possible
253 *
254 IF( N.EQ.0 )
255 $ RETURN
256 *
257 * Get machine constants
258 *
259 EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
260 SMLNUM = DLAMCH( 'S' )
261 BIGNUM = ONE / SMLNUM
262 CALL DLABAD( SMLNUM, BIGNUM )
263 SMLNUM = SQRT( SMLNUM ) / EPS
264 BIGNUM = ONE / SMLNUM
265 *
266 * Scale A if max element outside range [SMLNUM,BIGNUM]
267 *
268 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
269 ILASCL = .FALSE.
270 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
271 ANRMTO = SMLNUM
272 ILASCL = .TRUE.
273 ELSE IF( ANRM.GT.BIGNUM ) THEN
274 ANRMTO = BIGNUM
275 ILASCL = .TRUE.
276 END IF
277 IF( ILASCL )
278 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
279 *
280 * Scale B if max element outside range [SMLNUM,BIGNUM]
281 *
282 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
283 ILBSCL = .FALSE.
284 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
285 BNRMTO = SMLNUM
286 ILBSCL = .TRUE.
287 ELSE IF( BNRM.GT.BIGNUM ) THEN
288 BNRMTO = BIGNUM
289 ILBSCL = .TRUE.
290 END IF
291 IF( ILBSCL )
292 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
293 *
294 * Permute the matrices A, B to isolate eigenvalues if possible
295 * (Real Workspace: need 6*N)
296 *
297 ILEFT = 1
298 IRIGHT = N + 1
299 IRWRK = IRIGHT + N
300 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
301 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
302 *
303 * Reduce B to triangular form (QR decomposition of B)
304 * (Complex Workspace: need N, prefer N*NB)
305 *
306 IROWS = IHI + 1 - ILO
307 IF( ILV ) THEN
308 ICOLS = N + 1 - ILO
309 ELSE
310 ICOLS = IROWS
311 END IF
312 ITAU = 1
313 IWRK = ITAU + IROWS
314 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
315 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
316 *
317 * Apply the orthogonal transformation to matrix A
318 * (Complex Workspace: need N, prefer N*NB)
319 *
320 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
321 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
322 $ LWORK+1-IWRK, IERR )
323 *
324 * Initialize VL
325 * (Complex Workspace: need N, prefer N*NB)
326 *
327 IF( ILVL ) THEN
328 CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
329 IF( IROWS.GT.1 ) THEN
330 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
331 $ VL( ILO+1, ILO ), LDVL )
332 END IF
333 CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
334 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
335 END IF
336 *
337 * Initialize VR
338 *
339 IF( ILVR )
340 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
341 *
342 * Reduce to generalized Hessenberg form
343 *
344 IF( ILV ) THEN
345 *
346 * Eigenvectors requested -- work on whole matrix.
347 *
348 CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
349 $ LDVL, VR, LDVR, IERR )
350 ELSE
351 CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
352 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
353 END IF
354 *
355 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
356 * Schur form and Schur vectors)
357 * (Complex Workspace: need N)
358 * (Real Workspace: need N)
359 *
360 IWRK = ITAU
361 IF( ILV ) THEN
362 CHTEMP = 'S'
363 ELSE
364 CHTEMP = 'E'
365 END IF
366 CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
367 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
368 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
369 IF( IERR.NE.0 ) THEN
370 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
371 INFO = IERR
372 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
373 INFO = IERR - N
374 ELSE
375 INFO = N + 1
376 END IF
377 GO TO 70
378 END IF
379 *
380 * Compute Eigenvectors
381 * (Real Workspace: need 2*N)
382 * (Complex Workspace: need 2*N)
383 *
384 IF( ILV ) THEN
385 IF( ILVL ) THEN
386 IF( ILVR ) THEN
387 CHTEMP = 'B'
388 ELSE
389 CHTEMP = 'L'
390 END IF
391 ELSE
392 CHTEMP = 'R'
393 END IF
394 *
395 CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
396 $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
397 $ IERR )
398 IF( IERR.NE.0 ) THEN
399 INFO = N + 2
400 GO TO 70
401 END IF
402 *
403 * Undo balancing on VL and VR and normalization
404 * (Workspace: none needed)
405 *
406 IF( ILVL ) THEN
407 CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
408 $ RWORK( IRIGHT ), N, VL, LDVL, IERR )
409 DO 30 JC = 1, N
410 TEMP = ZERO
411 DO 10 JR = 1, N
412 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
413 10 CONTINUE
414 IF( TEMP.LT.SMLNUM )
415 $ GO TO 30
416 TEMP = ONE / TEMP
417 DO 20 JR = 1, N
418 VL( JR, JC ) = VL( JR, JC )*TEMP
419 20 CONTINUE
420 30 CONTINUE
421 END IF
422 IF( ILVR ) THEN
423 CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
424 $ RWORK( IRIGHT ), N, VR, LDVR, IERR )
425 DO 60 JC = 1, N
426 TEMP = ZERO
427 DO 40 JR = 1, N
428 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
429 40 CONTINUE
430 IF( TEMP.LT.SMLNUM )
431 $ GO TO 60
432 TEMP = ONE / TEMP
433 DO 50 JR = 1, N
434 VR( JR, JC ) = VR( JR, JC )*TEMP
435 50 CONTINUE
436 60 CONTINUE
437 END IF
438 END IF
439 *
440 * Undo scaling if necessary
441 *
442 IF( ILASCL )
443 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
444 *
445 IF( ILBSCL )
446 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
447 *
448 70 CONTINUE
449 WORK( 1 ) = LWKOPT
450 *
451 RETURN
452 *
453 * End of ZGGEV
454 *
455 END
2 $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBVL, JOBVR
11 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
16 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
17 $ WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
24 * (A,B), the generalized eigenvalues, and optionally, the left and/or
25 * right generalized eigenvectors.
26 *
27 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
28 * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
29 * singular. It is usually represented as the pair (alpha,beta), as
30 * there is a reasonable interpretation for beta=0, and even for both
31 * being zero.
32 *
33 * The right generalized eigenvector v(j) corresponding to the
34 * generalized eigenvalue lambda(j) of (A,B) satisfies
35 *
36 * A * v(j) = lambda(j) * B * v(j).
37 *
38 * The left generalized eigenvector u(j) corresponding to the
39 * generalized eigenvalues lambda(j) of (A,B) satisfies
40 *
41 * u(j)**H * A = lambda(j) * u(j)**H * B
42 *
43 * where u(j)**H is the conjugate-transpose of u(j).
44 *
45 * Arguments
46 * =========
47 *
48 * JOBVL (input) CHARACTER*1
49 * = 'N': do not compute the left generalized eigenvectors;
50 * = 'V': compute the left generalized eigenvectors.
51 *
52 * JOBVR (input) CHARACTER*1
53 * = 'N': do not compute the right generalized eigenvectors;
54 * = 'V': compute the right generalized eigenvectors.
55 *
56 * N (input) INTEGER
57 * The order of the matrices A, B, VL, and VR. N >= 0.
58 *
59 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
60 * On entry, the matrix A in the pair (A,B).
61 * On exit, A has been overwritten.
62 *
63 * LDA (input) INTEGER
64 * The leading dimension of A. LDA >= max(1,N).
65 *
66 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
67 * On entry, the matrix B in the pair (A,B).
68 * On exit, B has been overwritten.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of B. LDB >= max(1,N).
72 *
73 * ALPHA (output) COMPLEX*16 array, dimension (N)
74 * BETA (output) COMPLEX*16 array, dimension (N)
75 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
76 * generalized eigenvalues.
77 *
78 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
79 * underflow, and BETA(j) may even be zero. Thus, the user
80 * should avoid naively computing the ratio alpha/beta.
81 * However, ALPHA will be always less than and usually
82 * comparable with norm(A) in magnitude, and BETA always less
83 * than and usually comparable with norm(B).
84 *
85 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
86 * If JOBVL = 'V', the left generalized eigenvectors u(j) are
87 * stored one after another in the columns of VL, in the same
88 * order as their eigenvalues.
89 * Each eigenvector is scaled so the largest component has
90 * abs(real part) + abs(imag. part) = 1.
91 * Not referenced if JOBVL = 'N'.
92 *
93 * LDVL (input) INTEGER
94 * The leading dimension of the matrix VL. LDVL >= 1, and
95 * if JOBVL = 'V', LDVL >= N.
96 *
97 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
98 * If JOBVR = 'V', the right generalized eigenvectors v(j) are
99 * stored one after another in the columns of VR, in the same
100 * order as their eigenvalues.
101 * Each eigenvector is scaled so the largest component has
102 * abs(real part) + abs(imag. part) = 1.
103 * Not referenced if JOBVR = 'N'.
104 *
105 * LDVR (input) INTEGER
106 * The leading dimension of the matrix VR. LDVR >= 1, and
107 * if JOBVR = 'V', LDVR >= N.
108 *
109 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
110 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
111 *
112 * LWORK (input) INTEGER
113 * The dimension of the array WORK. LWORK >= max(1,2*N).
114 * For good performance, LWORK must generally be larger.
115 *
116 * If LWORK = -1, then a workspace query is assumed; the routine
117 * only calculates the optimal size of the WORK array, returns
118 * this value as the first entry of the WORK array, and no error
119 * message related to LWORK is issued by XERBLA.
120 *
121 * RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N)
122 *
123 * INFO (output) INTEGER
124 * = 0: successful exit
125 * < 0: if INFO = -i, the i-th argument had an illegal value.
126 * =1,...,N:
127 * The QZ iteration failed. No eigenvectors have been
128 * calculated, but ALPHA(j) and BETA(j) should be
129 * correct for j=INFO+1,...,N.
130 * > N: =N+1: other then QZ iteration failed in DHGEQZ,
131 * =N+2: error return from DTGEVC.
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136 DOUBLE PRECISION ZERO, ONE
137 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
138 COMPLEX*16 CZERO, CONE
139 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
140 $ CONE = ( 1.0D0, 0.0D0 ) )
141 * ..
142 * .. Local Scalars ..
143 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
144 CHARACTER CHTEMP
145 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
146 $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
147 $ LWKMIN, LWKOPT
148 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
149 $ SMLNUM, TEMP
150 COMPLEX*16 X
151 * ..
152 * .. Local Arrays ..
153 LOGICAL LDUMMA( 1 )
154 * ..
155 * .. External Subroutines ..
156 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
157 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
158 $ ZUNMQR
159 * ..
160 * .. External Functions ..
161 LOGICAL LSAME
162 INTEGER ILAENV
163 DOUBLE PRECISION DLAMCH, ZLANGE
164 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
165 * ..
166 * .. Intrinsic Functions ..
167 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
168 * ..
169 * .. Statement Functions ..
170 DOUBLE PRECISION ABS1
171 * ..
172 * .. Statement Function definitions ..
173 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
174 * ..
175 * .. Executable Statements ..
176 *
177 * Decode the input arguments
178 *
179 IF( LSAME( JOBVL, 'N' ) ) THEN
180 IJOBVL = 1
181 ILVL = .FALSE.
182 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
183 IJOBVL = 2
184 ILVL = .TRUE.
185 ELSE
186 IJOBVL = -1
187 ILVL = .FALSE.
188 END IF
189 *
190 IF( LSAME( JOBVR, 'N' ) ) THEN
191 IJOBVR = 1
192 ILVR = .FALSE.
193 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
194 IJOBVR = 2
195 ILVR = .TRUE.
196 ELSE
197 IJOBVR = -1
198 ILVR = .FALSE.
199 END IF
200 ILV = ILVL .OR. ILVR
201 *
202 * Test the input arguments
203 *
204 INFO = 0
205 LQUERY = ( LWORK.EQ.-1 )
206 IF( IJOBVL.LE.0 ) THEN
207 INFO = -1
208 ELSE IF( IJOBVR.LE.0 ) THEN
209 INFO = -2
210 ELSE IF( N.LT.0 ) THEN
211 INFO = -3
212 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
213 INFO = -5
214 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
215 INFO = -7
216 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
217 INFO = -11
218 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
219 INFO = -13
220 END IF
221 *
222 * Compute workspace
223 * (Note: Comments in the code beginning "Workspace:" describe the
224 * minimal amount of workspace needed at that point in the code,
225 * as well as the preferred amount for good performance.
226 * NB refers to the optimal block size for the immediately
227 * following subroutine, as returned by ILAENV. The workspace is
228 * computed assuming ILO = 1 and IHI = N, the worst case.)
229 *
230 IF( INFO.EQ.0 ) THEN
231 LWKMIN = MAX( 1, 2*N )
232 LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
233 LWKOPT = MAX( LWKOPT, N +
234 $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
235 IF( ILVL ) THEN
236 LWKOPT = MAX( LWKOPT, N +
237 $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
238 END IF
239 WORK( 1 ) = LWKOPT
240 *
241 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
242 $ INFO = -15
243 END IF
244 *
245 IF( INFO.NE.0 ) THEN
246 CALL XERBLA( 'ZGGEV ', -INFO )
247 RETURN
248 ELSE IF( LQUERY ) THEN
249 RETURN
250 END IF
251 *
252 * Quick return if possible
253 *
254 IF( N.EQ.0 )
255 $ RETURN
256 *
257 * Get machine constants
258 *
259 EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
260 SMLNUM = DLAMCH( 'S' )
261 BIGNUM = ONE / SMLNUM
262 CALL DLABAD( SMLNUM, BIGNUM )
263 SMLNUM = SQRT( SMLNUM ) / EPS
264 BIGNUM = ONE / SMLNUM
265 *
266 * Scale A if max element outside range [SMLNUM,BIGNUM]
267 *
268 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
269 ILASCL = .FALSE.
270 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
271 ANRMTO = SMLNUM
272 ILASCL = .TRUE.
273 ELSE IF( ANRM.GT.BIGNUM ) THEN
274 ANRMTO = BIGNUM
275 ILASCL = .TRUE.
276 END IF
277 IF( ILASCL )
278 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
279 *
280 * Scale B if max element outside range [SMLNUM,BIGNUM]
281 *
282 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
283 ILBSCL = .FALSE.
284 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
285 BNRMTO = SMLNUM
286 ILBSCL = .TRUE.
287 ELSE IF( BNRM.GT.BIGNUM ) THEN
288 BNRMTO = BIGNUM
289 ILBSCL = .TRUE.
290 END IF
291 IF( ILBSCL )
292 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
293 *
294 * Permute the matrices A, B to isolate eigenvalues if possible
295 * (Real Workspace: need 6*N)
296 *
297 ILEFT = 1
298 IRIGHT = N + 1
299 IRWRK = IRIGHT + N
300 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
301 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
302 *
303 * Reduce B to triangular form (QR decomposition of B)
304 * (Complex Workspace: need N, prefer N*NB)
305 *
306 IROWS = IHI + 1 - ILO
307 IF( ILV ) THEN
308 ICOLS = N + 1 - ILO
309 ELSE
310 ICOLS = IROWS
311 END IF
312 ITAU = 1
313 IWRK = ITAU + IROWS
314 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
315 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
316 *
317 * Apply the orthogonal transformation to matrix A
318 * (Complex Workspace: need N, prefer N*NB)
319 *
320 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
321 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
322 $ LWORK+1-IWRK, IERR )
323 *
324 * Initialize VL
325 * (Complex Workspace: need N, prefer N*NB)
326 *
327 IF( ILVL ) THEN
328 CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
329 IF( IROWS.GT.1 ) THEN
330 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
331 $ VL( ILO+1, ILO ), LDVL )
332 END IF
333 CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
334 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
335 END IF
336 *
337 * Initialize VR
338 *
339 IF( ILVR )
340 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
341 *
342 * Reduce to generalized Hessenberg form
343 *
344 IF( ILV ) THEN
345 *
346 * Eigenvectors requested -- work on whole matrix.
347 *
348 CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
349 $ LDVL, VR, LDVR, IERR )
350 ELSE
351 CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
352 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
353 END IF
354 *
355 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
356 * Schur form and Schur vectors)
357 * (Complex Workspace: need N)
358 * (Real Workspace: need N)
359 *
360 IWRK = ITAU
361 IF( ILV ) THEN
362 CHTEMP = 'S'
363 ELSE
364 CHTEMP = 'E'
365 END IF
366 CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
367 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
368 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
369 IF( IERR.NE.0 ) THEN
370 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
371 INFO = IERR
372 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
373 INFO = IERR - N
374 ELSE
375 INFO = N + 1
376 END IF
377 GO TO 70
378 END IF
379 *
380 * Compute Eigenvectors
381 * (Real Workspace: need 2*N)
382 * (Complex Workspace: need 2*N)
383 *
384 IF( ILV ) THEN
385 IF( ILVL ) THEN
386 IF( ILVR ) THEN
387 CHTEMP = 'B'
388 ELSE
389 CHTEMP = 'L'
390 END IF
391 ELSE
392 CHTEMP = 'R'
393 END IF
394 *
395 CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
396 $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
397 $ IERR )
398 IF( IERR.NE.0 ) THEN
399 INFO = N + 2
400 GO TO 70
401 END IF
402 *
403 * Undo balancing on VL and VR and normalization
404 * (Workspace: none needed)
405 *
406 IF( ILVL ) THEN
407 CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
408 $ RWORK( IRIGHT ), N, VL, LDVL, IERR )
409 DO 30 JC = 1, N
410 TEMP = ZERO
411 DO 10 JR = 1, N
412 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
413 10 CONTINUE
414 IF( TEMP.LT.SMLNUM )
415 $ GO TO 30
416 TEMP = ONE / TEMP
417 DO 20 JR = 1, N
418 VL( JR, JC ) = VL( JR, JC )*TEMP
419 20 CONTINUE
420 30 CONTINUE
421 END IF
422 IF( ILVR ) THEN
423 CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
424 $ RWORK( IRIGHT ), N, VR, LDVR, IERR )
425 DO 60 JC = 1, N
426 TEMP = ZERO
427 DO 40 JR = 1, N
428 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
429 40 CONTINUE
430 IF( TEMP.LT.SMLNUM )
431 $ GO TO 60
432 TEMP = ONE / TEMP
433 DO 50 JR = 1, N
434 VR( JR, JC ) = VR( JR, JC )*TEMP
435 50 CONTINUE
436 60 CONTINUE
437 END IF
438 END IF
439 *
440 * Undo scaling if necessary
441 *
442 IF( ILASCL )
443 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
444 *
445 IF( ILBSCL )
446 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
447 *
448 70 CONTINUE
449 WORK( 1 ) = LWKOPT
450 *
451 RETURN
452 *
453 * End of ZGGEV
454 *
455 END