1 SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
2 $ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
3 $ RCONDV, WORK, LWORK, RWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER BALANC, JOBVL, JOBVR, SENSE
12 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
13 DOUBLE PRECISION ABNRM
14 * ..
15 * .. Array Arguments ..
16 DOUBLE PRECISION RCONDE( * ), RCONDV( * ), RWORK( * ),
17 $ SCALE( * )
18 COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ W( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
26 * eigenvalues and, optionally, the left and/or right eigenvectors.
27 *
28 * Optionally also, it computes a balancing transformation to improve
29 * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
30 * SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
31 * (RCONDE), and reciprocal condition numbers for the right
32 * eigenvectors (RCONDV).
33 *
34 * The right eigenvector v(j) of A satisfies
35 * A * v(j) = lambda(j) * v(j)
36 * where lambda(j) is its eigenvalue.
37 * The left eigenvector u(j) of A satisfies
38 * u(j)**H * A = lambda(j) * u(j)**H
39 * where u(j)**H denotes the conjugate transpose of u(j).
40 *
41 * The computed eigenvectors are normalized to have Euclidean norm
42 * equal to 1 and largest component real.
43 *
44 * Balancing a matrix means permuting the rows and columns to make it
45 * more nearly upper triangular, and applying a diagonal similarity
46 * transformation D * A * D**(-1), where D is a diagonal matrix, to
47 * make its rows and columns closer in norm and the condition numbers
48 * of its eigenvalues and eigenvectors smaller. The computed
49 * reciprocal condition numbers correspond to the balanced matrix.
50 * Permuting rows and columns will not change the condition numbers
51 * (in exact arithmetic) but diagonal scaling will. For further
52 * explanation of balancing, see section 4.10.2 of the LAPACK
53 * Users' Guide.
54 *
55 * Arguments
56 * =========
57 *
58 * BALANC (input) CHARACTER*1
59 * Indicates how the input matrix should be diagonally scaled
60 * and/or permuted to improve the conditioning of its
61 * eigenvalues.
62 * = 'N': Do not diagonally scale or permute;
63 * = 'P': Perform permutations to make the matrix more nearly
64 * upper triangular. Do not diagonally scale;
65 * = 'S': Diagonally scale the matrix, ie. replace A by
66 * D*A*D**(-1), where D is a diagonal matrix chosen
67 * to make the rows and columns of A more equal in
68 * norm. Do not permute;
69 * = 'B': Both diagonally scale and permute A.
70 *
71 * Computed reciprocal condition numbers will be for the matrix
72 * after balancing and/or permuting. Permuting does not change
73 * condition numbers (in exact arithmetic), but balancing does.
74 *
75 * JOBVL (input) CHARACTER*1
76 * = 'N': left eigenvectors of A are not computed;
77 * = 'V': left eigenvectors of A are computed.
78 * If SENSE = 'E' or 'B', JOBVL must = 'V'.
79 *
80 * JOBVR (input) CHARACTER*1
81 * = 'N': right eigenvectors of A are not computed;
82 * = 'V': right eigenvectors of A are computed.
83 * If SENSE = 'E' or 'B', JOBVR must = 'V'.
84 *
85 * SENSE (input) CHARACTER*1
86 * Determines which reciprocal condition numbers are computed.
87 * = 'N': None are computed;
88 * = 'E': Computed for eigenvalues only;
89 * = 'V': Computed for right eigenvectors only;
90 * = 'B': Computed for eigenvalues and right eigenvectors.
91 *
92 * If SENSE = 'E' or 'B', both left and right eigenvectors
93 * must also be computed (JOBVL = 'V' and JOBVR = 'V').
94 *
95 * N (input) INTEGER
96 * The order of the matrix A. N >= 0.
97 *
98 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
99 * On entry, the N-by-N matrix A.
100 * On exit, A has been overwritten. If JOBVL = 'V' or
101 * JOBVR = 'V', A contains the Schur form of the balanced
102 * version of the matrix A.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of the array A. LDA >= max(1,N).
106 *
107 * W (output) COMPLEX*16 array, dimension (N)
108 * W contains the computed eigenvalues.
109 *
110 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
111 * If JOBVL = 'V', the left eigenvectors u(j) are stored one
112 * after another in the columns of VL, in the same order
113 * as their eigenvalues.
114 * If JOBVL = 'N', VL is not referenced.
115 * u(j) = VL(:,j), the j-th column of VL.
116 *
117 * LDVL (input) INTEGER
118 * The leading dimension of the array VL. LDVL >= 1; if
119 * JOBVL = 'V', LDVL >= N.
120 *
121 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
122 * If JOBVR = 'V', the right eigenvectors v(j) are stored one
123 * after another in the columns of VR, in the same order
124 * as their eigenvalues.
125 * If JOBVR = 'N', VR is not referenced.
126 * v(j) = VR(:,j), the j-th column of VR.
127 *
128 * LDVR (input) INTEGER
129 * The leading dimension of the array VR. LDVR >= 1; if
130 * JOBVR = 'V', LDVR >= N.
131 *
132 * ILO (output) INTEGER
133 * IHI (output) INTEGER
134 * ILO and IHI are integer values determined when A was
135 * balanced. The balanced A(i,j) = 0 if I > J and
136 * J = 1,...,ILO-1 or I = IHI+1,...,N.
137 *
138 * SCALE (output) DOUBLE PRECISION array, dimension (N)
139 * Details of the permutations and scaling factors applied
140 * when balancing A. If P(j) is the index of the row and column
141 * interchanged with row and column j, and D(j) is the scaling
142 * factor applied to row and column j, then
143 * SCALE(J) = P(J), for J = 1,...,ILO-1
144 * = D(J), for J = ILO,...,IHI
145 * = P(J) for J = IHI+1,...,N.
146 * The order in which the interchanges are made is N to IHI+1,
147 * then 1 to ILO-1.
148 *
149 * ABNRM (output) DOUBLE PRECISION
150 * The one-norm of the balanced matrix (the maximum
151 * of the sum of absolute values of elements of any column).
152 *
153 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
154 * RCONDE(j) is the reciprocal condition number of the j-th
155 * eigenvalue.
156 *
157 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
158 * RCONDV(j) is the reciprocal condition number of the j-th
159 * right eigenvector.
160 *
161 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
162 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
163 *
164 * LWORK (input) INTEGER
165 * The dimension of the array WORK. If SENSE = 'N' or 'E',
166 * LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
167 * LWORK >= N*N+2*N.
168 * For good performance, LWORK must generally be larger.
169 *
170 * If LWORK = -1, then a workspace query is assumed; the routine
171 * only calculates the optimal size of the WORK array, returns
172 * this value as the first entry of the WORK array, and no error
173 * message related to LWORK is issued by XERBLA.
174 *
175 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
176 *
177 * INFO (output) INTEGER
178 * = 0: successful exit
179 * < 0: if INFO = -i, the i-th argument had an illegal value.
180 * > 0: if INFO = i, the QR algorithm failed to compute all the
181 * eigenvalues, and no eigenvectors or condition numbers
182 * have been computed; elements 1:ILO-1 and i+1:N of W
183 * contain eigenvalues which have converged.
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188 DOUBLE PRECISION ZERO, ONE
189 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
190 * ..
191 * .. Local Scalars ..
192 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
193 $ WNTSNN, WNTSNV
194 CHARACTER JOB, SIDE
195 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
196 $ MINWRK, NOUT
197 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
198 COMPLEX*16 TMP
199 * ..
200 * .. Local Arrays ..
201 LOGICAL SELECT( 1 )
202 DOUBLE PRECISION DUM( 1 )
203 * ..
204 * .. External Subroutines ..
205 EXTERNAL DLABAD, DLASCL, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL,
206 $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC,
207 $ ZTRSNA, ZUNGHR
208 * ..
209 * .. External Functions ..
210 LOGICAL LSAME
211 INTEGER IDAMAX, ILAENV
212 DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE
213 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
214 * ..
215 * .. Intrinsic Functions ..
216 INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
217 * ..
218 * .. Executable Statements ..
219 *
220 * Test the input arguments
221 *
222 INFO = 0
223 LQUERY = ( LWORK.EQ.-1 )
224 WANTVL = LSAME( JOBVL, 'V' )
225 WANTVR = LSAME( JOBVR, 'V' )
226 WNTSNN = LSAME( SENSE, 'N' )
227 WNTSNE = LSAME( SENSE, 'E' )
228 WNTSNV = LSAME( SENSE, 'V' )
229 WNTSNB = LSAME( SENSE, 'B' )
230 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR.
231 $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN
232 INFO = -1
233 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
234 INFO = -2
235 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
236 INFO = -3
237 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
238 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
239 $ WANTVR ) ) ) THEN
240 INFO = -4
241 ELSE IF( N.LT.0 ) THEN
242 INFO = -5
243 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
244 INFO = -7
245 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
246 INFO = -10
247 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
248 INFO = -12
249 END IF
250 *
251 * Compute workspace
252 * (Note: Comments in the code beginning "Workspace:" describe the
253 * minimal amount of workspace needed at that point in the code,
254 * as well as the preferred amount for good performance.
255 * CWorkspace refers to complex workspace, and RWorkspace to real
256 * workspace. NB refers to the optimal block size for the
257 * immediately following subroutine, as returned by ILAENV.
258 * HSWORK refers to the workspace preferred by ZHSEQR, as
259 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
260 * the worst case.)
261 *
262 IF( INFO.EQ.0 ) THEN
263 IF( N.EQ.0 ) THEN
264 MINWRK = 1
265 MAXWRK = 1
266 ELSE
267 MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
268 *
269 IF( WANTVL ) THEN
270 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
271 $ WORK, -1, INFO )
272 ELSE IF( WANTVR ) THEN
273 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
274 $ WORK, -1, INFO )
275 ELSE
276 IF( WNTSNN ) THEN
277 CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
278 $ WORK, -1, INFO )
279 ELSE
280 CALL ZHSEQR( 'S', 'N', N, 1, N, A, LDA, W, VR, LDVR,
281 $ WORK, -1, INFO )
282 END IF
283 END IF
284 HSWORK = WORK( 1 )
285 *
286 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
287 MINWRK = 2*N
288 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
289 $ MINWRK = MAX( MINWRK, N*N + 2*N )
290 MAXWRK = MAX( MAXWRK, HSWORK )
291 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
292 $ MAXWRK = MAX( MAXWRK, N*N + 2*N )
293 ELSE
294 MINWRK = 2*N
295 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
296 $ MINWRK = MAX( MINWRK, N*N + 2*N )
297 MAXWRK = MAX( MAXWRK, HSWORK )
298 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
299 $ ' ', N, 1, N, -1 ) )
300 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
301 $ MAXWRK = MAX( MAXWRK, N*N + 2*N )
302 MAXWRK = MAX( MAXWRK, 2*N )
303 END IF
304 MAXWRK = MAX( MAXWRK, MINWRK )
305 END IF
306 WORK( 1 ) = MAXWRK
307 *
308 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
309 INFO = -20
310 END IF
311 END IF
312 *
313 IF( INFO.NE.0 ) THEN
314 CALL XERBLA( 'ZGEEVX', -INFO )
315 RETURN
316 ELSE IF( LQUERY ) THEN
317 RETURN
318 END IF
319 *
320 * Quick return if possible
321 *
322 IF( N.EQ.0 )
323 $ RETURN
324 *
325 * Get machine constants
326 *
327 EPS = DLAMCH( 'P' )
328 SMLNUM = DLAMCH( 'S' )
329 BIGNUM = ONE / SMLNUM
330 CALL DLABAD( SMLNUM, BIGNUM )
331 SMLNUM = SQRT( SMLNUM ) / EPS
332 BIGNUM = ONE / SMLNUM
333 *
334 * Scale A if max element outside range [SMLNUM,BIGNUM]
335 *
336 ICOND = 0
337 ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
338 SCALEA = .FALSE.
339 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
340 SCALEA = .TRUE.
341 CSCALE = SMLNUM
342 ELSE IF( ANRM.GT.BIGNUM ) THEN
343 SCALEA = .TRUE.
344 CSCALE = BIGNUM
345 END IF
346 IF( SCALEA )
347 $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
348 *
349 * Balance the matrix and compute ABNRM
350 *
351 CALL ZGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
352 ABNRM = ZLANGE( '1', N, N, A, LDA, DUM )
353 IF( SCALEA ) THEN
354 DUM( 1 ) = ABNRM
355 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
356 ABNRM = DUM( 1 )
357 END IF
358 *
359 * Reduce to upper Hessenberg form
360 * (CWorkspace: need 2*N, prefer N+N*NB)
361 * (RWorkspace: none)
362 *
363 ITAU = 1
364 IWRK = ITAU + N
365 CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
366 $ LWORK-IWRK+1, IERR )
367 *
368 IF( WANTVL ) THEN
369 *
370 * Want left eigenvectors
371 * Copy Householder vectors to VL
372 *
373 SIDE = 'L'
374 CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
375 *
376 * Generate unitary matrix in VL
377 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
378 * (RWorkspace: none)
379 *
380 CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
381 $ LWORK-IWRK+1, IERR )
382 *
383 * Perform QR iteration, accumulating Schur vectors in VL
384 * (CWorkspace: need 1, prefer HSWORK (see comments) )
385 * (RWorkspace: none)
386 *
387 IWRK = ITAU
388 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
389 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
390 *
391 IF( WANTVR ) THEN
392 *
393 * Want left and right eigenvectors
394 * Copy Schur vectors to VR
395 *
396 SIDE = 'B'
397 CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
398 END IF
399 *
400 ELSE IF( WANTVR ) THEN
401 *
402 * Want right eigenvectors
403 * Copy Householder vectors to VR
404 *
405 SIDE = 'R'
406 CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
407 *
408 * Generate unitary matrix in VR
409 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
410 * (RWorkspace: none)
411 *
412 CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
413 $ LWORK-IWRK+1, IERR )
414 *
415 * Perform QR iteration, accumulating Schur vectors in VR
416 * (CWorkspace: need 1, prefer HSWORK (see comments) )
417 * (RWorkspace: none)
418 *
419 IWRK = ITAU
420 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
421 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
422 *
423 ELSE
424 *
425 * Compute eigenvalues only
426 * If condition numbers desired, compute Schur form
427 *
428 IF( WNTSNN ) THEN
429 JOB = 'E'
430 ELSE
431 JOB = 'S'
432 END IF
433 *
434 * (CWorkspace: need 1, prefer HSWORK (see comments) )
435 * (RWorkspace: none)
436 *
437 IWRK = ITAU
438 CALL ZHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
439 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
440 END IF
441 *
442 * If INFO > 0 from ZHSEQR, then quit
443 *
444 IF( INFO.GT.0 )
445 $ GO TO 50
446 *
447 IF( WANTVL .OR. WANTVR ) THEN
448 *
449 * Compute left and/or right eigenvectors
450 * (CWorkspace: need 2*N)
451 * (RWorkspace: need N)
452 *
453 CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
454 $ N, NOUT, WORK( IWRK ), RWORK, IERR )
455 END IF
456 *
457 * Compute condition numbers if desired
458 * (CWorkspace: need N*N+2*N unless SENSE = 'E')
459 * (RWorkspace: need 2*N unless SENSE = 'E')
460 *
461 IF( .NOT.WNTSNN ) THEN
462 CALL ZTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
463 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, RWORK,
464 $ ICOND )
465 END IF
466 *
467 IF( WANTVL ) THEN
468 *
469 * Undo balancing of left eigenvectors
470 *
471 CALL ZGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
472 $ IERR )
473 *
474 * Normalize left eigenvectors and make largest component real
475 *
476 DO 20 I = 1, N
477 SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
478 CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
479 DO 10 K = 1, N
480 RWORK( K ) = DBLE( VL( K, I ) )**2 +
481 $ DIMAG( VL( K, I ) )**2
482 10 CONTINUE
483 K = IDAMAX( N, RWORK, 1 )
484 TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
485 CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
486 VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
487 20 CONTINUE
488 END IF
489 *
490 IF( WANTVR ) THEN
491 *
492 * Undo balancing of right eigenvectors
493 *
494 CALL ZGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
495 $ IERR )
496 *
497 * Normalize right eigenvectors and make largest component real
498 *
499 DO 40 I = 1, N
500 SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
501 CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
502 DO 30 K = 1, N
503 RWORK( K ) = DBLE( VR( K, I ) )**2 +
504 $ DIMAG( VR( K, I ) )**2
505 30 CONTINUE
506 K = IDAMAX( N, RWORK, 1 )
507 TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
508 CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
509 VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
510 40 CONTINUE
511 END IF
512 *
513 * Undo scaling if necessary
514 *
515 50 CONTINUE
516 IF( SCALEA ) THEN
517 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
518 $ MAX( N-INFO, 1 ), IERR )
519 IF( INFO.EQ.0 ) THEN
520 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
521 $ CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
522 $ IERR )
523 ELSE
524 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
525 END IF
526 END IF
527 *
528 WORK( 1 ) = MAXWRK
529 RETURN
530 *
531 * End of ZGEEVX
532 *
533 END
2 $ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
3 $ RCONDV, WORK, LWORK, RWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER BALANC, JOBVL, JOBVR, SENSE
12 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
13 DOUBLE PRECISION ABNRM
14 * ..
15 * .. Array Arguments ..
16 DOUBLE PRECISION RCONDE( * ), RCONDV( * ), RWORK( * ),
17 $ SCALE( * )
18 COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ W( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
26 * eigenvalues and, optionally, the left and/or right eigenvectors.
27 *
28 * Optionally also, it computes a balancing transformation to improve
29 * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
30 * SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
31 * (RCONDE), and reciprocal condition numbers for the right
32 * eigenvectors (RCONDV).
33 *
34 * The right eigenvector v(j) of A satisfies
35 * A * v(j) = lambda(j) * v(j)
36 * where lambda(j) is its eigenvalue.
37 * The left eigenvector u(j) of A satisfies
38 * u(j)**H * A = lambda(j) * u(j)**H
39 * where u(j)**H denotes the conjugate transpose of u(j).
40 *
41 * The computed eigenvectors are normalized to have Euclidean norm
42 * equal to 1 and largest component real.
43 *
44 * Balancing a matrix means permuting the rows and columns to make it
45 * more nearly upper triangular, and applying a diagonal similarity
46 * transformation D * A * D**(-1), where D is a diagonal matrix, to
47 * make its rows and columns closer in norm and the condition numbers
48 * of its eigenvalues and eigenvectors smaller. The computed
49 * reciprocal condition numbers correspond to the balanced matrix.
50 * Permuting rows and columns will not change the condition numbers
51 * (in exact arithmetic) but diagonal scaling will. For further
52 * explanation of balancing, see section 4.10.2 of the LAPACK
53 * Users' Guide.
54 *
55 * Arguments
56 * =========
57 *
58 * BALANC (input) CHARACTER*1
59 * Indicates how the input matrix should be diagonally scaled
60 * and/or permuted to improve the conditioning of its
61 * eigenvalues.
62 * = 'N': Do not diagonally scale or permute;
63 * = 'P': Perform permutations to make the matrix more nearly
64 * upper triangular. Do not diagonally scale;
65 * = 'S': Diagonally scale the matrix, ie. replace A by
66 * D*A*D**(-1), where D is a diagonal matrix chosen
67 * to make the rows and columns of A more equal in
68 * norm. Do not permute;
69 * = 'B': Both diagonally scale and permute A.
70 *
71 * Computed reciprocal condition numbers will be for the matrix
72 * after balancing and/or permuting. Permuting does not change
73 * condition numbers (in exact arithmetic), but balancing does.
74 *
75 * JOBVL (input) CHARACTER*1
76 * = 'N': left eigenvectors of A are not computed;
77 * = 'V': left eigenvectors of A are computed.
78 * If SENSE = 'E' or 'B', JOBVL must = 'V'.
79 *
80 * JOBVR (input) CHARACTER*1
81 * = 'N': right eigenvectors of A are not computed;
82 * = 'V': right eigenvectors of A are computed.
83 * If SENSE = 'E' or 'B', JOBVR must = 'V'.
84 *
85 * SENSE (input) CHARACTER*1
86 * Determines which reciprocal condition numbers are computed.
87 * = 'N': None are computed;
88 * = 'E': Computed for eigenvalues only;
89 * = 'V': Computed for right eigenvectors only;
90 * = 'B': Computed for eigenvalues and right eigenvectors.
91 *
92 * If SENSE = 'E' or 'B', both left and right eigenvectors
93 * must also be computed (JOBVL = 'V' and JOBVR = 'V').
94 *
95 * N (input) INTEGER
96 * The order of the matrix A. N >= 0.
97 *
98 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
99 * On entry, the N-by-N matrix A.
100 * On exit, A has been overwritten. If JOBVL = 'V' or
101 * JOBVR = 'V', A contains the Schur form of the balanced
102 * version of the matrix A.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of the array A. LDA >= max(1,N).
106 *
107 * W (output) COMPLEX*16 array, dimension (N)
108 * W contains the computed eigenvalues.
109 *
110 * VL (output) COMPLEX*16 array, dimension (LDVL,N)
111 * If JOBVL = 'V', the left eigenvectors u(j) are stored one
112 * after another in the columns of VL, in the same order
113 * as their eigenvalues.
114 * If JOBVL = 'N', VL is not referenced.
115 * u(j) = VL(:,j), the j-th column of VL.
116 *
117 * LDVL (input) INTEGER
118 * The leading dimension of the array VL. LDVL >= 1; if
119 * JOBVL = 'V', LDVL >= N.
120 *
121 * VR (output) COMPLEX*16 array, dimension (LDVR,N)
122 * If JOBVR = 'V', the right eigenvectors v(j) are stored one
123 * after another in the columns of VR, in the same order
124 * as their eigenvalues.
125 * If JOBVR = 'N', VR is not referenced.
126 * v(j) = VR(:,j), the j-th column of VR.
127 *
128 * LDVR (input) INTEGER
129 * The leading dimension of the array VR. LDVR >= 1; if
130 * JOBVR = 'V', LDVR >= N.
131 *
132 * ILO (output) INTEGER
133 * IHI (output) INTEGER
134 * ILO and IHI are integer values determined when A was
135 * balanced. The balanced A(i,j) = 0 if I > J and
136 * J = 1,...,ILO-1 or I = IHI+1,...,N.
137 *
138 * SCALE (output) DOUBLE PRECISION array, dimension (N)
139 * Details of the permutations and scaling factors applied
140 * when balancing A. If P(j) is the index of the row and column
141 * interchanged with row and column j, and D(j) is the scaling
142 * factor applied to row and column j, then
143 * SCALE(J) = P(J), for J = 1,...,ILO-1
144 * = D(J), for J = ILO,...,IHI
145 * = P(J) for J = IHI+1,...,N.
146 * The order in which the interchanges are made is N to IHI+1,
147 * then 1 to ILO-1.
148 *
149 * ABNRM (output) DOUBLE PRECISION
150 * The one-norm of the balanced matrix (the maximum
151 * of the sum of absolute values of elements of any column).
152 *
153 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
154 * RCONDE(j) is the reciprocal condition number of the j-th
155 * eigenvalue.
156 *
157 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
158 * RCONDV(j) is the reciprocal condition number of the j-th
159 * right eigenvector.
160 *
161 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
162 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
163 *
164 * LWORK (input) INTEGER
165 * The dimension of the array WORK. If SENSE = 'N' or 'E',
166 * LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
167 * LWORK >= N*N+2*N.
168 * For good performance, LWORK must generally be larger.
169 *
170 * If LWORK = -1, then a workspace query is assumed; the routine
171 * only calculates the optimal size of the WORK array, returns
172 * this value as the first entry of the WORK array, and no error
173 * message related to LWORK is issued by XERBLA.
174 *
175 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
176 *
177 * INFO (output) INTEGER
178 * = 0: successful exit
179 * < 0: if INFO = -i, the i-th argument had an illegal value.
180 * > 0: if INFO = i, the QR algorithm failed to compute all the
181 * eigenvalues, and no eigenvectors or condition numbers
182 * have been computed; elements 1:ILO-1 and i+1:N of W
183 * contain eigenvalues which have converged.
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188 DOUBLE PRECISION ZERO, ONE
189 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
190 * ..
191 * .. Local Scalars ..
192 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
193 $ WNTSNN, WNTSNV
194 CHARACTER JOB, SIDE
195 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
196 $ MINWRK, NOUT
197 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
198 COMPLEX*16 TMP
199 * ..
200 * .. Local Arrays ..
201 LOGICAL SELECT( 1 )
202 DOUBLE PRECISION DUM( 1 )
203 * ..
204 * .. External Subroutines ..
205 EXTERNAL DLABAD, DLASCL, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL,
206 $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC,
207 $ ZTRSNA, ZUNGHR
208 * ..
209 * .. External Functions ..
210 LOGICAL LSAME
211 INTEGER IDAMAX, ILAENV
212 DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE
213 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
214 * ..
215 * .. Intrinsic Functions ..
216 INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
217 * ..
218 * .. Executable Statements ..
219 *
220 * Test the input arguments
221 *
222 INFO = 0
223 LQUERY = ( LWORK.EQ.-1 )
224 WANTVL = LSAME( JOBVL, 'V' )
225 WANTVR = LSAME( JOBVR, 'V' )
226 WNTSNN = LSAME( SENSE, 'N' )
227 WNTSNE = LSAME( SENSE, 'E' )
228 WNTSNV = LSAME( SENSE, 'V' )
229 WNTSNB = LSAME( SENSE, 'B' )
230 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR.
231 $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN
232 INFO = -1
233 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
234 INFO = -2
235 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
236 INFO = -3
237 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
238 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
239 $ WANTVR ) ) ) THEN
240 INFO = -4
241 ELSE IF( N.LT.0 ) THEN
242 INFO = -5
243 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
244 INFO = -7
245 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
246 INFO = -10
247 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
248 INFO = -12
249 END IF
250 *
251 * Compute workspace
252 * (Note: Comments in the code beginning "Workspace:" describe the
253 * minimal amount of workspace needed at that point in the code,
254 * as well as the preferred amount for good performance.
255 * CWorkspace refers to complex workspace, and RWorkspace to real
256 * workspace. NB refers to the optimal block size for the
257 * immediately following subroutine, as returned by ILAENV.
258 * HSWORK refers to the workspace preferred by ZHSEQR, as
259 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
260 * the worst case.)
261 *
262 IF( INFO.EQ.0 ) THEN
263 IF( N.EQ.0 ) THEN
264 MINWRK = 1
265 MAXWRK = 1
266 ELSE
267 MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
268 *
269 IF( WANTVL ) THEN
270 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
271 $ WORK, -1, INFO )
272 ELSE IF( WANTVR ) THEN
273 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
274 $ WORK, -1, INFO )
275 ELSE
276 IF( WNTSNN ) THEN
277 CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
278 $ WORK, -1, INFO )
279 ELSE
280 CALL ZHSEQR( 'S', 'N', N, 1, N, A, LDA, W, VR, LDVR,
281 $ WORK, -1, INFO )
282 END IF
283 END IF
284 HSWORK = WORK( 1 )
285 *
286 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
287 MINWRK = 2*N
288 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
289 $ MINWRK = MAX( MINWRK, N*N + 2*N )
290 MAXWRK = MAX( MAXWRK, HSWORK )
291 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
292 $ MAXWRK = MAX( MAXWRK, N*N + 2*N )
293 ELSE
294 MINWRK = 2*N
295 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
296 $ MINWRK = MAX( MINWRK, N*N + 2*N )
297 MAXWRK = MAX( MAXWRK, HSWORK )
298 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
299 $ ' ', N, 1, N, -1 ) )
300 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
301 $ MAXWRK = MAX( MAXWRK, N*N + 2*N )
302 MAXWRK = MAX( MAXWRK, 2*N )
303 END IF
304 MAXWRK = MAX( MAXWRK, MINWRK )
305 END IF
306 WORK( 1 ) = MAXWRK
307 *
308 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
309 INFO = -20
310 END IF
311 END IF
312 *
313 IF( INFO.NE.0 ) THEN
314 CALL XERBLA( 'ZGEEVX', -INFO )
315 RETURN
316 ELSE IF( LQUERY ) THEN
317 RETURN
318 END IF
319 *
320 * Quick return if possible
321 *
322 IF( N.EQ.0 )
323 $ RETURN
324 *
325 * Get machine constants
326 *
327 EPS = DLAMCH( 'P' )
328 SMLNUM = DLAMCH( 'S' )
329 BIGNUM = ONE / SMLNUM
330 CALL DLABAD( SMLNUM, BIGNUM )
331 SMLNUM = SQRT( SMLNUM ) / EPS
332 BIGNUM = ONE / SMLNUM
333 *
334 * Scale A if max element outside range [SMLNUM,BIGNUM]
335 *
336 ICOND = 0
337 ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
338 SCALEA = .FALSE.
339 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
340 SCALEA = .TRUE.
341 CSCALE = SMLNUM
342 ELSE IF( ANRM.GT.BIGNUM ) THEN
343 SCALEA = .TRUE.
344 CSCALE = BIGNUM
345 END IF
346 IF( SCALEA )
347 $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
348 *
349 * Balance the matrix and compute ABNRM
350 *
351 CALL ZGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
352 ABNRM = ZLANGE( '1', N, N, A, LDA, DUM )
353 IF( SCALEA ) THEN
354 DUM( 1 ) = ABNRM
355 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
356 ABNRM = DUM( 1 )
357 END IF
358 *
359 * Reduce to upper Hessenberg form
360 * (CWorkspace: need 2*N, prefer N+N*NB)
361 * (RWorkspace: none)
362 *
363 ITAU = 1
364 IWRK = ITAU + N
365 CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
366 $ LWORK-IWRK+1, IERR )
367 *
368 IF( WANTVL ) THEN
369 *
370 * Want left eigenvectors
371 * Copy Householder vectors to VL
372 *
373 SIDE = 'L'
374 CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
375 *
376 * Generate unitary matrix in VL
377 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
378 * (RWorkspace: none)
379 *
380 CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
381 $ LWORK-IWRK+1, IERR )
382 *
383 * Perform QR iteration, accumulating Schur vectors in VL
384 * (CWorkspace: need 1, prefer HSWORK (see comments) )
385 * (RWorkspace: none)
386 *
387 IWRK = ITAU
388 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
389 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
390 *
391 IF( WANTVR ) THEN
392 *
393 * Want left and right eigenvectors
394 * Copy Schur vectors to VR
395 *
396 SIDE = 'B'
397 CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
398 END IF
399 *
400 ELSE IF( WANTVR ) THEN
401 *
402 * Want right eigenvectors
403 * Copy Householder vectors to VR
404 *
405 SIDE = 'R'
406 CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
407 *
408 * Generate unitary matrix in VR
409 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
410 * (RWorkspace: none)
411 *
412 CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
413 $ LWORK-IWRK+1, IERR )
414 *
415 * Perform QR iteration, accumulating Schur vectors in VR
416 * (CWorkspace: need 1, prefer HSWORK (see comments) )
417 * (RWorkspace: none)
418 *
419 IWRK = ITAU
420 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
421 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
422 *
423 ELSE
424 *
425 * Compute eigenvalues only
426 * If condition numbers desired, compute Schur form
427 *
428 IF( WNTSNN ) THEN
429 JOB = 'E'
430 ELSE
431 JOB = 'S'
432 END IF
433 *
434 * (CWorkspace: need 1, prefer HSWORK (see comments) )
435 * (RWorkspace: none)
436 *
437 IWRK = ITAU
438 CALL ZHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
439 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
440 END IF
441 *
442 * If INFO > 0 from ZHSEQR, then quit
443 *
444 IF( INFO.GT.0 )
445 $ GO TO 50
446 *
447 IF( WANTVL .OR. WANTVR ) THEN
448 *
449 * Compute left and/or right eigenvectors
450 * (CWorkspace: need 2*N)
451 * (RWorkspace: need N)
452 *
453 CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
454 $ N, NOUT, WORK( IWRK ), RWORK, IERR )
455 END IF
456 *
457 * Compute condition numbers if desired
458 * (CWorkspace: need N*N+2*N unless SENSE = 'E')
459 * (RWorkspace: need 2*N unless SENSE = 'E')
460 *
461 IF( .NOT.WNTSNN ) THEN
462 CALL ZTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
463 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, RWORK,
464 $ ICOND )
465 END IF
466 *
467 IF( WANTVL ) THEN
468 *
469 * Undo balancing of left eigenvectors
470 *
471 CALL ZGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
472 $ IERR )
473 *
474 * Normalize left eigenvectors and make largest component real
475 *
476 DO 20 I = 1, N
477 SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
478 CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
479 DO 10 K = 1, N
480 RWORK( K ) = DBLE( VL( K, I ) )**2 +
481 $ DIMAG( VL( K, I ) )**2
482 10 CONTINUE
483 K = IDAMAX( N, RWORK, 1 )
484 TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
485 CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
486 VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
487 20 CONTINUE
488 END IF
489 *
490 IF( WANTVR ) THEN
491 *
492 * Undo balancing of right eigenvectors
493 *
494 CALL ZGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
495 $ IERR )
496 *
497 * Normalize right eigenvectors and make largest component real
498 *
499 DO 40 I = 1, N
500 SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
501 CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
502 DO 30 K = 1, N
503 RWORK( K ) = DBLE( VR( K, I ) )**2 +
504 $ DIMAG( VR( K, I ) )**2
505 30 CONTINUE
506 K = IDAMAX( N, RWORK, 1 )
507 TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
508 CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
509 VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
510 40 CONTINUE
511 END IF
512 *
513 * Undo scaling if necessary
514 *
515 50 CONTINUE
516 IF( SCALEA ) THEN
517 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
518 $ MAX( N-INFO, 1 ), IERR )
519 IF( INFO.EQ.0 ) THEN
520 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
521 $ CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
522 $ IERR )
523 ELSE
524 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
525 END IF
526 END IF
527 *
528 WORK( 1 ) = MAXWRK
529 RETURN
530 *
531 * End of ZGEEVX
532 *
533 END