1 SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
2 $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
3 $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
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 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ),
18 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ WI( * ), WORK( * ), WR( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGEEVX computes for an N-by-N real 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)**T * A = lambda(j) * u(j)**T
39 * where u(j)**T denotes the 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, i.e. 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) DOUBLE PRECISION 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 real Schur form of the balanced
102 * version of the input matrix A.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of the array A. LDA >= max(1,N).
106 *
107 * WR (output) DOUBLE PRECISION array, dimension (N)
108 * WI (output) DOUBLE PRECISION array, dimension (N)
109 * WR and WI contain the real and imaginary parts,
110 * respectively, of the computed eigenvalues. Complex
111 * conjugate pairs of eigenvalues will appear consecutively
112 * with the eigenvalue having the positive imaginary part
113 * first.
114 *
115 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
116 * If JOBVL = 'V', the left eigenvectors u(j) are stored one
117 * after another in the columns of VL, in the same order
118 * as their eigenvalues.
119 * If JOBVL = 'N', VL is not referenced.
120 * If the j-th eigenvalue is real, then u(j) = VL(:,j),
121 * the j-th column of VL.
122 * If the j-th and (j+1)-st eigenvalues form a complex
123 * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
124 * u(j+1) = VL(:,j) - i*VL(:,j+1).
125 *
126 * LDVL (input) INTEGER
127 * The leading dimension of the array VL. LDVL >= 1; if
128 * JOBVL = 'V', LDVL >= N.
129 *
130 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
131 * If JOBVR = 'V', the right eigenvectors v(j) are stored one
132 * after another in the columns of VR, in the same order
133 * as their eigenvalues.
134 * If JOBVR = 'N', VR is not referenced.
135 * If the j-th eigenvalue is real, then v(j) = VR(:,j),
136 * the j-th column of VR.
137 * If the j-th and (j+1)-st eigenvalues form a complex
138 * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
139 * v(j+1) = VR(:,j) - i*VR(:,j+1).
140 *
141 * LDVR (input) INTEGER
142 * The leading dimension of the array VR. LDVR >= 1, and if
143 * JOBVR = 'V', LDVR >= N.
144 *
145 * ILO (output) INTEGER
146 * IHI (output) INTEGER
147 * ILO and IHI are integer values determined when A was
148 * balanced. The balanced A(i,j) = 0 if I > J and
149 * J = 1,...,ILO-1 or I = IHI+1,...,N.
150 *
151 * SCALE (output) DOUBLE PRECISION array, dimension (N)
152 * Details of the permutations and scaling factors applied
153 * when balancing A. If P(j) is the index of the row and column
154 * interchanged with row and column j, and D(j) is the scaling
155 * factor applied to row and column j, then
156 * SCALE(J) = P(J), for J = 1,...,ILO-1
157 * = D(J), for J = ILO,...,IHI
158 * = P(J) for J = IHI+1,...,N.
159 * The order in which the interchanges are made is N to IHI+1,
160 * then 1 to ILO-1.
161 *
162 * ABNRM (output) DOUBLE PRECISION
163 * The one-norm of the balanced matrix (the maximum
164 * of the sum of absolute values of elements of any column).
165 *
166 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
167 * RCONDE(j) is the reciprocal condition number of the j-th
168 * eigenvalue.
169 *
170 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
171 * RCONDV(j) is the reciprocal condition number of the j-th
172 * right eigenvector.
173 *
174 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
175 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
176 *
177 * LWORK (input) INTEGER
178 * The dimension of the array WORK. If SENSE = 'N' or 'E',
179 * LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
180 * LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6).
181 * For good performance, LWORK must generally be larger.
182 *
183 * If LWORK = -1, then a workspace query is assumed; the routine
184 * only calculates the optimal size of the WORK array, returns
185 * this value as the first entry of the WORK array, and no error
186 * message related to LWORK is issued by XERBLA.
187 *
188 * IWORK (workspace) INTEGER array, dimension (2*N-2)
189 * If SENSE = 'N' or 'E', not referenced.
190 *
191 * INFO (output) INTEGER
192 * = 0: successful exit
193 * < 0: if INFO = -i, the i-th argument had an illegal value.
194 * > 0: if INFO = i, the QR algorithm failed to compute all the
195 * eigenvalues, and no eigenvectors or condition numbers
196 * have been computed; elements 1:ILO-1 and i+1:N of WR
197 * and WI contain eigenvalues which have converged.
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202 DOUBLE PRECISION ZERO, ONE
203 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
204 * ..
205 * .. Local Scalars ..
206 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
207 $ WNTSNN, WNTSNV
208 CHARACTER JOB, SIDE
209 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
210 $ MINWRK, NOUT
211 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
212 $ SN
213 * ..
214 * .. Local Arrays ..
215 LOGICAL SELECT( 1 )
216 DOUBLE PRECISION DUM( 1 )
217 * ..
218 * .. External Subroutines ..
219 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
220 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
221 $ DTRSNA, XERBLA
222 * ..
223 * .. External Functions ..
224 LOGICAL LSAME
225 INTEGER IDAMAX, ILAENV
226 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
227 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
228 $ DNRM2
229 * ..
230 * .. Intrinsic Functions ..
231 INTRINSIC MAX, SQRT
232 * ..
233 * .. Executable Statements ..
234 *
235 * Test the input arguments
236 *
237 INFO = 0
238 LQUERY = ( LWORK.EQ.-1 )
239 WANTVL = LSAME( JOBVL, 'V' )
240 WANTVR = LSAME( JOBVR, 'V' )
241 WNTSNN = LSAME( SENSE, 'N' )
242 WNTSNE = LSAME( SENSE, 'E' )
243 WNTSNV = LSAME( SENSE, 'V' )
244 WNTSNB = LSAME( SENSE, 'B' )
245 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
246 $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
247 $ THEN
248 INFO = -1
249 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
250 INFO = -2
251 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
252 INFO = -3
253 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
254 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
255 $ WANTVR ) ) ) THEN
256 INFO = -4
257 ELSE IF( N.LT.0 ) THEN
258 INFO = -5
259 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
260 INFO = -7
261 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
262 INFO = -11
263 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
264 INFO = -13
265 END IF
266 *
267 * Compute workspace
268 * (Note: Comments in the code beginning "Workspace:" describe the
269 * minimal amount of workspace needed at that point in the code,
270 * as well as the preferred amount for good performance.
271 * NB refers to the optimal block size for the immediately
272 * following subroutine, as returned by ILAENV.
273 * HSWORK refers to the workspace preferred by DHSEQR, as
274 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
275 * the worst case.)
276 *
277 IF( INFO.EQ.0 ) THEN
278 IF( N.EQ.0 ) THEN
279 MINWRK = 1
280 MAXWRK = 1
281 ELSE
282 MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
283 *
284 IF( WANTVL ) THEN
285 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
286 $ WORK, -1, INFO )
287 ELSE IF( WANTVR ) THEN
288 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
289 $ WORK, -1, INFO )
290 ELSE
291 IF( WNTSNN ) THEN
292 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR,
293 $ LDVR, WORK, -1, INFO )
294 ELSE
295 CALL DHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR,
296 $ LDVR, WORK, -1, INFO )
297 END IF
298 END IF
299 HSWORK = WORK( 1 )
300 *
301 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
302 MINWRK = 2*N
303 IF( .NOT.WNTSNN )
304 $ MINWRK = MAX( MINWRK, N*N+6*N )
305 MAXWRK = MAX( MAXWRK, HSWORK )
306 IF( .NOT.WNTSNN )
307 $ MAXWRK = MAX( MAXWRK, N*N + 6*N )
308 ELSE
309 MINWRK = 3*N
310 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
311 $ MINWRK = MAX( MINWRK, N*N + 6*N )
312 MAXWRK = MAX( MAXWRK, HSWORK )
313 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'DORGHR',
314 $ ' ', N, 1, N, -1 ) )
315 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
316 $ MAXWRK = MAX( MAXWRK, N*N + 6*N )
317 MAXWRK = MAX( MAXWRK, 3*N )
318 END IF
319 MAXWRK = MAX( MAXWRK, MINWRK )
320 END IF
321 WORK( 1 ) = MAXWRK
322 *
323 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
324 INFO = -21
325 END IF
326 END IF
327 *
328 IF( INFO.NE.0 ) THEN
329 CALL XERBLA( 'DGEEVX', -INFO )
330 RETURN
331 ELSE IF( LQUERY ) THEN
332 RETURN
333 END IF
334 *
335 * Quick return if possible
336 *
337 IF( N.EQ.0 )
338 $ RETURN
339 *
340 * Get machine constants
341 *
342 EPS = DLAMCH( 'P' )
343 SMLNUM = DLAMCH( 'S' )
344 BIGNUM = ONE / SMLNUM
345 CALL DLABAD( SMLNUM, BIGNUM )
346 SMLNUM = SQRT( SMLNUM ) / EPS
347 BIGNUM = ONE / SMLNUM
348 *
349 * Scale A if max element outside range [SMLNUM,BIGNUM]
350 *
351 ICOND = 0
352 ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
353 SCALEA = .FALSE.
354 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
355 SCALEA = .TRUE.
356 CSCALE = SMLNUM
357 ELSE IF( ANRM.GT.BIGNUM ) THEN
358 SCALEA = .TRUE.
359 CSCALE = BIGNUM
360 END IF
361 IF( SCALEA )
362 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
363 *
364 * Balance the matrix and compute ABNRM
365 *
366 CALL DGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
367 ABNRM = DLANGE( '1', N, N, A, LDA, DUM )
368 IF( SCALEA ) THEN
369 DUM( 1 ) = ABNRM
370 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
371 ABNRM = DUM( 1 )
372 END IF
373 *
374 * Reduce to upper Hessenberg form
375 * (Workspace: need 2*N, prefer N+N*NB)
376 *
377 ITAU = 1
378 IWRK = ITAU + N
379 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
380 $ LWORK-IWRK+1, IERR )
381 *
382 IF( WANTVL ) THEN
383 *
384 * Want left eigenvectors
385 * Copy Householder vectors to VL
386 *
387 SIDE = 'L'
388 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
389 *
390 * Generate orthogonal matrix in VL
391 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
392 *
393 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
394 $ LWORK-IWRK+1, IERR )
395 *
396 * Perform QR iteration, accumulating Schur vectors in VL
397 * (Workspace: need 1, prefer HSWORK (see comments) )
398 *
399 IWRK = ITAU
400 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
401 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
402 *
403 IF( WANTVR ) THEN
404 *
405 * Want left and right eigenvectors
406 * Copy Schur vectors to VR
407 *
408 SIDE = 'B'
409 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
410 END IF
411 *
412 ELSE IF( WANTVR ) THEN
413 *
414 * Want right eigenvectors
415 * Copy Householder vectors to VR
416 *
417 SIDE = 'R'
418 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
419 *
420 * Generate orthogonal matrix in VR
421 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
422 *
423 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
424 $ LWORK-IWRK+1, IERR )
425 *
426 * Perform QR iteration, accumulating Schur vectors in VR
427 * (Workspace: need 1, prefer HSWORK (see comments) )
428 *
429 IWRK = ITAU
430 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
431 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
432 *
433 ELSE
434 *
435 * Compute eigenvalues only
436 * If condition numbers desired, compute Schur form
437 *
438 IF( WNTSNN ) THEN
439 JOB = 'E'
440 ELSE
441 JOB = 'S'
442 END IF
443 *
444 * (Workspace: need 1, prefer HSWORK (see comments) )
445 *
446 IWRK = ITAU
447 CALL DHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
448 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
449 END IF
450 *
451 * If INFO > 0 from DHSEQR, then quit
452 *
453 IF( INFO.GT.0 )
454 $ GO TO 50
455 *
456 IF( WANTVL .OR. WANTVR ) THEN
457 *
458 * Compute left and/or right eigenvectors
459 * (Workspace: need 3*N)
460 *
461 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
462 $ N, NOUT, WORK( IWRK ), IERR )
463 END IF
464 *
465 * Compute condition numbers if desired
466 * (Workspace: need N*N+6*N unless SENSE = 'E')
467 *
468 IF( .NOT.WNTSNN ) THEN
469 CALL DTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
470 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK,
471 $ ICOND )
472 END IF
473 *
474 IF( WANTVL ) THEN
475 *
476 * Undo balancing of left eigenvectors
477 *
478 CALL DGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
479 $ IERR )
480 *
481 * Normalize left eigenvectors and make largest component real
482 *
483 DO 20 I = 1, N
484 IF( WI( I ).EQ.ZERO ) THEN
485 SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
486 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
487 ELSE IF( WI( I ).GT.ZERO ) THEN
488 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
489 $ DNRM2( N, VL( 1, I+1 ), 1 ) )
490 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
491 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
492 DO 10 K = 1, N
493 WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2
494 10 CONTINUE
495 K = IDAMAX( N, WORK, 1 )
496 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
497 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
498 VL( K, I+1 ) = ZERO
499 END IF
500 20 CONTINUE
501 END IF
502 *
503 IF( WANTVR ) THEN
504 *
505 * Undo balancing of right eigenvectors
506 *
507 CALL DGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
508 $ IERR )
509 *
510 * Normalize right eigenvectors and make largest component real
511 *
512 DO 40 I = 1, N
513 IF( WI( I ).EQ.ZERO ) THEN
514 SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
515 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
516 ELSE IF( WI( I ).GT.ZERO ) THEN
517 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
518 $ DNRM2( N, VR( 1, I+1 ), 1 ) )
519 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
520 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
521 DO 30 K = 1, N
522 WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2
523 30 CONTINUE
524 K = IDAMAX( N, WORK, 1 )
525 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
526 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
527 VR( K, I+1 ) = ZERO
528 END IF
529 40 CONTINUE
530 END IF
531 *
532 * Undo scaling if necessary
533 *
534 50 CONTINUE
535 IF( SCALEA ) THEN
536 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
537 $ MAX( N-INFO, 1 ), IERR )
538 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
539 $ MAX( N-INFO, 1 ), IERR )
540 IF( INFO.EQ.0 ) THEN
541 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
542 $ CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
543 $ IERR )
544 ELSE
545 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
546 $ IERR )
547 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
548 $ IERR )
549 END IF
550 END IF
551 *
552 WORK( 1 ) = MAXWRK
553 RETURN
554 *
555 * End of DGEEVX
556 *
557 END
2 $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
3 $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
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 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ),
18 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
19 $ WI( * ), WORK( * ), WR( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGEEVX computes for an N-by-N real 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)**T * A = lambda(j) * u(j)**T
39 * where u(j)**T denotes the 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, i.e. 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) DOUBLE PRECISION 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 real Schur form of the balanced
102 * version of the input matrix A.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of the array A. LDA >= max(1,N).
106 *
107 * WR (output) DOUBLE PRECISION array, dimension (N)
108 * WI (output) DOUBLE PRECISION array, dimension (N)
109 * WR and WI contain the real and imaginary parts,
110 * respectively, of the computed eigenvalues. Complex
111 * conjugate pairs of eigenvalues will appear consecutively
112 * with the eigenvalue having the positive imaginary part
113 * first.
114 *
115 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
116 * If JOBVL = 'V', the left eigenvectors u(j) are stored one
117 * after another in the columns of VL, in the same order
118 * as their eigenvalues.
119 * If JOBVL = 'N', VL is not referenced.
120 * If the j-th eigenvalue is real, then u(j) = VL(:,j),
121 * the j-th column of VL.
122 * If the j-th and (j+1)-st eigenvalues form a complex
123 * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
124 * u(j+1) = VL(:,j) - i*VL(:,j+1).
125 *
126 * LDVL (input) INTEGER
127 * The leading dimension of the array VL. LDVL >= 1; if
128 * JOBVL = 'V', LDVL >= N.
129 *
130 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
131 * If JOBVR = 'V', the right eigenvectors v(j) are stored one
132 * after another in the columns of VR, in the same order
133 * as their eigenvalues.
134 * If JOBVR = 'N', VR is not referenced.
135 * If the j-th eigenvalue is real, then v(j) = VR(:,j),
136 * the j-th column of VR.
137 * If the j-th and (j+1)-st eigenvalues form a complex
138 * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
139 * v(j+1) = VR(:,j) - i*VR(:,j+1).
140 *
141 * LDVR (input) INTEGER
142 * The leading dimension of the array VR. LDVR >= 1, and if
143 * JOBVR = 'V', LDVR >= N.
144 *
145 * ILO (output) INTEGER
146 * IHI (output) INTEGER
147 * ILO and IHI are integer values determined when A was
148 * balanced. The balanced A(i,j) = 0 if I > J and
149 * J = 1,...,ILO-1 or I = IHI+1,...,N.
150 *
151 * SCALE (output) DOUBLE PRECISION array, dimension (N)
152 * Details of the permutations and scaling factors applied
153 * when balancing A. If P(j) is the index of the row and column
154 * interchanged with row and column j, and D(j) is the scaling
155 * factor applied to row and column j, then
156 * SCALE(J) = P(J), for J = 1,...,ILO-1
157 * = D(J), for J = ILO,...,IHI
158 * = P(J) for J = IHI+1,...,N.
159 * The order in which the interchanges are made is N to IHI+1,
160 * then 1 to ILO-1.
161 *
162 * ABNRM (output) DOUBLE PRECISION
163 * The one-norm of the balanced matrix (the maximum
164 * of the sum of absolute values of elements of any column).
165 *
166 * RCONDE (output) DOUBLE PRECISION array, dimension (N)
167 * RCONDE(j) is the reciprocal condition number of the j-th
168 * eigenvalue.
169 *
170 * RCONDV (output) DOUBLE PRECISION array, dimension (N)
171 * RCONDV(j) is the reciprocal condition number of the j-th
172 * right eigenvector.
173 *
174 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
175 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
176 *
177 * LWORK (input) INTEGER
178 * The dimension of the array WORK. If SENSE = 'N' or 'E',
179 * LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
180 * LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6).
181 * For good performance, LWORK must generally be larger.
182 *
183 * If LWORK = -1, then a workspace query is assumed; the routine
184 * only calculates the optimal size of the WORK array, returns
185 * this value as the first entry of the WORK array, and no error
186 * message related to LWORK is issued by XERBLA.
187 *
188 * IWORK (workspace) INTEGER array, dimension (2*N-2)
189 * If SENSE = 'N' or 'E', not referenced.
190 *
191 * INFO (output) INTEGER
192 * = 0: successful exit
193 * < 0: if INFO = -i, the i-th argument had an illegal value.
194 * > 0: if INFO = i, the QR algorithm failed to compute all the
195 * eigenvalues, and no eigenvectors or condition numbers
196 * have been computed; elements 1:ILO-1 and i+1:N of WR
197 * and WI contain eigenvalues which have converged.
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202 DOUBLE PRECISION ZERO, ONE
203 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
204 * ..
205 * .. Local Scalars ..
206 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
207 $ WNTSNN, WNTSNV
208 CHARACTER JOB, SIDE
209 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
210 $ MINWRK, NOUT
211 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
212 $ SN
213 * ..
214 * .. Local Arrays ..
215 LOGICAL SELECT( 1 )
216 DOUBLE PRECISION DUM( 1 )
217 * ..
218 * .. External Subroutines ..
219 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
220 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
221 $ DTRSNA, XERBLA
222 * ..
223 * .. External Functions ..
224 LOGICAL LSAME
225 INTEGER IDAMAX, ILAENV
226 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
227 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
228 $ DNRM2
229 * ..
230 * .. Intrinsic Functions ..
231 INTRINSIC MAX, SQRT
232 * ..
233 * .. Executable Statements ..
234 *
235 * Test the input arguments
236 *
237 INFO = 0
238 LQUERY = ( LWORK.EQ.-1 )
239 WANTVL = LSAME( JOBVL, 'V' )
240 WANTVR = LSAME( JOBVR, 'V' )
241 WNTSNN = LSAME( SENSE, 'N' )
242 WNTSNE = LSAME( SENSE, 'E' )
243 WNTSNV = LSAME( SENSE, 'V' )
244 WNTSNB = LSAME( SENSE, 'B' )
245 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
246 $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
247 $ THEN
248 INFO = -1
249 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
250 INFO = -2
251 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
252 INFO = -3
253 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
254 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
255 $ WANTVR ) ) ) THEN
256 INFO = -4
257 ELSE IF( N.LT.0 ) THEN
258 INFO = -5
259 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
260 INFO = -7
261 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
262 INFO = -11
263 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
264 INFO = -13
265 END IF
266 *
267 * Compute workspace
268 * (Note: Comments in the code beginning "Workspace:" describe the
269 * minimal amount of workspace needed at that point in the code,
270 * as well as the preferred amount for good performance.
271 * NB refers to the optimal block size for the immediately
272 * following subroutine, as returned by ILAENV.
273 * HSWORK refers to the workspace preferred by DHSEQR, as
274 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
275 * the worst case.)
276 *
277 IF( INFO.EQ.0 ) THEN
278 IF( N.EQ.0 ) THEN
279 MINWRK = 1
280 MAXWRK = 1
281 ELSE
282 MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
283 *
284 IF( WANTVL ) THEN
285 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
286 $ WORK, -1, INFO )
287 ELSE IF( WANTVR ) THEN
288 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
289 $ WORK, -1, INFO )
290 ELSE
291 IF( WNTSNN ) THEN
292 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR,
293 $ LDVR, WORK, -1, INFO )
294 ELSE
295 CALL DHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR,
296 $ LDVR, WORK, -1, INFO )
297 END IF
298 END IF
299 HSWORK = WORK( 1 )
300 *
301 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
302 MINWRK = 2*N
303 IF( .NOT.WNTSNN )
304 $ MINWRK = MAX( MINWRK, N*N+6*N )
305 MAXWRK = MAX( MAXWRK, HSWORK )
306 IF( .NOT.WNTSNN )
307 $ MAXWRK = MAX( MAXWRK, N*N + 6*N )
308 ELSE
309 MINWRK = 3*N
310 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
311 $ MINWRK = MAX( MINWRK, N*N + 6*N )
312 MAXWRK = MAX( MAXWRK, HSWORK )
313 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'DORGHR',
314 $ ' ', N, 1, N, -1 ) )
315 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
316 $ MAXWRK = MAX( MAXWRK, N*N + 6*N )
317 MAXWRK = MAX( MAXWRK, 3*N )
318 END IF
319 MAXWRK = MAX( MAXWRK, MINWRK )
320 END IF
321 WORK( 1 ) = MAXWRK
322 *
323 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
324 INFO = -21
325 END IF
326 END IF
327 *
328 IF( INFO.NE.0 ) THEN
329 CALL XERBLA( 'DGEEVX', -INFO )
330 RETURN
331 ELSE IF( LQUERY ) THEN
332 RETURN
333 END IF
334 *
335 * Quick return if possible
336 *
337 IF( N.EQ.0 )
338 $ RETURN
339 *
340 * Get machine constants
341 *
342 EPS = DLAMCH( 'P' )
343 SMLNUM = DLAMCH( 'S' )
344 BIGNUM = ONE / SMLNUM
345 CALL DLABAD( SMLNUM, BIGNUM )
346 SMLNUM = SQRT( SMLNUM ) / EPS
347 BIGNUM = ONE / SMLNUM
348 *
349 * Scale A if max element outside range [SMLNUM,BIGNUM]
350 *
351 ICOND = 0
352 ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
353 SCALEA = .FALSE.
354 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
355 SCALEA = .TRUE.
356 CSCALE = SMLNUM
357 ELSE IF( ANRM.GT.BIGNUM ) THEN
358 SCALEA = .TRUE.
359 CSCALE = BIGNUM
360 END IF
361 IF( SCALEA )
362 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
363 *
364 * Balance the matrix and compute ABNRM
365 *
366 CALL DGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
367 ABNRM = DLANGE( '1', N, N, A, LDA, DUM )
368 IF( SCALEA ) THEN
369 DUM( 1 ) = ABNRM
370 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
371 ABNRM = DUM( 1 )
372 END IF
373 *
374 * Reduce to upper Hessenberg form
375 * (Workspace: need 2*N, prefer N+N*NB)
376 *
377 ITAU = 1
378 IWRK = ITAU + N
379 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
380 $ LWORK-IWRK+1, IERR )
381 *
382 IF( WANTVL ) THEN
383 *
384 * Want left eigenvectors
385 * Copy Householder vectors to VL
386 *
387 SIDE = 'L'
388 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
389 *
390 * Generate orthogonal matrix in VL
391 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
392 *
393 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
394 $ LWORK-IWRK+1, IERR )
395 *
396 * Perform QR iteration, accumulating Schur vectors in VL
397 * (Workspace: need 1, prefer HSWORK (see comments) )
398 *
399 IWRK = ITAU
400 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
401 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
402 *
403 IF( WANTVR ) THEN
404 *
405 * Want left and right eigenvectors
406 * Copy Schur vectors to VR
407 *
408 SIDE = 'B'
409 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
410 END IF
411 *
412 ELSE IF( WANTVR ) THEN
413 *
414 * Want right eigenvectors
415 * Copy Householder vectors to VR
416 *
417 SIDE = 'R'
418 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
419 *
420 * Generate orthogonal matrix in VR
421 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
422 *
423 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
424 $ LWORK-IWRK+1, IERR )
425 *
426 * Perform QR iteration, accumulating Schur vectors in VR
427 * (Workspace: need 1, prefer HSWORK (see comments) )
428 *
429 IWRK = ITAU
430 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
431 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
432 *
433 ELSE
434 *
435 * Compute eigenvalues only
436 * If condition numbers desired, compute Schur form
437 *
438 IF( WNTSNN ) THEN
439 JOB = 'E'
440 ELSE
441 JOB = 'S'
442 END IF
443 *
444 * (Workspace: need 1, prefer HSWORK (see comments) )
445 *
446 IWRK = ITAU
447 CALL DHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
448 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
449 END IF
450 *
451 * If INFO > 0 from DHSEQR, then quit
452 *
453 IF( INFO.GT.0 )
454 $ GO TO 50
455 *
456 IF( WANTVL .OR. WANTVR ) THEN
457 *
458 * Compute left and/or right eigenvectors
459 * (Workspace: need 3*N)
460 *
461 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
462 $ N, NOUT, WORK( IWRK ), IERR )
463 END IF
464 *
465 * Compute condition numbers if desired
466 * (Workspace: need N*N+6*N unless SENSE = 'E')
467 *
468 IF( .NOT.WNTSNN ) THEN
469 CALL DTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
470 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK,
471 $ ICOND )
472 END IF
473 *
474 IF( WANTVL ) THEN
475 *
476 * Undo balancing of left eigenvectors
477 *
478 CALL DGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
479 $ IERR )
480 *
481 * Normalize left eigenvectors and make largest component real
482 *
483 DO 20 I = 1, N
484 IF( WI( I ).EQ.ZERO ) THEN
485 SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
486 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
487 ELSE IF( WI( I ).GT.ZERO ) THEN
488 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
489 $ DNRM2( N, VL( 1, I+1 ), 1 ) )
490 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
491 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
492 DO 10 K = 1, N
493 WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2
494 10 CONTINUE
495 K = IDAMAX( N, WORK, 1 )
496 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
497 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
498 VL( K, I+1 ) = ZERO
499 END IF
500 20 CONTINUE
501 END IF
502 *
503 IF( WANTVR ) THEN
504 *
505 * Undo balancing of right eigenvectors
506 *
507 CALL DGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
508 $ IERR )
509 *
510 * Normalize right eigenvectors and make largest component real
511 *
512 DO 40 I = 1, N
513 IF( WI( I ).EQ.ZERO ) THEN
514 SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
515 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
516 ELSE IF( WI( I ).GT.ZERO ) THEN
517 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
518 $ DNRM2( N, VR( 1, I+1 ), 1 ) )
519 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
520 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
521 DO 30 K = 1, N
522 WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2
523 30 CONTINUE
524 K = IDAMAX( N, WORK, 1 )
525 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
526 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
527 VR( K, I+1 ) = ZERO
528 END IF
529 40 CONTINUE
530 END IF
531 *
532 * Undo scaling if necessary
533 *
534 50 CONTINUE
535 IF( SCALEA ) THEN
536 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
537 $ MAX( N-INFO, 1 ), IERR )
538 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
539 $ MAX( N-INFO, 1 ), IERR )
540 IF( INFO.EQ.0 ) THEN
541 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
542 $ CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
543 $ IERR )
544 ELSE
545 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
546 $ IERR )
547 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
548 $ IERR )
549 END IF
550 END IF
551 *
552 WORK( 1 ) = MAXWRK
553 RETURN
554 *
555 * End of DGEEVX
556 *
557 END