LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
sgeevx.f
Go to the documentation of this file.
1 *> \brief <b> SGEEVX 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 SGEEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGEEVX( 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 * REAL ABNRM
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IWORK( * )
32 * REAL 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 *> SGEEVX 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 REAL 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 REAL array, dimension (N)
150 *> \endverbatim
151 *>
152 *> \param[out] WI
153 *> \verbatim
154 *> WI is REAL 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 REAL 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 REAL 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 REAL 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 REAL
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 REAL 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 REAL 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 REAL 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 realGEeigen
300 *
301 * =====================================================================
302  SUBROUTINE sgeevx( 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  REAL ABNRM
315 * ..
316 * .. Array Arguments ..
317  INTEGER IWORK( * )
318  REAL A( lda, * ), RCONDE( * ), RCONDV( * ),
319  $ scale( * ), vl( ldvl, * ), vr( ldvr, * ),
320  $ wi( * ), work( * ), wr( * )
321 * ..
322 *
323 * =====================================================================
324 *
325 * .. Parameters ..
326  REAL ZERO, ONE
327  parameter( zero = 0.0e0, one = 1.0e0 )
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  REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
336  $ sn
337 * ..
338 * .. Local Arrays ..
339  LOGICAL SELECT( 1 )
340  REAL DUM( 1 )
341 * ..
342 * .. External Subroutines ..
343  EXTERNAL sgebak, sgebal, sgehrd, shseqr, slabad, slacpy,
345  $ strsna, xerbla
346 * ..
347 * .. External Functions ..
348  LOGICAL LSAME
349  INTEGER ILAENV, ISAMAX
350  REAL SLAMCH, SLANGE, SLAPY2, SNRM2
351  EXTERNAL lsame, ilaenv, isamax, slamch, slange, slapy2,
352  $ snrm2
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, 'S' ) .OR.
370  $ lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) ) THEN
371  info = -1
372  ELSE IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
373  info = -2
374  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
375  info = -3
376  ELSE IF( .NOT.( wntsnn .OR. wntsne .OR. wntsnb .OR. wntsnv ) .OR.
377  $ ( ( wntsne .OR. wntsnb ) .AND. .NOT.( wantvl .AND.
378  $ wantvr ) ) ) THEN
379  info = -4
380  ELSE IF( n.LT.0 ) THEN
381  info = -5
382  ELSE IF( lda.LT.max( 1, n ) ) THEN
383  info = -7
384  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
385  info = -11
386  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
387  info = -13
388  END IF
389 *
390 * Compute workspace
391 * (Note: Comments in the code beginning "Workspace:" describe the
392 * minimal amount of workspace needed at that point in the code,
393 * as well as the preferred amount for good performance.
394 * NB refers to the optimal block size for the immediately
395 * following subroutine, as returned by ILAENV.
396 * HSWORK refers to the workspace preferred by SHSEQR, as
397 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
398 * the worst case.)
399 *
400  IF( info.EQ.0 ) THEN
401  IF( n.EQ.0 ) THEN
402  minwrk = 1
403  maxwrk = 1
404  ELSE
405  maxwrk = n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
406 *
407  IF( wantvl ) THEN
408  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
409  $ work, -1, info )
410  ELSE IF( wantvr ) THEN
411  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
412  $ work, -1, info )
413  ELSE
414  IF( wntsnn ) THEN
415  CALL shseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr,
416  $ ldvr, work, -1, info )
417  ELSE
418  CALL shseqr( 'S', 'N', n, 1, n, a, lda, wr, wi, vr,
419  $ ldvr, work, -1, info )
420  END IF
421  END IF
422  hswork = work( 1 )
423 *
424  IF( ( .NOT.wantvl ) .AND. ( .NOT.wantvr ) ) THEN
425  minwrk = 2*n
426  IF( .NOT.wntsnn )
427  $ minwrk = max( minwrk, n*n+6*n )
428  maxwrk = max( maxwrk, hswork )
429  IF( .NOT.wntsnn )
430  $ maxwrk = max( maxwrk, n*n + 6*n )
431  ELSE
432  minwrk = 3*n
433  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
434  $ minwrk = max( minwrk, n*n + 6*n )
435  maxwrk = max( maxwrk, hswork )
436  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'SORGHR',
437  $ ' ', n, 1, n, -1 ) )
438  IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
439  $ maxwrk = max( maxwrk, n*n + 6*n )
440  maxwrk = max( maxwrk, 3*n )
441  END IF
442  maxwrk = max( maxwrk, minwrk )
443  END IF
444  work( 1 ) = maxwrk
445 *
446  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
447  info = -21
448  END IF
449  END IF
450 *
451  IF( info.NE.0 ) THEN
452  CALL xerbla( 'SGEEVX', -info )
453  RETURN
454  ELSE IF( lquery ) THEN
455  RETURN
456  END IF
457 *
458 * Quick return if possible
459 *
460  IF( n.EQ.0 )
461  $ RETURN
462 *
463 * Get machine constants
464 *
465  eps = slamch( 'P' )
466  smlnum = slamch( 'S' )
467  bignum = one / smlnum
468  CALL slabad( smlnum, bignum )
469  smlnum = sqrt( smlnum ) / eps
470  bignum = one / smlnum
471 *
472 * Scale A if max element outside range [SMLNUM,BIGNUM]
473 *
474  icond = 0
475  anrm = slange( 'M', n, n, a, lda, dum )
476  scalea = .false.
477  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
478  scalea = .true.
479  cscale = smlnum
480  ELSE IF( anrm.GT.bignum ) THEN
481  scalea = .true.
482  cscale = bignum
483  END IF
484  IF( scalea )
485  $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
486 *
487 * Balance the matrix and compute ABNRM
488 *
489  CALL sgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
490  abnrm = slange( '1', n, n, a, lda, dum )
491  IF( scalea ) THEN
492  dum( 1 ) = abnrm
493  CALL slascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
494  abnrm = dum( 1 )
495  END IF
496 *
497 * Reduce to upper Hessenberg form
498 * (Workspace: need 2*N, prefer N+N*NB)
499 *
500  itau = 1
501  iwrk = itau + n
502  CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
503  $ lwork-iwrk+1, ierr )
504 *
505  IF( wantvl ) THEN
506 *
507 * Want left eigenvectors
508 * Copy Householder vectors to VL
509 *
510  side = 'L'
511  CALL slacpy( 'L', n, n, a, lda, vl, ldvl )
512 *
513 * Generate orthogonal matrix in VL
514 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
515 *
516  CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
517  $ lwork-iwrk+1, ierr )
518 *
519 * Perform QR iteration, accumulating Schur vectors in VL
520 * (Workspace: need 1, prefer HSWORK (see comments) )
521 *
522  iwrk = itau
523  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
524  $ work( iwrk ), lwork-iwrk+1, info )
525 *
526  IF( wantvr ) THEN
527 *
528 * Want left and right eigenvectors
529 * Copy Schur vectors to VR
530 *
531  side = 'B'
532  CALL slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
533  END IF
534 *
535  ELSE IF( wantvr ) THEN
536 *
537 * Want right eigenvectors
538 * Copy Householder vectors to VR
539 *
540  side = 'R'
541  CALL slacpy( 'L', n, n, a, lda, vr, ldvr )
542 *
543 * Generate orthogonal matrix in VR
544 * (Workspace: need 2*N-1, prefer N+(N-1)*NB)
545 *
546  CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
547  $ lwork-iwrk+1, ierr )
548 *
549 * Perform QR iteration, accumulating Schur vectors in VR
550 * (Workspace: need 1, prefer HSWORK (see comments) )
551 *
552  iwrk = itau
553  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
554  $ work( iwrk ), lwork-iwrk+1, info )
555 *
556  ELSE
557 *
558 * Compute eigenvalues only
559 * If condition numbers desired, compute Schur form
560 *
561  IF( wntsnn ) THEN
562  job = 'E'
563  ELSE
564  job = 'S'
565  END IF
566 *
567 * (Workspace: need 1, prefer HSWORK (see comments) )
568 *
569  iwrk = itau
570  CALL shseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
571  $ work( iwrk ), lwork-iwrk+1, info )
572  END IF
573 *
574 * If INFO > 0 from SHSEQR, then quit
575 *
576  IF( info.GT.0 )
577  $ GO TO 50
578 *
579  IF( wantvl .OR. wantvr ) THEN
580 *
581 * Compute left and/or right eigenvectors
582 * (Workspace: need 3*N)
583 *
584  CALL strevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
585  $ n, nout, work( iwrk ), ierr )
586  END IF
587 *
588 * Compute condition numbers if desired
589 * (Workspace: need N*N+6*N unless SENSE = 'E')
590 *
591  IF( .NOT.wntsnn ) THEN
592  CALL strsna( sense, 'A', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
593  $ rconde, rcondv, n, nout, work( iwrk ), n, iwork,
594  $ icond )
595  END IF
596 *
597  IF( wantvl ) THEN
598 *
599 * Undo balancing of left eigenvectors
600 *
601  CALL sgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,
602  $ ierr )
603 *
604 * Normalize left eigenvectors and make largest component real
605 *
606  DO 20 i = 1, n
607  IF( wi( i ).EQ.zero ) THEN
608  scl = one / snrm2( n, vl( 1, i ), 1 )
609  CALL sscal( n, scl, vl( 1, i ), 1 )
610  ELSE IF( wi( i ).GT.zero ) THEN
611  scl = one / slapy2( snrm2( n, vl( 1, i ), 1 ),
612  $ snrm2( n, vl( 1, i+1 ), 1 ) )
613  CALL sscal( n, scl, vl( 1, i ), 1 )
614  CALL sscal( n, scl, vl( 1, i+1 ), 1 )
615  DO 10 k = 1, n
616  work( k ) = vl( k, i )**2 + vl( k, i+1 )**2
617  10 CONTINUE
618  k = isamax( n, work, 1 )
619  CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
620  CALL srot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
621  vl( k, i+1 ) = zero
622  END IF
623  20 CONTINUE
624  END IF
625 *
626  IF( wantvr ) THEN
627 *
628 * Undo balancing of right eigenvectors
629 *
630  CALL sgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,
631  $ ierr )
632 *
633 * Normalize right eigenvectors and make largest component real
634 *
635  DO 40 i = 1, n
636  IF( wi( i ).EQ.zero ) THEN
637  scl = one / snrm2( n, vr( 1, i ), 1 )
638  CALL sscal( n, scl, vr( 1, i ), 1 )
639  ELSE IF( wi( i ).GT.zero ) THEN
640  scl = one / slapy2( snrm2( n, vr( 1, i ), 1 ),
641  $ snrm2( n, vr( 1, i+1 ), 1 ) )
642  CALL sscal( n, scl, vr( 1, i ), 1 )
643  CALL sscal( n, scl, vr( 1, i+1 ), 1 )
644  DO 30 k = 1, n
645  work( k ) = vr( k, i )**2 + vr( k, i+1 )**2
646  30 CONTINUE
647  k = isamax( n, work, 1 )
648  CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
649  CALL srot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
650  vr( k, i+1 ) = zero
651  END IF
652  40 CONTINUE
653  END IF
654 *
655 * Undo scaling if necessary
656 *
657  50 CONTINUE
658  IF( scalea ) THEN
659  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
660  $ max( n-info, 1 ), ierr )
661  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
662  $ max( n-info, 1 ), ierr )
663  IF( info.EQ.0 ) THEN
664  IF( ( wntsnv .OR. wntsnb ) .AND. icond.EQ.0 )
665  $ CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
666  $ ierr )
667  ELSE
668  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
669  $ ierr )
670  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
671  $ ierr )
672  END IF
673  END IF
674 *
675  work( 1 ) = maxwrk
676  RETURN
677 *
678 * End of SGEEVX
679 *
680  END
subroutine strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
Definition: strsna.f:267
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:132
subroutine sgeevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: sgeevx.f:305
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:162
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:141
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
Definition: strevc.f:224