LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgegv.f
Go to the documentation of this file.
1*> \brief <b> SGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B).</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGEGV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgegv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgegv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgegv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
22* BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBVL, JOBVR
26* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27* ..
28* .. Array Arguments ..
29* REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
30* $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
31* $ VR( LDVR, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> This routine is deprecated and has been replaced by routine SGGEV.
41*>
42*> SGEGV computes the eigenvalues and, optionally, the left and/or right
43*> eigenvectors of a real matrix pair (A,B).
44*> Given two square matrices A and B,
45*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
46*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
47*> that
48*>
49*> A*x = lambda*B*x.
50*>
51*> An alternate form is to find the eigenvalues mu and corresponding
52*> eigenvectors y such that
53*>
54*> mu*A*y = B*y.
55*>
56*> These two forms are equivalent with mu = 1/lambda and x = y if
57*> neither lambda nor mu is zero. In order to deal with the case that
58*> lambda or mu is zero or small, two values alpha and beta are returned
59*> for each eigenvalue, such that lambda = alpha/beta and
60*> mu = beta/alpha.
61*>
62*> The vectors x and y in the above equations are right eigenvectors of
63*> the matrix pair (A,B). Vectors u and v satisfying
64*>
65*> u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
66*>
67*> are left eigenvectors of (A,B).
68*>
69*> Note: this routine performs "full balancing" on A and B
70*> \endverbatim
71*
72* Arguments:
73* ==========
74*
75*> \param[in] JOBVL
76*> \verbatim
77*> JOBVL is CHARACTER*1
78*> = 'N': do not compute the left generalized eigenvectors;
79*> = 'V': compute the left generalized eigenvectors (returned
80*> in VL).
81*> \endverbatim
82*>
83*> \param[in] JOBVR
84*> \verbatim
85*> JOBVR is CHARACTER*1
86*> = 'N': do not compute the right generalized eigenvectors;
87*> = 'V': compute the right generalized eigenvectors (returned
88*> in VR).
89*> \endverbatim
90*>
91*> \param[in] N
92*> \verbatim
93*> N is INTEGER
94*> The order of the matrices A, B, VL, and VR. N >= 0.
95*> \endverbatim
96*>
97*> \param[in,out] A
98*> \verbatim
99*> A is REAL array, dimension (LDA, N)
100*> On entry, the matrix A.
101*> If JOBVL = 'V' or JOBVR = 'V', then on exit A
102*> contains the real Schur form of A from the generalized Schur
103*> factorization of the pair (A,B) after balancing.
104*> If no eigenvectors were computed, then only the diagonal
105*> blocks from the Schur form will be correct. See SGGHRD and
106*> SHGEQZ for details.
107*> \endverbatim
108*>
109*> \param[in] LDA
110*> \verbatim
111*> LDA is INTEGER
112*> The leading dimension of A. LDA >= max(1,N).
113*> \endverbatim
114*>
115*> \param[in,out] B
116*> \verbatim
117*> B is REAL array, dimension (LDB, N)
118*> On entry, the matrix B.
119*> If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
120*> upper triangular matrix obtained from B in the generalized
121*> Schur factorization of the pair (A,B) after balancing.
122*> If no eigenvectors were computed, then only those elements of
123*> B corresponding to the diagonal blocks from the Schur form of
124*> A will be correct. See SGGHRD and SHGEQZ for details.
125*> \endverbatim
126*>
127*> \param[in] LDB
128*> \verbatim
129*> LDB is INTEGER
130*> The leading dimension of B. LDB >= max(1,N).
131*> \endverbatim
132*>
133*> \param[out] ALPHAR
134*> \verbatim
135*> ALPHAR is REAL array, dimension (N)
136*> The real parts of each scalar alpha defining an eigenvalue of
137*> GNEP.
138*> \endverbatim
139*>
140*> \param[out] ALPHAI
141*> \verbatim
142*> ALPHAI is REAL array, dimension (N)
143*> The imaginary parts of each scalar alpha defining an
144*> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
145*> eigenvalue is real; if positive, then the j-th and
146*> (j+1)-st eigenvalues are a complex conjugate pair, with
147*> ALPHAI(j+1) = -ALPHAI(j).
148*> \endverbatim
149*>
150*> \param[out] BETA
151*> \verbatim
152*> BETA is REAL array, dimension (N)
153*> The scalars beta that define the eigenvalues of GNEP.
154*>
155*> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
156*> beta = BETA(j) represent the j-th eigenvalue of the matrix
157*> pair (A,B), in one of the forms lambda = alpha/beta or
158*> mu = beta/alpha. Since either lambda or mu may overflow,
159*> they should not, in general, be computed.
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
166*> in the columns of VL, in the same order as their eigenvalues.
167*> If the j-th eigenvalue is real, then u(j) = VL(:,j).
168*> If the j-th and (j+1)-st eigenvalues form a complex conjugate
169*> pair, then
170*> u(j) = VL(:,j) + i*VL(:,j+1)
171*> and
172*> u(j+1) = VL(:,j) - i*VL(:,j+1).
173*>
174*> Each eigenvector is scaled so that its largest component has
175*> abs(real part) + abs(imag. part) = 1, except for eigenvectors
176*> corresponding to an eigenvalue with alpha = beta = 0, which
177*> are set to zero.
178*> Not referenced if JOBVL = 'N'.
179*> \endverbatim
180*>
181*> \param[in] LDVL
182*> \verbatim
183*> LDVL is INTEGER
184*> The leading dimension of the matrix VL. LDVL >= 1, and
185*> if JOBVL = 'V', LDVL >= N.
186*> \endverbatim
187*>
188*> \param[out] VR
189*> \verbatim
190*> VR is REAL array, dimension (LDVR,N)
191*> If JOBVR = 'V', the right eigenvectors x(j) are stored
192*> in the columns of VR, in the same order as their eigenvalues.
193*> If the j-th eigenvalue is real, then x(j) = VR(:,j).
194*> If the j-th and (j+1)-st eigenvalues form a complex conjugate
195*> pair, then
196*> x(j) = VR(:,j) + i*VR(:,j+1)
197*> and
198*> x(j+1) = VR(:,j) - i*VR(:,j+1).
199*>
200*> Each eigenvector is scaled so that its largest component has
201*> abs(real part) + abs(imag. part) = 1, except for eigenvalues
202*> corresponding to an eigenvalue with alpha = beta = 0, which
203*> are set to zero.
204*> Not referenced if JOBVR = 'N'.
205*> \endverbatim
206*>
207*> \param[in] LDVR
208*> \verbatim
209*> LDVR is INTEGER
210*> The leading dimension of the matrix VR. LDVR >= 1, and
211*> if JOBVR = 'V', LDVR >= N.
212*> \endverbatim
213*>
214*> \param[out] WORK
215*> \verbatim
216*> WORK is REAL array, dimension (MAX(1,LWORK))
217*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
218*> \endverbatim
219*>
220*> \param[in] LWORK
221*> \verbatim
222*> LWORK is INTEGER
223*> The dimension of the array WORK. LWORK >= max(1,8*N).
224*> For good performance, LWORK must generally be larger.
225*> To compute the optimal value of LWORK, call ILAENV to get
226*> blocksizes (for SGEQRF, SORMQR, and SORGQR.) Then compute:
227*> NB -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR;
228*> The optimal LWORK is:
229*> 2*N + MAX( 6*N, N*(NB+1) ).
230*>
231*> If LWORK = -1, then a workspace query is assumed; the routine
232*> only calculates the optimal size of the WORK array, returns
233*> this value as the first entry of the WORK array, and no error
234*> message related to LWORK is issued by XERBLA.
235*> \endverbatim
236*>
237*> \param[out] INFO
238*> \verbatim
239*> INFO is INTEGER
240*> = 0: successful exit
241*> < 0: if INFO = -i, the i-th argument had an illegal value.
242*> = 1,...,N:
243*> The QZ iteration failed. No eigenvectors have been
244*> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
245*> should be correct for j=INFO+1,...,N.
246*> > N: errors that usually indicate LAPACK problems:
247*> =N+1: error return from SGGBAL
248*> =N+2: error return from SGEQRF
249*> =N+3: error return from SORMQR
250*> =N+4: error return from SORGQR
251*> =N+5: error return from SGGHRD
252*> =N+6: error return from SHGEQZ (other than failed
253*> iteration)
254*> =N+7: error return from STGEVC
255*> =N+8: error return from SGGBAK (computing VL)
256*> =N+9: error return from SGGBAK (computing VR)
257*> =N+10: error return from SLASCL (various calls)
258*> \endverbatim
259*
260* Authors:
261* ========
262*
263*> \author Univ. of Tennessee
264*> \author Univ. of California Berkeley
265*> \author Univ. of Colorado Denver
266*> \author NAG Ltd.
267*
268*> \ingroup realGEeigen
269*
270*> \par Further Details:
271* =====================
272*>
273*> \verbatim
274*>
275*> Balancing
276*> ---------
277*>
278*> This driver calls SGGBAL to both permute and scale rows and columns
279*> of A and B. The permutations PL and PR are chosen so that PL*A*PR
280*> and PL*B*R will be upper triangular except for the diagonal blocks
281*> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
282*> possible. The diagonal scaling matrices DL and DR are chosen so
283*> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
284*> one (except for the elements that start out zero.)
285*>
286*> After the eigenvalues and eigenvectors of the balanced matrices
287*> have been computed, SGGBAK transforms the eigenvectors back to what
288*> they would have been (in perfect arithmetic) if they had not been
289*> balanced.
290*>
291*> Contents of A and B on Exit
292*> -------- -- - --- - -- ----
293*>
294*> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
295*> both), then on exit the arrays A and B will contain the real Schur
296*> form[*] of the "balanced" versions of A and B. If no eigenvectors
297*> are computed, then only the diagonal blocks will be correct.
298*>
299*> [*] See SHGEQZ, SGEGS, or read the book "Matrix Computations",
300*> by Golub & van Loan, pub. by Johns Hopkins U. Press.
301*> \endverbatim
302*>
303* =====================================================================
304 SUBROUTINE sgegv( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
305 $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
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 JOBVL, JOBVR
313 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314* ..
315* .. Array Arguments ..
316 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317 $ b( ldb, * ), beta( * ), vl( ldvl, * ),
318 $ vr( ldvr, * ), work( * )
319* ..
320*
321* =====================================================================
322*
323* .. Parameters ..
324 REAL ZERO, ONE
325 parameter( zero = 0.0e0, one = 1.0e0 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329 CHARACTER CHTEMP
330 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331 $ in, iright, irows, itau, iwork, jc, jr, lopt,
332 $ lwkmin, lwkopt, nb, nb1, nb2, nb3
333 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334 $ bnrm1, bnrm2, eps, onepls, safmax, safmin,
335 $ salfai, salfar, sbeta, scale, temp
336* ..
337* .. Local Arrays ..
338 LOGICAL LDUMMA( 1 )
339* ..
340* .. External Subroutines ..
341 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slacpy,
343* ..
344* .. External Functions ..
345 LOGICAL LSAME
346 INTEGER ILAENV
347 REAL SLAMCH, SLANGE
348 EXTERNAL ilaenv, lsame, slamch, slange
349* ..
350* .. Intrinsic Functions ..
351 INTRINSIC abs, int, max
352* ..
353* .. Executable Statements ..
354*
355* Decode the input arguments
356*
357 IF( lsame( jobvl, 'N' ) ) THEN
358 ijobvl = 1
359 ilvl = .false.
360 ELSE IF( lsame( jobvl, 'V' ) ) THEN
361 ijobvl = 2
362 ilvl = .true.
363 ELSE
364 ijobvl = -1
365 ilvl = .false.
366 END IF
367*
368 IF( lsame( jobvr, 'N' ) ) THEN
369 ijobvr = 1
370 ilvr = .false.
371 ELSE IF( lsame( jobvr, 'V' ) ) THEN
372 ijobvr = 2
373 ilvr = .true.
374 ELSE
375 ijobvr = -1
376 ilvr = .false.
377 END IF
378 ilv = ilvl .OR. ilvr
379*
380* Test the input arguments
381*
382 lwkmin = max( 8*n, 1 )
383 lwkopt = lwkmin
384 work( 1 ) = lwkopt
385 lquery = ( lwork.EQ.-1 )
386 info = 0
387 IF( ijobvl.LE.0 ) THEN
388 info = -1
389 ELSE IF( ijobvr.LE.0 ) THEN
390 info = -2
391 ELSE IF( n.LT.0 ) THEN
392 info = -3
393 ELSE IF( lda.LT.max( 1, n ) ) THEN
394 info = -5
395 ELSE IF( ldb.LT.max( 1, n ) ) THEN
396 info = -7
397 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
398 info = -12
399 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
400 info = -14
401 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
402 info = -16
403 END IF
404*
405 IF( info.EQ.0 ) THEN
406 nb1 = ilaenv( 1, 'SGEQRF', ' ', n, n, -1, -1 )
407 nb2 = ilaenv( 1, 'SORMQR', ' ', n, n, n, -1 )
408 nb3 = ilaenv( 1, 'SORGQR', ' ', n, n, n, -1 )
409 nb = max( nb1, nb2, nb3 )
410 lopt = 2*n + max( 6*n, n*(nb+1) )
411 work( 1 ) = lopt
412 END IF
413*
414 IF( info.NE.0 ) THEN
415 CALL xerbla( 'SGEGV ', -info )
416 RETURN
417 ELSE IF( lquery ) THEN
418 RETURN
419 END IF
420*
421* Quick return if possible
422*
423 IF( n.EQ.0 )
424 $ RETURN
425*
426* Get machine constants
427*
428 eps = slamch( 'E' )*slamch( 'B' )
429 safmin = slamch( 'S' )
430 safmin = safmin + safmin
431 safmax = one / safmin
432 onepls = one + ( 4*eps )
433*
434* Scale A
435*
436 anrm = slange( 'M', n, n, a, lda, work )
437 anrm1 = anrm
438 anrm2 = one
439 IF( anrm.LT.one ) THEN
440 IF( safmax*anrm.LT.one ) THEN
441 anrm1 = safmin
442 anrm2 = safmax*anrm
443 END IF
444 END IF
445*
446 IF( anrm.GT.zero ) THEN
447 CALL slascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
448 IF( iinfo.NE.0 ) THEN
449 info = n + 10
450 RETURN
451 END IF
452 END IF
453*
454* Scale B
455*
456 bnrm = slange( 'M', n, n, b, ldb, work )
457 bnrm1 = bnrm
458 bnrm2 = one
459 IF( bnrm.LT.one ) THEN
460 IF( safmax*bnrm.LT.one ) THEN
461 bnrm1 = safmin
462 bnrm2 = safmax*bnrm
463 END IF
464 END IF
465*
466 IF( bnrm.GT.zero ) THEN
467 CALL slascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
468 IF( iinfo.NE.0 ) THEN
469 info = n + 10
470 RETURN
471 END IF
472 END IF
473*
474* Permute the matrix to make it more nearly triangular
475* Workspace layout: (8*N words -- "work" requires 6*N words)
476* left_permutation, right_permutation, work...
477*
478 ileft = 1
479 iright = n + 1
480 iwork = iright + n
481 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
482 $ work( iright ), work( iwork ), iinfo )
483 IF( iinfo.NE.0 ) THEN
484 info = n + 1
485 GO TO 120
486 END IF
487*
488* Reduce B to triangular form, and initialize VL and/or VR
489* Workspace layout: ("work..." must have at least N words)
490* left_permutation, right_permutation, tau, work...
491*
492 irows = ihi + 1 - ilo
493 IF( ilv ) THEN
494 icols = n + 1 - ilo
495 ELSE
496 icols = irows
497 END IF
498 itau = iwork
499 iwork = itau + irows
500 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
501 $ work( iwork ), lwork+1-iwork, iinfo )
502 IF( iinfo.GE.0 )
503 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
504 IF( iinfo.NE.0 ) THEN
505 info = n + 2
506 GO TO 120
507 END IF
508*
509 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
510 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
511 $ lwork+1-iwork, iinfo )
512 IF( iinfo.GE.0 )
513 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514 IF( iinfo.NE.0 ) THEN
515 info = n + 3
516 GO TO 120
517 END IF
518*
519 IF( ilvl ) THEN
520 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
521 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
522 $ vl( ilo+1, ilo ), ldvl )
523 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
524 $ work( itau ), work( iwork ), lwork+1-iwork,
525 $ iinfo )
526 IF( iinfo.GE.0 )
527 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
528 IF( iinfo.NE.0 ) THEN
529 info = n + 4
530 GO TO 120
531 END IF
532 END IF
533*
534 IF( ilvr )
535 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
536*
537* Reduce to generalized Hessenberg form
538*
539 IF( ilv ) THEN
540*
541* Eigenvectors requested -- work on whole matrix.
542*
543 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
544 $ ldvl, vr, ldvr, iinfo )
545 ELSE
546 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
547 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548 END IF
549 IF( iinfo.NE.0 ) THEN
550 info = n + 5
551 GO TO 120
552 END IF
553*
554* Perform QZ algorithm
555* Workspace layout: ("work..." must have at least 1 word)
556* left_permutation, right_permutation, work...
557*
558 iwork = itau
559 IF( ilv ) THEN
560 chtemp = 'S'
561 ELSE
562 chtemp = 'E'
563 END IF
564 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
565 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
566 $ work( iwork ), lwork+1-iwork, iinfo )
567 IF( iinfo.GE.0 )
568 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
569 IF( iinfo.NE.0 ) THEN
570 IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
571 info = iinfo
572 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
573 info = iinfo - n
574 ELSE
575 info = n + 6
576 END IF
577 GO TO 120
578 END IF
579*
580 IF( ilv ) THEN
581*
582* Compute Eigenvectors (STGEVC requires 6*N words of workspace)
583*
584 IF( ilvl ) THEN
585 IF( ilvr ) THEN
586 chtemp = 'B'
587 ELSE
588 chtemp = 'L'
589 END IF
590 ELSE
591 chtemp = 'R'
592 END IF
593*
594 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
595 $ vr, ldvr, n, in, work( iwork ), iinfo )
596 IF( iinfo.NE.0 ) THEN
597 info = n + 7
598 GO TO 120
599 END IF
600*
601* Undo balancing on VL and VR, rescale
602*
603 IF( ilvl ) THEN
604 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
605 $ work( iright ), n, vl, ldvl, iinfo )
606 IF( iinfo.NE.0 ) THEN
607 info = n + 8
608 GO TO 120
609 END IF
610 DO 50 jc = 1, n
611 IF( alphai( jc ).LT.zero )
612 $ GO TO 50
613 temp = zero
614 IF( alphai( jc ).EQ.zero ) THEN
615 DO 10 jr = 1, n
616 temp = max( temp, abs( vl( jr, jc ) ) )
617 10 CONTINUE
618 ELSE
619 DO 20 jr = 1, n
620 temp = max( temp, abs( vl( jr, jc ) )+
621 $ abs( vl( jr, jc+1 ) ) )
622 20 CONTINUE
623 END IF
624 IF( temp.LT.safmin )
625 $ GO TO 50
626 temp = one / temp
627 IF( alphai( jc ).EQ.zero ) THEN
628 DO 30 jr = 1, n
629 vl( jr, jc ) = vl( jr, jc )*temp
630 30 CONTINUE
631 ELSE
632 DO 40 jr = 1, n
633 vl( jr, jc ) = vl( jr, jc )*temp
634 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
635 40 CONTINUE
636 END IF
637 50 CONTINUE
638 END IF
639 IF( ilvr ) THEN
640 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
641 $ work( iright ), n, vr, ldvr, iinfo )
642 IF( iinfo.NE.0 ) THEN
643 info = n + 9
644 GO TO 120
645 END IF
646 DO 100 jc = 1, n
647 IF( alphai( jc ).LT.zero )
648 $ GO TO 100
649 temp = zero
650 IF( alphai( jc ).EQ.zero ) THEN
651 DO 60 jr = 1, n
652 temp = max( temp, abs( vr( jr, jc ) ) )
653 60 CONTINUE
654 ELSE
655 DO 70 jr = 1, n
656 temp = max( temp, abs( vr( jr, jc ) )+
657 $ abs( vr( jr, jc+1 ) ) )
658 70 CONTINUE
659 END IF
660 IF( temp.LT.safmin )
661 $ GO TO 100
662 temp = one / temp
663 IF( alphai( jc ).EQ.zero ) THEN
664 DO 80 jr = 1, n
665 vr( jr, jc ) = vr( jr, jc )*temp
666 80 CONTINUE
667 ELSE
668 DO 90 jr = 1, n
669 vr( jr, jc ) = vr( jr, jc )*temp
670 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
671 90 CONTINUE
672 END IF
673 100 CONTINUE
674 END IF
675*
676* End of eigenvector calculation
677*
678 END IF
679*
680* Undo scaling in alpha, beta
681*
682* Note: this does not give the alpha and beta for the unscaled
683* problem.
684*
685* Un-scaling is limited to avoid underflow in alpha and beta
686* if they are significant.
687*
688 DO 110 jc = 1, n
689 absar = abs( alphar( jc ) )
690 absai = abs( alphai( jc ) )
691 absb = abs( beta( jc ) )
692 salfar = anrm*alphar( jc )
693 salfai = anrm*alphai( jc )
694 sbeta = bnrm*beta( jc )
695 ilimit = .false.
696 scale = one
697*
698* Check for significant underflow in ALPHAI
699*
700 IF( abs( salfai ).LT.safmin .AND. absai.GE.
701 $ max( safmin, eps*absar, eps*absb ) ) THEN
702 ilimit = .true.
703 scale = ( onepls*safmin / anrm1 ) /
704 $ max( onepls*safmin, anrm2*absai )
705*
706 ELSE IF( salfai.EQ.zero ) THEN
707*
708* If insignificant underflow in ALPHAI, then make the
709* conjugate eigenvalue real.
710*
711 IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
712 alphai( jc-1 ) = zero
713 ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
714 alphai( jc+1 ) = zero
715 END IF
716 END IF
717*
718* Check for significant underflow in ALPHAR
719*
720 IF( abs( salfar ).LT.safmin .AND. absar.GE.
721 $ max( safmin, eps*absai, eps*absb ) ) THEN
722 ilimit = .true.
723 scale = max( scale, ( onepls*safmin / anrm1 ) /
724 $ max( onepls*safmin, anrm2*absar ) )
725 END IF
726*
727* Check for significant underflow in BETA
728*
729 IF( abs( sbeta ).LT.safmin .AND. absb.GE.
730 $ max( safmin, eps*absar, eps*absai ) ) THEN
731 ilimit = .true.
732 scale = max( scale, ( onepls*safmin / bnrm1 ) /
733 $ max( onepls*safmin, bnrm2*absb ) )
734 END IF
735*
736* Check for possible overflow when limiting scaling
737*
738 IF( ilimit ) THEN
739 temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
740 $ abs( sbeta ) )
741 IF( temp.GT.one )
742 $ scale = scale / temp
743 IF( scale.LT.one )
744 $ ilimit = .false.
745 END IF
746*
747* Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748*
749 IF( ilimit ) THEN
750 salfar = ( scale*alphar( jc ) )*anrm
751 salfai = ( scale*alphai( jc ) )*anrm
752 sbeta = ( scale*beta( jc ) )*bnrm
753 END IF
754 alphar( jc ) = salfar
755 alphai( jc ) = salfai
756 beta( jc ) = sbeta
757 110 CONTINUE
758*
759 120 CONTINUE
760 work( 1 ) = lwkopt
761*
762 RETURN
763*
764* End of SGEGV
765*
766 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
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 sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
Definition sgghrd.f:207
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 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 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 stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC
Definition stgevc.f:295
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 sgegv(jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
SGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix p...
Definition sgegv.f:306