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