LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sggevx.f
Go to the documentation of this file.
1*> \brief <b> SGGEVX 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 SGGEVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggevx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggevx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggevx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
22* ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
23* IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
24* RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
25*
26* .. Scalar Arguments ..
27* CHARACTER BALANC, JOBVL, JOBVR, SENSE
28* INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
29* REAL ABNRM, BBNRM
30* ..
31* .. Array Arguments ..
32* LOGICAL BWORK( * )
33* INTEGER IWORK( * )
34* REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
35* $ B( LDB, * ), BETA( * ), LSCALE( * ),
36* $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
37* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
38* ..
39*
40*
41*> \par Purpose:
42* =============
43*>
44*> \verbatim
45*>
46*> SGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
47*> the generalized eigenvalues, and optionally, the left and/or right
48*> generalized eigenvectors.
49*>
50*> Optionally also, it computes a balancing transformation to improve
51*> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
52*> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
53*> the eigenvalues (RCONDE), and reciprocal condition numbers for the
54*> right eigenvectors (RCONDV).
55*>
56*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
57*> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
58*> singular. It is usually represented as the pair (alpha,beta), as
59*> there is a reasonable interpretation for beta=0, and even for both
60*> being zero.
61*>
62*> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
63*> of (A,B) satisfies
64*>
65*> A * v(j) = lambda(j) * B * v(j) .
66*>
67*> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
68*> of (A,B) satisfies
69*>
70*> u(j)**H * A = lambda(j) * u(j)**H * B.
71*>
72*> where u(j)**H is the conjugate-transpose of u(j).
73*>
74*> \endverbatim
75*
76* Arguments:
77* ==========
78*
79*> \param[in] BALANC
80*> \verbatim
81*> BALANC is CHARACTER*1
82*> Specifies the balance option to be performed.
83*> = 'N': do not diagonally scale or permute;
84*> = 'P': permute only;
85*> = 'S': scale only;
86*> = 'B': both permute and scale.
87*> Computed reciprocal condition numbers will be for the
88*> matrices after permuting and/or balancing. Permuting does
89*> not change condition numbers (in exact arithmetic), but
90*> balancing does.
91*> \endverbatim
92*>
93*> \param[in] JOBVL
94*> \verbatim
95*> JOBVL is CHARACTER*1
96*> = 'N': do not compute the left generalized eigenvectors;
97*> = 'V': compute the left generalized eigenvectors.
98*> \endverbatim
99*>
100*> \param[in] JOBVR
101*> \verbatim
102*> JOBVR is CHARACTER*1
103*> = 'N': do not compute the right generalized eigenvectors;
104*> = 'V': compute the right generalized eigenvectors.
105*> \endverbatim
106*>
107*> \param[in] SENSE
108*> \verbatim
109*> SENSE is CHARACTER*1
110*> Determines which reciprocal condition numbers are computed.
111*> = 'N': none are computed;
112*> = 'E': computed for eigenvalues only;
113*> = 'V': computed for eigenvectors only;
114*> = 'B': computed for eigenvalues and eigenvectors.
115*> \endverbatim
116*>
117*> \param[in] N
118*> \verbatim
119*> N is INTEGER
120*> The order of the matrices A, B, VL, and VR. N >= 0.
121*> \endverbatim
122*>
123*> \param[in,out] A
124*> \verbatim
125*> A is REAL array, dimension (LDA, N)
126*> On entry, the matrix A in the pair (A,B).
127*> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
128*> or both, then A contains the first part of the real Schur
129*> form of the "balanced" versions of the input A and B.
130*> \endverbatim
131*>
132*> \param[in] LDA
133*> \verbatim
134*> LDA is INTEGER
135*> The leading dimension of A. LDA >= max(1,N).
136*> \endverbatim
137*>
138*> \param[in,out] B
139*> \verbatim
140*> B is REAL array, dimension (LDB, N)
141*> On entry, the matrix B in the pair (A,B).
142*> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
143*> or both, then B contains the second part of the real Schur
144*> form of the "balanced" versions of the input A and B.
145*> \endverbatim
146*>
147*> \param[in] LDB
148*> \verbatim
149*> LDB is INTEGER
150*> The leading dimension of B. LDB >= max(1,N).
151*> \endverbatim
152*>
153*> \param[out] ALPHAR
154*> \verbatim
155*> ALPHAR is REAL array, dimension (N)
156*> \endverbatim
157*>
158*> \param[out] ALPHAI
159*> \verbatim
160*> ALPHAI is REAL array, dimension (N)
161*> \endverbatim
162*>
163*> \param[out] BETA
164*> \verbatim
165*> BETA is REAL array, dimension (N)
166*> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
167*> be the generalized eigenvalues. If ALPHAI(j) is zero, then
168*> the j-th eigenvalue is real; if positive, then the j-th and
169*> (j+1)-st eigenvalues are a complex conjugate pair, with
170*> ALPHAI(j+1) negative.
171*>
172*> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
173*> may easily over- or underflow, and BETA(j) may even be zero.
174*> Thus, the user should avoid naively computing the ratio
175*> ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
176*> than and usually comparable with norm(A) in magnitude, and
177*> BETA always less than and usually comparable with norm(B).
178*> \endverbatim
179*>
180*> \param[out] VL
181*> \verbatim
182*> VL is REAL array, dimension (LDVL,N)
183*> If JOBVL = 'V', the left eigenvectors u(j) are stored one
184*> after another in the columns of VL, in the same order as
185*> their eigenvalues. If the j-th eigenvalue is real, then
186*> u(j) = VL(:,j), the j-th column of VL. If the j-th and
187*> (j+1)-th eigenvalues form a complex conjugate pair, then
188*> u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
189*> Each eigenvector will be scaled so the largest component have
190*> abs(real part) + abs(imag. part) = 1.
191*> Not referenced if JOBVL = 'N'.
192*> \endverbatim
193*>
194*> \param[in] LDVL
195*> \verbatim
196*> LDVL is INTEGER
197*> The leading dimension of the matrix VL. LDVL >= 1, and
198*> if JOBVL = 'V', LDVL >= N.
199*> \endverbatim
200*>
201*> \param[out] VR
202*> \verbatim
203*> VR is REAL array, dimension (LDVR,N)
204*> If JOBVR = 'V', the right eigenvectors v(j) are stored one
205*> after another in the columns of VR, in the same order as
206*> their eigenvalues. If the j-th eigenvalue is real, then
207*> v(j) = VR(:,j), the j-th column of VR. If the j-th and
208*> (j+1)-th eigenvalues form a complex conjugate pair, then
209*> v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
210*> Each eigenvector will be scaled so the largest component have
211*> abs(real part) + abs(imag. part) = 1.
212*> Not referenced if JOBVR = 'N'.
213*> \endverbatim
214*>
215*> \param[in] LDVR
216*> \verbatim
217*> LDVR is INTEGER
218*> The leading dimension of the matrix VR. LDVR >= 1, and
219*> if JOBVR = 'V', LDVR >= N.
220*> \endverbatim
221*>
222*> \param[out] ILO
223*> \verbatim
224*> ILO is INTEGER
225*> \endverbatim
226*>
227*> \param[out] IHI
228*> \verbatim
229*> IHI is INTEGER
230*> ILO and IHI are integer values such that on exit
231*> A(i,j) = 0 and B(i,j) = 0 if i > j and
232*> j = 1,...,ILO-1 or i = IHI+1,...,N.
233*> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
234*> \endverbatim
235*>
236*> \param[out] LSCALE
237*> \verbatim
238*> LSCALE is REAL array, dimension (N)
239*> Details of the permutations and scaling factors applied
240*> to the left side of A and B. If PL(j) is the index of the
241*> row interchanged with row j, and DL(j) is the scaling
242*> factor applied to row j, then
243*> LSCALE(j) = PL(j) for j = 1,...,ILO-1
244*> = DL(j) for j = ILO,...,IHI
245*> = PL(j) for j = IHI+1,...,N.
246*> The order in which the interchanges are made is N to IHI+1,
247*> then 1 to ILO-1.
248*> \endverbatim
249*>
250*> \param[out] RSCALE
251*> \verbatim
252*> RSCALE is REAL array, dimension (N)
253*> Details of the permutations and scaling factors applied
254*> to the right side of A and B. If PR(j) is the index of the
255*> column interchanged with column j, and DR(j) is the scaling
256*> factor applied to column j, then
257*> RSCALE(j) = PR(j) for j = 1,...,ILO-1
258*> = DR(j) for j = ILO,...,IHI
259*> = PR(j) for j = IHI+1,...,N
260*> The order in which the interchanges are made is N to IHI+1,
261*> then 1 to ILO-1.
262*> \endverbatim
263*>
264*> \param[out] ABNRM
265*> \verbatim
266*> ABNRM is REAL
267*> The one-norm of the balanced matrix A.
268*> \endverbatim
269*>
270*> \param[out] BBNRM
271*> \verbatim
272*> BBNRM is REAL
273*> The one-norm of the balanced matrix B.
274*> \endverbatim
275*>
276*> \param[out] RCONDE
277*> \verbatim
278*> RCONDE is REAL array, dimension (N)
279*> If SENSE = 'E' or 'B', the reciprocal condition numbers of
280*> the eigenvalues, stored in consecutive elements of the array.
281*> For a complex conjugate pair of eigenvalues two consecutive
282*> elements of RCONDE are set to the same value. Thus RCONDE(j),
283*> RCONDV(j), and the j-th columns of VL and VR all correspond
284*> to the j-th eigenpair.
285*> If SENSE = 'N' or 'V', RCONDE is not referenced.
286*> \endverbatim
287*>
288*> \param[out] RCONDV
289*> \verbatim
290*> RCONDV is REAL array, dimension (N)
291*> If SENSE = 'V' or 'B', the estimated reciprocal condition
292*> numbers of the eigenvectors, stored in consecutive elements
293*> of the array. For a complex eigenvector two consecutive
294*> elements of RCONDV are set to the same value. If the
295*> eigenvalues cannot be reordered to compute RCONDV(j),
296*> RCONDV(j) is set to 0; this can only occur when the true
297*> value would be very small anyway.
298*> If SENSE = 'N' or 'E', RCONDV is not referenced.
299*> \endverbatim
300*>
301*> \param[out] WORK
302*> \verbatim
303*> WORK is REAL array, dimension (MAX(1,LWORK))
304*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
305*> \endverbatim
306*>
307*> \param[in] LWORK
308*> \verbatim
309*> LWORK is INTEGER
310*> The dimension of the array WORK. LWORK >= max(1,2*N).
311*> If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
312*> LWORK >= max(1,6*N).
313*> If SENSE = 'E', LWORK >= max(1,10*N).
314*> If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
315*>
316*> If LWORK = -1, then a workspace query is assumed; the routine
317*> only calculates the optimal size of the WORK array, returns
318*> this value as the first entry of the WORK array, and no error
319*> message related to LWORK is issued by XERBLA.
320*> \endverbatim
321*>
322*> \param[out] IWORK
323*> \verbatim
324*> IWORK is INTEGER array, dimension (N+6)
325*> If SENSE = 'E', IWORK is not referenced.
326*> \endverbatim
327*>
328*> \param[out] BWORK
329*> \verbatim
330*> BWORK is LOGICAL array, dimension (N)
331*> If SENSE = 'N', BWORK is not referenced.
332*> \endverbatim
333*>
334*> \param[out] INFO
335*> \verbatim
336*> INFO is INTEGER
337*> = 0: successful exit
338*> < 0: if INFO = -i, the i-th argument had an illegal value.
339*> = 1,...,N:
340*> The QZ iteration failed. No eigenvectors have been
341*> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
342*> should be correct for j=INFO+1,...,N.
343*> > N: =N+1: other than QZ iteration failed in SHGEQZ.
344*> =N+2: error return from STGEVC.
345*> \endverbatim
346*
347* Authors:
348* ========
349*
350*> \author Univ. of Tennessee
351*> \author Univ. of California Berkeley
352*> \author Univ. of Colorado Denver
353*> \author NAG Ltd.
354*
355*> \ingroup realGEeigen
356*
357*> \par Further Details:
358* =====================
359*>
360*> \verbatim
361*>
362*> Balancing a matrix pair (A,B) includes, first, permuting rows and
363*> columns to isolate eigenvalues, second, applying diagonal similarity
364*> transformation to the rows and columns to make the rows and columns
365*> as close in norm as possible. The computed reciprocal condition
366*> numbers correspond to the balanced matrix. Permuting rows and columns
367*> will not change the condition numbers (in exact arithmetic) but
368*> diagonal scaling will. For further explanation of balancing, see
369*> section 4.11.1.2 of LAPACK Users' Guide.
370*>
371*> An approximate error bound on the chordal distance between the i-th
372*> computed generalized eigenvalue w and the corresponding exact
373*> eigenvalue lambda is
374*>
375*> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
376*>
377*> An approximate error bound for the angle between the i-th computed
378*> eigenvector VL(i) or VR(i) is given by
379*>
380*> EPS * norm(ABNRM, BBNRM) / DIF(i).
381*>
382*> For further explanation of the reciprocal condition numbers RCONDE
383*> and RCONDV, see section 4.11 of LAPACK User's Guide.
384*> \endverbatim
385*>
386* =====================================================================
387 SUBROUTINE sggevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
388 $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
389 $ IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
390 $ RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
391*
392* -- LAPACK driver routine --
393* -- LAPACK is a software package provided by Univ. of Tennessee, --
394* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
395*
396* .. Scalar Arguments ..
397 CHARACTER BALANC, JOBVL, JOBVR, SENSE
398 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
399 REAL ABNRM, BBNRM
400* ..
401* .. Array Arguments ..
402 LOGICAL BWORK( * )
403 INTEGER IWORK( * )
404 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
405 $ b( ldb, * ), beta( * ), lscale( * ),
406 $ rconde( * ), rcondv( * ), rscale( * ),
407 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
408* ..
409*
410* =====================================================================
411*
412* .. Parameters ..
413 REAL ZERO, ONE
414 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
415* ..
416* .. Local Scalars ..
417 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
418 $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
419 CHARACTER CHTEMP
420 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
421 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
422 $ minwrk, mm
423 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
424 $ SMLNUM, TEMP
425* ..
426* .. Local Arrays ..
427 LOGICAL LDUMMA( 1 )
428* ..
429* .. External Subroutines ..
430 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
432 $ stgsna, xerbla
433* ..
434* .. External Functions ..
435 LOGICAL LSAME
436 INTEGER ILAENV
437 REAL SLAMCH, SLANGE
438 EXTERNAL lsame, ilaenv, slamch, slange
439* ..
440* .. Intrinsic Functions ..
441 INTRINSIC abs, max, sqrt
442* ..
443* .. Executable Statements ..
444*
445* Decode the input arguments
446*
447 IF( lsame( jobvl, 'N' ) ) THEN
448 ijobvl = 1
449 ilvl = .false.
450 ELSE IF( lsame( jobvl, 'V' ) ) THEN
451 ijobvl = 2
452 ilvl = .true.
453 ELSE
454 ijobvl = -1
455 ilvl = .false.
456 END IF
457*
458 IF( lsame( jobvr, 'N' ) ) THEN
459 ijobvr = 1
460 ilvr = .false.
461 ELSE IF( lsame( jobvr, 'V' ) ) THEN
462 ijobvr = 2
463 ilvr = .true.
464 ELSE
465 ijobvr = -1
466 ilvr = .false.
467 END IF
468 ilv = ilvl .OR. ilvr
469*
470 noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
471 wantsn = lsame( sense, 'N' )
472 wantse = lsame( sense, 'E' )
473 wantsv = lsame( sense, 'V' )
474 wantsb = lsame( sense, 'B' )
475*
476* Test the input arguments
477*
478 info = 0
479 lquery = ( lwork.EQ.-1 )
480 IF( .NOT.( noscl .OR. lsame( balanc, 'S' ) .OR.
481 $ lsame( balanc, 'B' ) ) ) THEN
482 info = -1
483 ELSE IF( ijobvl.LE.0 ) THEN
484 info = -2
485 ELSE IF( ijobvr.LE.0 ) THEN
486 info = -3
487 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
488 $ THEN
489 info = -4
490 ELSE IF( n.LT.0 ) THEN
491 info = -5
492 ELSE IF( lda.LT.max( 1, n ) ) THEN
493 info = -7
494 ELSE IF( ldb.LT.max( 1, n ) ) THEN
495 info = -9
496 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
497 info = -14
498 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
499 info = -16
500 END IF
501*
502* Compute workspace
503* (Note: Comments in the code beginning "Workspace:" describe the
504* minimal amount of workspace needed at that point in the code,
505* as well as the preferred amount for good performance.
506* NB refers to the optimal block size for the immediately
507* following subroutine, as returned by ILAENV. The workspace is
508* computed assuming ILO = 1 and IHI = N, the worst case.)
509*
510 IF( info.EQ.0 ) THEN
511 IF( n.EQ.0 ) THEN
512 minwrk = 1
513 maxwrk = 1
514 ELSE
515 IF( noscl .AND. .NOT.ilv ) THEN
516 minwrk = 2*n
517 ELSE
518 minwrk = 6*n
519 END IF
520 IF( wantse ) THEN
521 minwrk = 10*n
522 ELSE IF( wantsv .OR. wantsb ) THEN
523 minwrk = 2*n*( n + 4 ) + 16
524 END IF
525 maxwrk = minwrk
526 maxwrk = max( maxwrk,
527 $ n + n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) )
528 maxwrk = max( maxwrk,
529 $ n + n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) )
530 IF( ilvl ) THEN
531 maxwrk = max( maxwrk, n +
532 $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, 0 ) )
533 END IF
534 END IF
535 work( 1 ) = maxwrk
536*
537 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
538 info = -26
539 END IF
540 END IF
541*
542 IF( info.NE.0 ) THEN
543 CALL xerbla( 'SGGEVX', -info )
544 RETURN
545 ELSE IF( lquery ) THEN
546 RETURN
547 END IF
548*
549* Quick return if possible
550*
551 IF( n.EQ.0 )
552 $ RETURN
553*
554*
555* Get machine constants
556*
557 eps = slamch( 'P' )
558 smlnum = slamch( 'S' )
559 bignum = one / smlnum
560 CALL slabad( smlnum, bignum )
561 smlnum = sqrt( smlnum ) / eps
562 bignum = one / smlnum
563*
564* Scale A if max element outside range [SMLNUM,BIGNUM]
565*
566 anrm = slange( 'M', n, n, a, lda, work )
567 ilascl = .false.
568 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
569 anrmto = smlnum
570 ilascl = .true.
571 ELSE IF( anrm.GT.bignum ) THEN
572 anrmto = bignum
573 ilascl = .true.
574 END IF
575 IF( ilascl )
576 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
577*
578* Scale B if max element outside range [SMLNUM,BIGNUM]
579*
580 bnrm = slange( 'M', n, n, b, ldb, work )
581 ilbscl = .false.
582 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
583 bnrmto = smlnum
584 ilbscl = .true.
585 ELSE IF( bnrm.GT.bignum ) THEN
586 bnrmto = bignum
587 ilbscl = .true.
588 END IF
589 IF( ilbscl )
590 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
591*
592* Permute and/or balance the matrix pair (A,B)
593* (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
594*
595 CALL sggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
596 $ work, ierr )
597*
598* Compute ABNRM and BBNRM
599*
600 abnrm = slange( '1', n, n, a, lda, work( 1 ) )
601 IF( ilascl ) THEN
602 work( 1 ) = abnrm
603 CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
604 $ ierr )
605 abnrm = work( 1 )
606 END IF
607*
608 bbnrm = slange( '1', n, n, b, ldb, work( 1 ) )
609 IF( ilbscl ) THEN
610 work( 1 ) = bbnrm
611 CALL slascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
612 $ ierr )
613 bbnrm = work( 1 )
614 END IF
615*
616* Reduce B to triangular form (QR decomposition of B)
617* (Workspace: need N, prefer N*NB )
618*
619 irows = ihi + 1 - ilo
620 IF( ilv .OR. .NOT.wantsn ) THEN
621 icols = n + 1 - ilo
622 ELSE
623 icols = irows
624 END IF
625 itau = 1
626 iwrk = itau + irows
627 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
628 $ work( iwrk ), lwork+1-iwrk, ierr )
629*
630* Apply the orthogonal transformation to A
631* (Workspace: need N, prefer N*NB)
632*
633 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
634 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
635 $ lwork+1-iwrk, ierr )
636*
637* Initialize VL and/or VR
638* (Workspace: need N, prefer N*NB)
639*
640 IF( ilvl ) THEN
641 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
642 IF( irows.GT.1 ) THEN
643 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
644 $ vl( ilo+1, ilo ), ldvl )
645 END IF
646 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
647 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
648 END IF
649*
650 IF( ilvr )
651 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
652*
653* Reduce to generalized Hessenberg form
654* (Workspace: none needed)
655*
656 IF( ilv .OR. .NOT.wantsn ) THEN
657*
658* Eigenvectors requested -- work on whole matrix.
659*
660 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
661 $ ldvl, vr, ldvr, ierr )
662 ELSE
663 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
664 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
665 END IF
666*
667* Perform QZ algorithm (Compute eigenvalues, and optionally, the
668* Schur forms and Schur vectors)
669* (Workspace: need N)
670*
671 IF( ilv .OR. .NOT.wantsn ) THEN
672 chtemp = 'S'
673 ELSE
674 chtemp = 'E'
675 END IF
676*
677 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
678 $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
679 $ lwork, ierr )
680 IF( ierr.NE.0 ) THEN
681 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
682 info = ierr
683 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
684 info = ierr - n
685 ELSE
686 info = n + 1
687 END IF
688 GO TO 130
689 END IF
690*
691* Compute Eigenvectors and estimate condition numbers if desired
692* (Workspace: STGEVC: need 6*N
693* STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
694* need N otherwise )
695*
696 IF( ilv .OR. .NOT.wantsn ) THEN
697 IF( ilv ) THEN
698 IF( ilvl ) THEN
699 IF( ilvr ) THEN
700 chtemp = 'B'
701 ELSE
702 chtemp = 'L'
703 END IF
704 ELSE
705 chtemp = 'R'
706 END IF
707*
708 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
709 $ ldvl, vr, ldvr, n, in, work, ierr )
710 IF( ierr.NE.0 ) THEN
711 info = n + 2
712 GO TO 130
713 END IF
714 END IF
715*
716 IF( .NOT.wantsn ) THEN
717*
718* compute eigenvectors (STGEVC) and estimate condition
719* numbers (STGSNA). Note that the definition of the condition
720* number is not invariant under transformation (u,v) to
721* (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
722* Schur form (S,T), Q and Z are orthogonal matrices. In order
723* to avoid using extra 2*N*N workspace, we have to recalculate
724* eigenvectors and estimate one condition numbers at a time.
725*
726 pair = .false.
727 DO 20 i = 1, n
728*
729 IF( pair ) THEN
730 pair = .false.
731 GO TO 20
732 END IF
733 mm = 1
734 IF( i.LT.n ) THEN
735 IF( a( i+1, i ).NE.zero ) THEN
736 pair = .true.
737 mm = 2
738 END IF
739 END IF
740*
741 DO 10 j = 1, n
742 bwork( j ) = .false.
743 10 CONTINUE
744 IF( mm.EQ.1 ) THEN
745 bwork( i ) = .true.
746 ELSE IF( mm.EQ.2 ) THEN
747 bwork( i ) = .true.
748 bwork( i+1 ) = .true.
749 END IF
750*
751 iwrk = mm*n + 1
752 iwrk1 = iwrk + mm*n
753*
754* Compute a pair of left and right eigenvectors.
755* (compute workspace: need up to 4*N + 6*N)
756*
757 IF( wantse .OR. wantsb ) THEN
758 CALL stgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
759 $ work( 1 ), n, work( iwrk ), n, mm, m,
760 $ work( iwrk1 ), ierr )
761 IF( ierr.NE.0 ) THEN
762 info = n + 2
763 GO TO 130
764 END IF
765 END IF
766*
767 CALL stgsna( sense, 'S', bwork, n, a, lda, b, ldb,
768 $ work( 1 ), n, work( iwrk ), n, rconde( i ),
769 $ rcondv( i ), mm, m, work( iwrk1 ),
770 $ lwork-iwrk1+1, iwork, ierr )
771*
772 20 CONTINUE
773 END IF
774 END IF
775*
776* Undo balancing on VL and VR and normalization
777* (Workspace: none needed)
778*
779 IF( ilvl ) THEN
780 CALL sggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
781 $ ldvl, ierr )
782*
783 DO 70 jc = 1, n
784 IF( alphai( jc ).LT.zero )
785 $ GO TO 70
786 temp = zero
787 IF( alphai( jc ).EQ.zero ) THEN
788 DO 30 jr = 1, n
789 temp = max( temp, abs( vl( jr, jc ) ) )
790 30 CONTINUE
791 ELSE
792 DO 40 jr = 1, n
793 temp = max( temp, abs( vl( jr, jc ) )+
794 $ abs( vl( jr, jc+1 ) ) )
795 40 CONTINUE
796 END IF
797 IF( temp.LT.smlnum )
798 $ GO TO 70
799 temp = one / temp
800 IF( alphai( jc ).EQ.zero ) THEN
801 DO 50 jr = 1, n
802 vl( jr, jc ) = vl( jr, jc )*temp
803 50 CONTINUE
804 ELSE
805 DO 60 jr = 1, n
806 vl( jr, jc ) = vl( jr, jc )*temp
807 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
808 60 CONTINUE
809 END IF
810 70 CONTINUE
811 END IF
812 IF( ilvr ) THEN
813 CALL sggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
814 $ ldvr, ierr )
815 DO 120 jc = 1, n
816 IF( alphai( jc ).LT.zero )
817 $ GO TO 120
818 temp = zero
819 IF( alphai( jc ).EQ.zero ) THEN
820 DO 80 jr = 1, n
821 temp = max( temp, abs( vr( jr, jc ) ) )
822 80 CONTINUE
823 ELSE
824 DO 90 jr = 1, n
825 temp = max( temp, abs( vr( jr, jc ) )+
826 $ abs( vr( jr, jc+1 ) ) )
827 90 CONTINUE
828 END IF
829 IF( temp.LT.smlnum )
830 $ GO TO 120
831 temp = one / temp
832 IF( alphai( jc ).EQ.zero ) THEN
833 DO 100 jr = 1, n
834 vr( jr, jc ) = vr( jr, jc )*temp
835 100 CONTINUE
836 ELSE
837 DO 110 jr = 1, n
838 vr( jr, jc ) = vr( jr, jc )*temp
839 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
840 110 CONTINUE
841 END IF
842 120 CONTINUE
843 END IF
844*
845* Undo scaling if necessary
846*
847 130 CONTINUE
848*
849 IF( ilascl ) THEN
850 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
851 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
852 END IF
853*
854 IF( ilbscl ) THEN
855 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
856 END IF
857*
858 work( 1 ) = maxwrk
859 RETURN
860*
861* End of SGGEVX
862*
863 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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
Definition: sggbak.f:147
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
Definition: sggbal.f:177
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:295
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:146
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
Definition: shgeqz.f:304
subroutine sggevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, BWORK, INFO)
SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: sggevx.f:391
subroutine stgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
STGSNA
Definition: stgsna.f:381
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:128
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
Definition: sgghrd.f:207