LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgeevx.f
Go to the documentation of this file.
1 *> \brief <b> DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGEEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
22 * VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
23 * RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER BALANC, JOBVL, JOBVR, SENSE
27 * INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
28 * DOUBLE PRECISION ABNRM
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IWORK( * )
32 * DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ),
33 * $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
34 * $ WI( * ), WORK( * ), WR( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> DGEEVX computes for an N-by-N real nonsymmetric matrix A, the
44 *> eigenvalues and, optionally, the left and/or right eigenvectors.
45 *>
46 *> Optionally also, it computes a balancing transformation to improve
47 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
48 *> SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
49 *> (RCONDE), and reciprocal condition numbers for the right
50 *> eigenvectors (RCONDV).
51 *>
52 *> The right eigenvector v(j) of A satisfies
53 *> A * v(j) = lambda(j) * v(j)
54 *> where lambda(j) is its eigenvalue.
55 *> The left eigenvector u(j) of A satisfies
56 *> u(j)**H * A = lambda(j) * u(j)**H
57 *> where u(j)**H denotes the conjugate-transpose of u(j).
58 *>
59 *> The computed eigenvectors are normalized to have Euclidean norm
60 *> equal to 1 and largest component real.
61 *>
62 *> Balancing a matrix means permuting the rows and columns to make it
63 *> more nearly upper triangular, and applying a diagonal similarity
64 *> transformation D * A * D**(-1), where D is a diagonal matrix, to
65 *> make its rows and columns closer in norm and the condition numbers
66 *> of its eigenvalues and eigenvectors smaller. The computed
67 *> reciprocal condition numbers correspond to the balanced matrix.
68 *> Permuting rows and columns will not change the condition numbers
69 *> (in exact arithmetic) but diagonal scaling will. For further
70 *> explanation of balancing, see section 4.10.2 of the LAPACK
71 *> Users' Guide.
72 *> \endverbatim
73 *
74 * Arguments:
75 * ==========
76 *
77 *> \param[in] BALANC
78 *> \verbatim
79 *> BALANC is CHARACTER*1
80 *> Indicates how the input matrix should be diagonally scaled
81 *> and/or permuted to improve the conditioning of its
82 *> eigenvalues.
83 *> = 'N': Do not diagonally scale or permute;
84 *> = 'P': Perform permutations to make the matrix more nearly
85 *> upper triangular. Do not diagonally scale;
86 *> = 'S': Diagonally scale the matrix, i.e. replace A by
87 *> D*A*D**(-1), where D is a diagonal matrix chosen
88 *> to make the rows and columns of A more equal in
89 *> norm. Do not permute;
90 *> = 'B': Both diagonally scale and permute A.
91 *>
92 *> Computed reciprocal condition numbers will be for the matrix
93 *> after balancing and/or permuting. Permuting does not change
94 *> condition numbers (in exact arithmetic), but balancing does.
95 *> \endverbatim
96 *>
97 *> \param[in] JOBVL
98 *> \verbatim
99 *> JOBVL is CHARACTER*1
100 *> = 'N': left eigenvectors of A are not computed;
101 *> = 'V': left eigenvectors of A are computed.
102 *> If SENSE = 'E' or 'B', JOBVL must = 'V'.
103 *> \endverbatim
104 *>
105 *> \param[in] JOBVR
106 *> \verbatim
107 *> JOBVR is CHARACTER*1
108 *> = 'N': right eigenvectors of A are not computed;
109 *> = 'V': right eigenvectors of A are computed.
110 *> If SENSE = 'E' or 'B', JOBVR must = 'V'.
111 *> \endverbatim
112 *>
113 *> \param[in] SENSE
114 *> \verbatim
115 *> SENSE is CHARACTER*1
116 *> Determines which reciprocal condition numbers are computed.
117 *> = 'N': None are computed;
118 *> = 'E': Computed for eigenvalues only;
119 *> = 'V': Computed for right eigenvectors only;
120 *> = 'B': Computed for eigenvalues and right eigenvectors.
121 *>
122 *> If SENSE = 'E' or 'B', both left and right eigenvectors
123 *> must also be computed (JOBVL = 'V' and JOBVR = 'V').
124 *> \endverbatim
125 *>
126 *> \param[in] N
127 *> \verbatim
128 *> N is INTEGER
129 *> The order of the matrix A. N >= 0.
130 *> \endverbatim
131 *>
132 *> \param[in,out] A
133 *> \verbatim
134 *> A is DOUBLE PRECISION array, dimension (LDA,N)
135 *> On entry, the N-by-N matrix A.
136 *> On exit, A has been overwritten. If JOBVL = 'V' or
137 *> JOBVR = 'V', A contains the real Schur form of the balanced
138 *> version of the input matrix A.
139 *> \endverbatim
140 *>
141 *> \param[in] LDA
142 *> \verbatim
143 *> LDA is INTEGER
144 *> The leading dimension of the array A. LDA >= max(1,N).
145 *> \endverbatim
146 *>
147 *> \param[out] WR
148 *> \verbatim
149 *> WR is DOUBLE PRECISION array, dimension (N)
150 *> \endverbatim
151 *>
152 *> \param[out] WI
153 *> \verbatim
154 *> WI is DOUBLE PRECISION array, dimension (N)
155 *> WR and WI contain the real and imaginary parts,
156 *> respectively, of the computed eigenvalues. Complex
157 *> conjugate pairs of eigenvalues will appear consecutively
158 *> with the eigenvalue having the positive imaginary part
159 *> first.
160 *> \endverbatim
161 *>
162 *> \param[out] VL
163 *> \verbatim
164 *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
165 *> If JOBVL = 'V', the left eigenvectors u(j) are stored one
166 *> after another in the columns of VL, in the same order
167 *> as their eigenvalues.
168 *> If JOBVL = 'N', VL is not referenced.
169 *> If the j-th eigenvalue is real, then u(j) = VL(:,j),
170 *> the j-th column of VL.
171 *> If the j-th and (j+1)-st eigenvalues form a complex
172 *> conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
173 *> u(j+1) = VL(:,j) - i*VL(:,j+1).
174 *> \endverbatim
175 *>
176 *> \param[in] LDVL
177 *> \verbatim
178 *> LDVL is INTEGER
179 *> The leading dimension of the array VL. LDVL >= 1; if
180 *> JOBVL = 'V', LDVL >= N.
181 *> \endverbatim
182 *>
183 *> \param[out] VR
184 *> \verbatim
185 *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
186 *> If JOBVR = 'V', the right eigenvectors v(j) are stored one
187 *> after another in the columns of VR, in the same order
188 *> as their eigenvalues.
189 *> If JOBVR = 'N', VR is not referenced.
190 *> If the j-th eigenvalue is real, then v(j) = VR(:,j),
191 *> the j-th column of VR.
192 *> If the j-th and (j+1)-st eigenvalues form a complex
193 *> conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
194 *> v(j+1) = VR(:,j) - i*VR(:,j+1).
195 *> \endverbatim
196 *>
197 *> \param[in] LDVR
198 *> \verbatim
199 *> LDVR is INTEGER
200 *> The leading dimension of the array VR. LDVR >= 1, and if
201 *> JOBVR = 'V', LDVR >= N.
202 *> \endverbatim
203 *>
204 *> \param[out] ILO
205 *> \verbatim
206 *> ILO is INTEGER
207 *> \endverbatim
208 *>
209 *> \param[out] IHI
210 *> \verbatim
211 *> IHI is INTEGER
212 *> ILO and IHI are integer values determined when A was
213 *> balanced. The balanced A(i,j) = 0 if I > J and
214 *> J = 1,...,ILO-1 or I = IHI+1,...,N.
215 *> \endverbatim
216 *>
217 *> \param[out] SCALE
218 *> \verbatim
219 *> SCALE is DOUBLE PRECISION array, dimension (N)
220 *> Details of the permutations and scaling factors applied
221 *> when balancing A. If P(j) is the index of the row and column
222 *> interchanged with row and column j, and D(j) is the scaling
223 *> factor applied to row and column j, then
224 *> SCALE(J) = P(J), for J = 1,...,ILO-1
225 *> = D(J), for J = ILO,...,IHI
226 *> = P(J) for J = IHI+1,...,N.
227 *> The order in which the interchanges are made is N to IHI+1,
228 *> then 1 to ILO-1.
229 *> \endverbatim
230 *>
231 *> \param[out] ABNRM
232 *> \verbatim
233 *> ABNRM is DOUBLE PRECISION
234 *> The one-norm of the balanced matrix (the maximum
235 *> of the sum of absolute values of elements of any column).
236 *> \endverbatim
237 *>
238 *> \param[out] RCONDE
239 *> \verbatim
240 *> RCONDE is DOUBLE PRECISION array, dimension (N)
241 *> RCONDE(j) is the reciprocal condition number of the j-th
242 *> eigenvalue.
243 *> \endverbatim
244 *>
245 *> \param[out] RCONDV
246 *> \verbatim
247 *> RCONDV is DOUBLE PRECISION array, dimension (N)
248 *> RCONDV(j) is the reciprocal condition number of the j-th
249 *> right eigenvector.
250 *> \endverbatim
251 *>
252 *> \param[out] WORK
253 *> \verbatim
254 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
255 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
256 *> \endverbatim
257 *>
258 *> \param[in] LWORK
259 *> \verbatim
260 *> LWORK is INTEGER
261 *> The dimension of the array WORK. If SENSE = 'N' or 'E',
262 *> LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
263 *> LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6).
264 *> For good performance, LWORK must generally be larger.
265 *>
266 *> If LWORK = -1, then a workspace query is assumed; the routine
267 *> only calculates the optimal size of the WORK array, returns
268 *> this value as the first entry of the WORK array, and no error
269 *> message related to LWORK is issued by XERBLA.
270 *> \endverbatim
271 *>
272 *> \param[out] IWORK
273 *> \verbatim
274 *> IWORK is INTEGER array, dimension (2*N-2)
275 *> If SENSE = 'N' or 'E', not referenced.
276 *> \endverbatim
277 *>
278 *> \param[out] INFO
279 *> \verbatim
280 *> INFO is INTEGER
281 *> = 0: successful exit
282 *> < 0: if INFO = -i, the i-th argument had an illegal value.
283 *> > 0: if INFO = i, the QR algorithm failed to compute all the
284 *> eigenvalues, and no eigenvectors or condition numbers
285 *> have been computed; elements 1:ILO-1 and i+1:N of WR
286 *> and WI contain eigenvalues which have converged.
287 *> \endverbatim
288 *
289 * Authors:
290 * ========
291 *
292 *> \author Univ. of Tennessee
293 *> \author Univ. of California Berkeley
294 *> \author Univ. of Colorado Denver
295 *> \author NAG Ltd.
296 *
297 *> \date September 2012
298 *
299 *> \ingroup doubleGEeigen
300 *
301 * =====================================================================
302  SUBROUTINE dgeevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
303  $ vl, ldvl, vr, ldvr, ilo, ihi, scale, abnrm,
304  $ rconde, rcondv, work, lwork, iwork, info )
305 *
306 * -- LAPACK driver routine (version 3.4.2) --
307 * -- LAPACK is a software package provided by Univ. of Tennessee, --
308 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
309 * September 2012
310 *
311 * .. Scalar Arguments ..
312  CHARACTER balanc, jobvl, jobvr, sense
313  INTEGER ihi, ilo, info, lda, ldvl, ldvr, lwork, n
314  DOUBLE PRECISION abnrm
315 * ..
316 * .. Array Arguments ..
317  INTEGER iwork( * )
318  DOUBLE PRECISION a( lda, * ), rconde( * ), rcondv( * ),
319  $ scale( * ), vl( ldvl, * ), vr( ldvr, * ),
320  $ wi( * ), work( * ), wr( * )
321 * ..
322 *
323 * =====================================================================
324 *
325 * .. Parameters ..
326  DOUBLE PRECISION zero, one
327  parameter( zero = 0.0d0, one = 1.0d0 )
328 * ..
329 * .. Local Scalars ..
330  LOGICAL lquery, scalea, wantvl, wantvr, wntsnb, wntsne,
331  $ wntsnn, wntsnv
332  CHARACTER job, side
333  INTEGER hswork, i, icond, ierr, itau, iwrk, k, maxwrk,
334  $ minwrk, nout
335  DOUBLE PRECISION anrm, bignum, cs, cscale, eps, r, scl, smlnum,
336  $ sn
337 * ..
338 * .. Local Arrays ..
339  LOGICAL select( 1 )
340  DOUBLE PRECISION dum( 1 )
341 * ..
342 * .. External Subroutines ..
343  EXTERNAL dgebak, dgebal, dgehrd, dhseqr, dlabad, dlacpy,
345  $ dtrsna, xerbla
346 * ..
347 * .. External Functions ..
348  LOGICAL lsame
349  INTEGER idamax, ilaenv
350  DOUBLE PRECISION dlamch, dlange, dlapy2, dnrm2
351  EXTERNAL lsame, idamax, ilaenv, dlamch, dlange, dlapy2,
352  $ dnrm2
353 * ..
354 * .. Intrinsic Functions ..
355  INTRINSIC max, sqrt
356 * ..
357 * .. Executable Statements ..
358 *
359 * Test the input arguments
360 *
361  info = 0
362  lquery = ( lwork.EQ.-1 )
363  wantvl = lsame( jobvl, 'V' )
364  wantvr = lsame( jobvr, 'V' )
365  wntsnn = lsame( sense, 'N' )
366  wntsne = lsame( sense, 'E' )
367  wntsnv = lsame( sense, 'V' )
368  wntsnb = lsame( sense, 'B' )
369  IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc,
370  $ 'S' ) .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
371  $ THEN
372  info = -1
373  ELSE IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
374  info = -2
375  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
376  info = -3
377  ELSE IF( .NOT.( wntsnn .OR. wntsne .OR. wntsnb .OR. wntsnv ) .OR.
378  $ ( ( wntsne .OR. wntsnb ) .AND. .NOT.( wantvl .AND.
379  $ wantvr ) ) ) THEN
380  info = -4
381  ELSE IF( n.LT.0 ) THEN
382  info = -5
383  ELSE IF( lda.LT.max( 1, n ) ) THEN
384  info = -7
385  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
386  info = -11
387  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
388  info = -13
389  END IF
390 *
391 * Compute workspace
392 * (Note: Comments in the code beginning "Workspace:" describe the
393 * minimal amount of workspace needed at that point in the code,
394 * as well as the preferred amount for good performance.
395 * NB refers to the optimal block size for the immediately
396 * following subroutine, as returned by ILAENV.
397 * HSWORK refers to the workspace preferred by DHSEQR, as
398 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
399 * the worst case.)
400 *
401  IF( info.EQ.0 ) THEN
402  IF( n.EQ.0 ) THEN
403  minwrk = 1
404  maxwrk = 1
405  ELSE
406  maxwrk = n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
407 *
408  IF( wantvl ) THEN
409  CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
410  $ work, -1, info )
411  ELSE IF( wantvr ) THEN
412  CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
413  $ work, -1, info )
414  ELSE
415  IF( wntsnn ) THEN
416  CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr,
417  $ ldvr, work, -1, info )
418  ELSE
419  CALL dhseqr( 'S', 'N', n, 1, n, a, lda, wr, wi, vr,
420  $ ldvr, work, -1, info )
421  END IF
422  END IF
423  hswork = work( 1 )
424 *
425  IF( ( .NOT.wantvl ) .AND. ( .NOT.wantvr ) ) THEN
426  minwrk = 2*n
427  IF( .NOT.wntsnn )
428  $ minwrk = max( minwrk, n*n+6*n )
429  maxwrk = max( maxwrk, hswork )
430  IF( .NOT.wntsnn )
431  $ maxwrk = max( maxwrk, n*n + 6*n )
432  ELSE
433  minwrk = 3*n
434  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
435  $ minwrk = max( minwrk, n*n + 6*n )
436  maxwrk = max( maxwrk, hswork )
437  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'DORGHR',
438  $ ' ', n, 1, n, -1 ) )
439  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
440  $ maxwrk = max( maxwrk, n*n + 6*n )
441  maxwrk = max( maxwrk, 3*n )
442  END IF
443  maxwrk = max( maxwrk, minwrk )
444  END IF
445  work( 1 ) = maxwrk
446 *
447  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
448  info = -21
449  END IF
450  END IF
451 *
452  IF( info.NE.0 ) THEN
453  CALL xerbla( 'DGEEVX', -info )
454  RETURN
455  ELSE IF( lquery ) THEN
456  RETURN
457  END IF
458 *
459 * Quick return if possible
460 *
461  IF( n.EQ.0 )
462  $ RETURN
463 *
464 * Get machine constants
465 *
466  eps = dlamch( 'P' )
467  smlnum = dlamch( 'S' )
468  bignum = one / smlnum
469  CALL dlabad( smlnum, bignum )
470  smlnum = sqrt( smlnum ) / eps
471  bignum = one / smlnum
472 *
473 * Scale A if max element outside range [SMLNUM,BIGNUM]
474 *
475  icond = 0
476  anrm = dlange( 'M', n, n, a, lda, dum )
477  scalea = .false.
478  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
479  scalea = .true.
480  cscale = smlnum
481  ELSE IF( anrm.GT.bignum ) THEN
482  scalea = .true.
483  cscale = bignum
484  END IF
485  IF( scalea )
486  $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
487 *
488 * Balance the matrix and compute ABNRM
489 *
490  CALL dgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
491  abnrm = dlange( '1', n, n, a, lda, dum )
492  IF( scalea ) THEN
493  dum( 1 ) = abnrm
494  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
495  abnrm = dum( 1 )
496  END IF
497 *
498 * Reduce to upper Hessenberg form
499 * (Workspace: need 2*N, prefer N+N*NB)
500 *
501  itau = 1
502  iwrk = itau + n
503  CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
504  $ lwork-iwrk+1, ierr )
505 *
506  IF( wantvl ) THEN
507 *
508 * Want left eigenvectors
509 * Copy Householder vectors to VL
510 *
511  side = 'L'
512  CALL dlacpy( 'L', n, n, a, lda, vl, ldvl )
513 *
514 * Generate orthogonal matrix in VL
515 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
516 *
517  CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
518  $ lwork-iwrk+1, ierr )
519 *
520 * Perform QR iteration, accumulating Schur vectors in VL
521 * (Workspace: need 1, prefer HSWORK (see comments) )
522 *
523  iwrk = itau
524  CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
525  $ work( iwrk ), lwork-iwrk+1, info )
526 *
527  IF( wantvr ) THEN
528 *
529 * Want left and right eigenvectors
530 * Copy Schur vectors to VR
531 *
532  side = 'B'
533  CALL dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
534  END IF
535 *
536  ELSE IF( wantvr ) THEN
537 *
538 * Want right eigenvectors
539 * Copy Householder vectors to VR
540 *
541  side = 'R'
542  CALL dlacpy( 'L', n, n, a, lda, vr, ldvr )
543 *
544 * Generate orthogonal matrix in VR
545 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
546 *
547  CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
548  $ lwork-iwrk+1, ierr )
549 *
550 * Perform QR iteration, accumulating Schur vectors in VR
551 * (Workspace: need 1, prefer HSWORK (see comments) )
552 *
553  iwrk = itau
554  CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
555  $ work( iwrk ), lwork-iwrk+1, info )
556 *
557  ELSE
558 *
559 * Compute eigenvalues only
560 * If condition numbers desired, compute Schur form
561 *
562  IF( wntsnn ) THEN
563  job = 'E'
564  ELSE
565  job = 'S'
566  END IF
567 *
568 * (Workspace: need 1, prefer HSWORK (see comments) )
569 *
570  iwrk = itau
571  CALL dhseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
572  $ work( iwrk ), lwork-iwrk+1, info )
573  END IF
574 *
575 * If INFO > 0 from DHSEQR, then quit
576 *
577  IF( info.GT.0 )
578  $ go to 50
579 *
580  IF( wantvl .OR. wantvr ) THEN
581 *
582 * Compute left and/or right eigenvectors
583 * (Workspace: need 3*N)
584 *
585  CALL dtrevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
586  $ n, nout, work( iwrk ), ierr )
587  END IF
588 *
589 * Compute condition numbers if desired
590 * (Workspace: need N*N+6*N unless SENSE = 'E')
591 *
592  IF( .NOT.wntsnn ) THEN
593  CALL dtrsna( sense, 'A', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
594  $ rconde, rcondv, n, nout, work( iwrk ), n, iwork,
595  $ icond )
596  END IF
597 *
598  IF( wantvl ) THEN
599 *
600 * Undo balancing of left eigenvectors
601 *
602  CALL dgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,
603  $ ierr )
604 *
605 * Normalize left eigenvectors and make largest component real
606 *
607  DO 20 i = 1, n
608  IF( wi( i ).EQ.zero ) THEN
609  scl = one / dnrm2( n, vl( 1, i ), 1 )
610  CALL dscal( n, scl, vl( 1, i ), 1 )
611  ELSE IF( wi( i ).GT.zero ) THEN
612  scl = one / dlapy2( dnrm2( n, vl( 1, i ), 1 ),
613  $ dnrm2( n, vl( 1, i+1 ), 1 ) )
614  CALL dscal( n, scl, vl( 1, i ), 1 )
615  CALL dscal( n, scl, vl( 1, i+1 ), 1 )
616  DO 10 k = 1, n
617  work( k ) = vl( k, i )**2 + vl( k, i+1 )**2
618  10 CONTINUE
619  k = idamax( n, work, 1 )
620  CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
621  CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
622  vl( k, i+1 ) = zero
623  END IF
624  20 CONTINUE
625  END IF
626 *
627  IF( wantvr ) THEN
628 *
629 * Undo balancing of right eigenvectors
630 *
631  CALL dgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,
632  $ ierr )
633 *
634 * Normalize right eigenvectors and make largest component real
635 *
636  DO 40 i = 1, n
637  IF( wi( i ).EQ.zero ) THEN
638  scl = one / dnrm2( n, vr( 1, i ), 1 )
639  CALL dscal( n, scl, vr( 1, i ), 1 )
640  ELSE IF( wi( i ).GT.zero ) THEN
641  scl = one / dlapy2( dnrm2( n, vr( 1, i ), 1 ),
642  $ dnrm2( n, vr( 1, i+1 ), 1 ) )
643  CALL dscal( n, scl, vr( 1, i ), 1 )
644  CALL dscal( n, scl, vr( 1, i+1 ), 1 )
645  DO 30 k = 1, n
646  work( k ) = vr( k, i )**2 + vr( k, i+1 )**2
647  30 CONTINUE
648  k = idamax( n, work, 1 )
649  CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
650  CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
651  vr( k, i+1 ) = zero
652  END IF
653  40 CONTINUE
654  END IF
655 *
656 * Undo scaling if necessary
657 *
658  50 CONTINUE
659  IF( scalea ) THEN
660  CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
661  $ max( n-info, 1 ), ierr )
662  CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
663  $ max( n-info, 1 ), ierr )
664  IF( info.EQ.0 ) THEN
665  IF( ( wntsnv .OR. wntsnb ) .AND. icond.EQ.0 )
666  $ CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
667  $ ierr )
668  ELSE
669  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
670  $ ierr )
671  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
672  $ ierr )
673  END IF
674  END IF
675 *
676  work( 1 ) = maxwrk
677  RETURN
678 *
679 * End of DGEEVX
680 *
681  END