LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:318
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dtrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTREVC
Definition: dtrevc.f:224
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:162
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dtrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
DTRSNA
Definition: dtrsna.f:267
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:141
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:169
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:132
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:128
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dgeevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO)
DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: dgeevx.f:305