LAPACK  3.9.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 June 2016
298 *
299 * @precisions fortran d -> s
300 *
301 *> \ingroup doubleGEeigen
302 *
303 * =====================================================================
304  SUBROUTINE dgeevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
305  $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
306  $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
307  implicit none
308 *
309 * -- LAPACK driver routine (version 3.7.0) --
310 * -- LAPACK is a software package provided by Univ. of Tennessee, --
311 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
312 * June 2016
313 *
314 * .. Scalar Arguments ..
315  CHARACTER BALANC, JOBVL, JOBVR, SENSE
316  INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
317  DOUBLE PRECISION ABNRM
318 * ..
319 * .. Array Arguments ..
320  INTEGER IWORK( * )
321  DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ),
322  $ scale( * ), vl( ldvl, * ), vr( ldvr, * ),
323  $ wi( * ), work( * ), wr( * )
324 * ..
325 *
326 * =====================================================================
327 *
328 * .. Parameters ..
329  DOUBLE PRECISION ZERO, ONE
330  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
331 * ..
332 * .. Local Scalars ..
333  LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
334  $ WNTSNN, WNTSNV
335  CHARACTER JOB, SIDE
336  INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
337  $ lwork_trevc, maxwrk, minwrk, nout
338  DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
339  $ SN
340 * ..
341 * .. Local Arrays ..
342  LOGICAL SELECT( 1 )
343  DOUBLE PRECISION DUM( 1 )
344 * ..
345 * .. External Subroutines ..
346  EXTERNAL dgebak, dgebal, dgehrd, dhseqr, dlabad, dlacpy,
347  $ dlartg, dlascl, dorghr, drot, dscal, dtrevc3,
348  $ dtrsna, xerbla
349 * ..
350 * .. External Functions ..
351  LOGICAL LSAME
352  INTEGER IDAMAX, ILAENV
353  DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
354  EXTERNAL lsame, idamax, ilaenv, dlamch, dlange, dlapy2,
355  $ dnrm2
356 * ..
357 * .. Intrinsic Functions ..
358  INTRINSIC max, sqrt
359 * ..
360 * .. Executable Statements ..
361 *
362 * Test the input arguments
363 *
364  info = 0
365  lquery = ( lwork.EQ.-1 )
366  wantvl = lsame( jobvl, 'V' )
367  wantvr = lsame( jobvr, 'V' )
368  wntsnn = lsame( sense, 'N' )
369  wntsne = lsame( sense, 'E' )
370  wntsnv = lsame( sense, 'V' )
371  wntsnb = lsame( sense, 'B' )
372  IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc, 'S' )
373  $ .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
374  $ THEN
375  info = -1
376  ELSE IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
377  info = -2
378  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
379  info = -3
380  ELSE IF( .NOT.( wntsnn .OR. wntsne .OR. wntsnb .OR. wntsnv ) .OR.
381  $ ( ( wntsne .OR. wntsnb ) .AND. .NOT.( wantvl .AND.
382  $ wantvr ) ) ) THEN
383  info = -4
384  ELSE IF( n.LT.0 ) THEN
385  info = -5
386  ELSE IF( lda.LT.max( 1, n ) ) THEN
387  info = -7
388  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
389  info = -11
390  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
391  info = -13
392  END IF
393 *
394 * Compute workspace
395 * (Note: Comments in the code beginning "Workspace:" describe the
396 * minimal amount of workspace needed at that point in the code,
397 * as well as the preferred amount for good performance.
398 * NB refers to the optimal block size for the immediately
399 * following subroutine, as returned by ILAENV.
400 * HSWORK refers to the workspace preferred by DHSEQR, as
401 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
402 * the worst case.)
403 *
404  IF( info.EQ.0 ) THEN
405  IF( n.EQ.0 ) THEN
406  minwrk = 1
407  maxwrk = 1
408  ELSE
409  maxwrk = n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
410 *
411  IF( wantvl ) THEN
412  CALL dtrevc3( 'L', 'B', SELECT, n, a, lda,
413  $ vl, ldvl, vr, ldvr,
414  $ n, nout, work, -1, ierr )
415  lwork_trevc = int( work(1) )
416  maxwrk = max( maxwrk, n + lwork_trevc )
417  CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
418  $ work, -1, info )
419  ELSE IF( wantvr ) THEN
420  CALL dtrevc3( 'R', 'B', SELECT, n, a, lda,
421  $ vl, ldvl, vr, ldvr,
422  $ n, nout, work, -1, ierr )
423  lwork_trevc = int( work(1) )
424  maxwrk = max( maxwrk, n + lwork_trevc )
425  CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
426  $ work, -1, info )
427  ELSE
428  IF( wntsnn ) THEN
429  CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr,
430  $ ldvr, work, -1, info )
431  ELSE
432  CALL dhseqr( 'S', 'N', n, 1, n, a, lda, wr, wi, vr,
433  $ ldvr, work, -1, info )
434  END IF
435  END IF
436  hswork = int( work(1) )
437 *
438  IF( ( .NOT.wantvl ) .AND. ( .NOT.wantvr ) ) THEN
439  minwrk = 2*n
440  IF( .NOT.wntsnn )
441  $ minwrk = max( minwrk, n*n+6*n )
442  maxwrk = max( maxwrk, hswork )
443  IF( .NOT.wntsnn )
444  $ maxwrk = max( maxwrk, n*n + 6*n )
445  ELSE
446  minwrk = 3*n
447  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
448  $ minwrk = max( minwrk, n*n + 6*n )
449  maxwrk = max( maxwrk, hswork )
450  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'DORGHR',
451  $ ' ', n, 1, n, -1 ) )
452  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
453  $ maxwrk = max( maxwrk, n*n + 6*n )
454  maxwrk = max( maxwrk, 3*n )
455  END IF
456  maxwrk = max( maxwrk, minwrk )
457  END IF
458  work( 1 ) = maxwrk
459 *
460  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
461  info = -21
462  END IF
463  END IF
464 *
465  IF( info.NE.0 ) THEN
466  CALL xerbla( 'DGEEVX', -info )
467  RETURN
468  ELSE IF( lquery ) THEN
469  RETURN
470  END IF
471 *
472 * Quick return if possible
473 *
474  IF( n.EQ.0 )
475  $ RETURN
476 *
477 * Get machine constants
478 *
479  eps = dlamch( 'P' )
480  smlnum = dlamch( 'S' )
481  bignum = one / smlnum
482  CALL dlabad( smlnum, bignum )
483  smlnum = sqrt( smlnum ) / eps
484  bignum = one / smlnum
485 *
486 * Scale A if max element outside range [SMLNUM,BIGNUM]
487 *
488  icond = 0
489  anrm = dlange( 'M', n, n, a, lda, dum )
490  scalea = .false.
491  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
492  scalea = .true.
493  cscale = smlnum
494  ELSE IF( anrm.GT.bignum ) THEN
495  scalea = .true.
496  cscale = bignum
497  END IF
498  IF( scalea )
499  $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
500 *
501 * Balance the matrix and compute ABNRM
502 *
503  CALL dgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
504  abnrm = dlange( '1', n, n, a, lda, dum )
505  IF( scalea ) THEN
506  dum( 1 ) = abnrm
507  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
508  abnrm = dum( 1 )
509  END IF
510 *
511 * Reduce to upper Hessenberg form
512 * (Workspace: need 2*N, prefer N+N*NB)
513 *
514  itau = 1
515  iwrk = itau + n
516  CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
517  $ lwork-iwrk+1, ierr )
518 *
519  IF( wantvl ) THEN
520 *
521 * Want left eigenvectors
522 * Copy Householder vectors to VL
523 *
524  side = 'L'
525  CALL dlacpy( 'L', n, n, a, lda, vl, ldvl )
526 *
527 * Generate orthogonal matrix in VL
528 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
529 *
530  CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
531  $ lwork-iwrk+1, ierr )
532 *
533 * Perform QR iteration, accumulating Schur vectors in VL
534 * (Workspace: need 1, prefer HSWORK (see comments) )
535 *
536  iwrk = itau
537  CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
538  $ work( iwrk ), lwork-iwrk+1, info )
539 *
540  IF( wantvr ) THEN
541 *
542 * Want left and right eigenvectors
543 * Copy Schur vectors to VR
544 *
545  side = 'B'
546  CALL dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
547  END IF
548 *
549  ELSE IF( wantvr ) THEN
550 *
551 * Want right eigenvectors
552 * Copy Householder vectors to VR
553 *
554  side = 'R'
555  CALL dlacpy( 'L', n, n, a, lda, vr, ldvr )
556 *
557 * Generate orthogonal matrix in VR
558 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
559 *
560  CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
561  $ lwork-iwrk+1, ierr )
562 *
563 * Perform QR iteration, accumulating Schur vectors in VR
564 * (Workspace: need 1, prefer HSWORK (see comments) )
565 *
566  iwrk = itau
567  CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
568  $ work( iwrk ), lwork-iwrk+1, info )
569 *
570  ELSE
571 *
572 * Compute eigenvalues only
573 * If condition numbers desired, compute Schur form
574 *
575  IF( wntsnn ) THEN
576  job = 'E'
577  ELSE
578  job = 'S'
579  END IF
580 *
581 * (Workspace: need 1, prefer HSWORK (see comments) )
582 *
583  iwrk = itau
584  CALL dhseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
585  $ work( iwrk ), lwork-iwrk+1, info )
586  END IF
587 *
588 * If INFO .NE. 0 from DHSEQR, then quit
589 *
590  IF( info.NE.0 )
591  $ GO TO 50
592 *
593  IF( wantvl .OR. wantvr ) THEN
594 *
595 * Compute left and/or right eigenvectors
596 * (Workspace: need 3*N, prefer N + 2*N*NB)
597 *
598  CALL dtrevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
599  $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
600  END IF
601 *
602 * Compute condition numbers if desired
603 * (Workspace: need N*N+6*N unless SENSE = 'E')
604 *
605  IF( .NOT.wntsnn ) THEN
606  CALL dtrsna( sense, 'A', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
607  $ rconde, rcondv, n, nout, work( iwrk ), n, iwork,
608  $ icond )
609  END IF
610 *
611  IF( wantvl ) THEN
612 *
613 * Undo balancing of left eigenvectors
614 *
615  CALL dgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,
616  $ ierr )
617 *
618 * Normalize left eigenvectors and make largest component real
619 *
620  DO 20 i = 1, n
621  IF( wi( i ).EQ.zero ) THEN
622  scl = one / dnrm2( n, vl( 1, i ), 1 )
623  CALL dscal( n, scl, vl( 1, i ), 1 )
624  ELSE IF( wi( i ).GT.zero ) THEN
625  scl = one / dlapy2( dnrm2( n, vl( 1, i ), 1 ),
626  $ dnrm2( n, vl( 1, i+1 ), 1 ) )
627  CALL dscal( n, scl, vl( 1, i ), 1 )
628  CALL dscal( n, scl, vl( 1, i+1 ), 1 )
629  DO 10 k = 1, n
630  work( k ) = vl( k, i )**2 + vl( k, i+1 )**2
631  10 CONTINUE
632  k = idamax( n, work, 1 )
633  CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
634  CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
635  vl( k, i+1 ) = zero
636  END IF
637  20 CONTINUE
638  END IF
639 *
640  IF( wantvr ) THEN
641 *
642 * Undo balancing of right eigenvectors
643 *
644  CALL dgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,
645  $ ierr )
646 *
647 * Normalize right eigenvectors and make largest component real
648 *
649  DO 40 i = 1, n
650  IF( wi( i ).EQ.zero ) THEN
651  scl = one / dnrm2( n, vr( 1, i ), 1 )
652  CALL dscal( n, scl, vr( 1, i ), 1 )
653  ELSE IF( wi( i ).GT.zero ) THEN
654  scl = one / dlapy2( dnrm2( n, vr( 1, i ), 1 ),
655  $ dnrm2( n, vr( 1, i+1 ), 1 ) )
656  CALL dscal( n, scl, vr( 1, i ), 1 )
657  CALL dscal( n, scl, vr( 1, i+1 ), 1 )
658  DO 30 k = 1, n
659  work( k ) = vr( k, i )**2 + vr( k, i+1 )**2
660  30 CONTINUE
661  k = idamax( n, work, 1 )
662  CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
663  CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
664  vr( k, i+1 ) = zero
665  END IF
666  40 CONTINUE
667  END IF
668 *
669 * Undo scaling if necessary
670 *
671  50 CONTINUE
672  IF( scalea ) THEN
673  CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
674  $ max( n-info, 1 ), ierr )
675  CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
676  $ max( n-info, 1 ), ierr )
677  IF( info.EQ.0 ) THEN
678  IF( ( wntsnv .OR. wntsnb ) .AND. icond.EQ.0 )
679  $ CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
680  $ ierr )
681  ELSE
682  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
683  $ ierr )
684  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
685  $ ierr )
686  END IF
687  END IF
688 *
689  work( 1 ) = maxwrk
690  RETURN
691 *
692 * End of DGEEVX
693 *
694  END
dlascl
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:145
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
dgehrd
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:169
dtrevc3
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
Definition: dtrevc3.f:241
dgebal
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:162
dhseqr
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:318
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dgeevx
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:307
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:94
dgebak
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:132
dtrsna
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
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
dorghr
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:128
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
dlacpy
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