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