LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
stgevc.f
Go to the documentation of this file.
1*> \brief \b STGEVC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download STGEVC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgevc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgevc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgevc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
22* LDVL, VR, LDVR, MM, M, WORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER HOWMNY, SIDE
26* INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
27* ..
28* .. Array Arguments ..
29* LOGICAL SELECT( * )
30* REAL P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
31* $ VR( LDVR, * ), WORK( * )
32* ..
33*
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> STGEVC computes some or all of the right and/or left eigenvectors of
42*> a pair of real matrices (S,P), where S is a quasi-triangular matrix
43*> and P is upper triangular. Matrix pairs of this type are produced by
44*> the generalized Schur factorization of a matrix pair (A,B):
45*>
46*> A = Q*S*Z**T, B = Q*P*Z**T
47*>
48*> as computed by SGGHRD + SHGEQZ.
49*>
50*> The right eigenvector x and the left eigenvector y of (S,P)
51*> corresponding to an eigenvalue w are defined by:
52*>
53*> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
54*>
55*> where y**H denotes the conjugate transpose of y.
56*> The eigenvalues are not input to this routine, but are computed
57*> directly from the diagonal blocks of S and P.
58*>
59*> This routine returns the matrices X and/or Y of right and left
60*> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
61*> where Z and Q are input matrices.
62*> If Q and Z are the orthogonal factors from the generalized Schur
63*> factorization of a matrix pair (A,B), then Z*X and Q*Y
64*> are the matrices of right and left eigenvectors of (A,B).
65*>
66*> \endverbatim
67*
68* Arguments:
69* ==========
70*
71*> \param[in] SIDE
72*> \verbatim
73*> SIDE is CHARACTER*1
74*> = 'R': compute right eigenvectors only;
75*> = 'L': compute left eigenvectors only;
76*> = 'B': compute both right and left eigenvectors.
77*> \endverbatim
78*>
79*> \param[in] HOWMNY
80*> \verbatim
81*> HOWMNY is CHARACTER*1
82*> = 'A': compute all right and/or left eigenvectors;
83*> = 'B': compute all right and/or left eigenvectors,
84*> backtransformed by the matrices in VR and/or VL;
85*> = 'S': compute selected right and/or left eigenvectors,
86*> specified by the logical array SELECT.
87*> \endverbatim
88*>
89*> \param[in] SELECT
90*> \verbatim
91*> SELECT is LOGICAL array, dimension (N)
92*> If HOWMNY='S', SELECT specifies the eigenvectors to be
93*> computed. If w(j) is a real eigenvalue, the corresponding
94*> real eigenvector is computed if SELECT(j) is .TRUE..
95*> If w(j) and w(j+1) are the real and imaginary parts of a
96*> complex eigenvalue, the corresponding complex eigenvector
97*> is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
98*> and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
99*> set to .FALSE..
100*> Not referenced if HOWMNY = 'A' or 'B'.
101*> \endverbatim
102*>
103*> \param[in] N
104*> \verbatim
105*> N is INTEGER
106*> The order of the matrices S and P. N >= 0.
107*> \endverbatim
108*>
109*> \param[in] S
110*> \verbatim
111*> S is REAL array, dimension (LDS,N)
112*> The upper quasi-triangular matrix S from a generalized Schur
113*> factorization, as computed by SHGEQZ.
114*> \endverbatim
115*>
116*> \param[in] LDS
117*> \verbatim
118*> LDS is INTEGER
119*> The leading dimension of array S. LDS >= max(1,N).
120*> \endverbatim
121*>
122*> \param[in] P
123*> \verbatim
124*> P is REAL array, dimension (LDP,N)
125*> The upper triangular matrix P from a generalized Schur
126*> factorization, as computed by SHGEQZ.
127*> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
128*> of S must be in positive diagonal form.
129*> \endverbatim
130*>
131*> \param[in] LDP
132*> \verbatim
133*> LDP is INTEGER
134*> The leading dimension of array P. LDP >= max(1,N).
135*> \endverbatim
136*>
137*> \param[in,out] VL
138*> \verbatim
139*> VL is REAL array, dimension (LDVL,MM)
140*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
141*> contain an N-by-N matrix Q (usually the orthogonal matrix Q
142*> of left Schur vectors returned by SHGEQZ).
143*> On exit, if SIDE = 'L' or 'B', VL contains:
144*> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
145*> if HOWMNY = 'B', the matrix Q*Y;
146*> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
147*> SELECT, stored consecutively in the columns of
148*> VL, in the same order as their eigenvalues.
149*>
150*> A complex eigenvector corresponding to a complex eigenvalue
151*> is stored in two consecutive columns, the first holding the
152*> real part, and the second the imaginary part.
153*>
154*> Not referenced if SIDE = 'R'.
155*> \endverbatim
156*>
157*> \param[in] LDVL
158*> \verbatim
159*> LDVL is INTEGER
160*> The leading dimension of array VL. LDVL >= 1, and if
161*> SIDE = 'L' or 'B', LDVL >= N.
162*> \endverbatim
163*>
164*> \param[in,out] VR
165*> \verbatim
166*> VR is REAL array, dimension (LDVR,MM)
167*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
168*> contain an N-by-N matrix Z (usually the orthogonal matrix Z
169*> of right Schur vectors returned by SHGEQZ).
170*>
171*> On exit, if SIDE = 'R' or 'B', VR contains:
172*> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
173*> if HOWMNY = 'B' or 'b', the matrix Z*X;
174*> if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
175*> specified by SELECT, stored consecutively in the
176*> columns of VR, in the same order as their
177*> eigenvalues.
178*>
179*> A complex eigenvector corresponding to a complex eigenvalue
180*> is stored in two consecutive columns, the first holding the
181*> real part and the second the imaginary part.
182*>
183*> Not referenced if SIDE = 'L'.
184*> \endverbatim
185*>
186*> \param[in] LDVR
187*> \verbatim
188*> LDVR is INTEGER
189*> The leading dimension of the array VR. LDVR >= 1, and if
190*> SIDE = 'R' or 'B', LDVR >= N.
191*> \endverbatim
192*>
193*> \param[in] MM
194*> \verbatim
195*> MM is INTEGER
196*> The number of columns in the arrays VL and/or VR. MM >= M.
197*> \endverbatim
198*>
199*> \param[out] M
200*> \verbatim
201*> M is INTEGER
202*> The number of columns in the arrays VL and/or VR actually
203*> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
204*> is set to N. Each selected real eigenvector occupies one
205*> column and each selected complex eigenvector occupies two
206*> columns.
207*> \endverbatim
208*>
209*> \param[out] WORK
210*> \verbatim
211*> WORK is REAL array, dimension (6*N)
212*> \endverbatim
213*>
214*> \param[out] INFO
215*> \verbatim
216*> INFO is INTEGER
217*> = 0: successful exit.
218*> < 0: if INFO = -i, the i-th argument had an illegal value.
219*> > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex
220*> eigenvalue.
221*> \endverbatim
222*
223* Authors:
224* ========
225*
226*> \author Univ. of Tennessee
227*> \author Univ. of California Berkeley
228*> \author Univ. of Colorado Denver
229*> \author NAG Ltd.
230*
231*> \ingroup tgevc
232*
233*> \par Further Details:
234* =====================
235*>
236*> \verbatim
237*>
238*> Allocation of workspace:
239*> ---------- -- ---------
240*>
241*> WORK( j ) = 1-norm of j-th column of A, above the diagonal
242*> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
243*> WORK( 2*N+1:3*N ) = real part of eigenvector
244*> WORK( 3*N+1:4*N ) = imaginary part of eigenvector
245*> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
246*> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
247*>
248*> Rowwise vs. columnwise solution methods:
249*> ------- -- ---------- -------- -------
250*>
251*> Finding a generalized eigenvector consists basically of solving the
252*> singular triangular system
253*>
254*> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
255*>
256*> Consider finding the i-th right eigenvector (assume all eigenvalues
257*> are real). The equation to be solved is:
258*> n i
259*> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
260*> k=j k=j
261*>
262*> where C = (A - w B) (The components v(i+1:n) are 0.)
263*>
264*> The "rowwise" method is:
265*>
266*> (1) v(i) := 1
267*> for j = i-1,. . .,1:
268*> i
269*> (2) compute s = - sum C(j,k) v(k) and
270*> k=j+1
271*>
272*> (3) v(j) := s / C(j,j)
273*>
274*> Step 2 is sometimes called the "dot product" step, since it is an
275*> inner product between the j-th row and the portion of the eigenvector
276*> that has been computed so far.
277*>
278*> The "columnwise" method consists basically in doing the sums
279*> for all the rows in parallel. As each v(j) is computed, the
280*> contribution of v(j) times the j-th column of C is added to the
281*> partial sums. Since FORTRAN arrays are stored columnwise, this has
282*> the advantage that at each step, the elements of C that are accessed
283*> are adjacent to one another, whereas with the rowwise method, the
284*> elements accessed at a step are spaced LDS (and LDP) words apart.
285*>
286*> When finding left eigenvectors, the matrix in question is the
287*> transpose of the one in storage, so the rowwise method then
288*> actually accesses columns of A and B at each step, and so is the
289*> preferred method.
290*> \endverbatim
291*>
292* =====================================================================
293 SUBROUTINE stgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
294 $ LDVL, VR, LDVR, MM, M, WORK, INFO )
295*
296* -- LAPACK computational routine --
297* -- LAPACK is a software package provided by Univ. of Tennessee, --
298* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299*
300* .. Scalar Arguments ..
301 CHARACTER HOWMNY, SIDE
302 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
303* ..
304* .. Array Arguments ..
305 LOGICAL SELECT( * )
306 REAL P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
307 $ vr( ldvr, * ), work( * )
308* ..
309*
310*
311* =====================================================================
312*
313* .. Parameters ..
314 REAL ZERO, ONE, SAFETY
315 parameter( zero = 0.0e+0, one = 1.0e+0,
316 $ safety = 1.0e+2 )
317* ..
318* .. Local Scalars ..
319 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
320 $ ilbbad, ilcomp, ilcplx, lsa, lsb
321 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
322 $ j, ja, jc, je, jr, jw, na, nw
323 REAL ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
324 $ bcoefr, big, bignum, bnorm, bscale, cim2a,
325 $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
326 $ crealb, dmin, safmin, salfar, sbeta, scale,
327 $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
328 $ xscale
329* ..
330* .. Local Arrays ..
331 REAL BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
332 $ sump( 2, 2 )
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 REAL SLAMCH
337 EXTERNAL lsame, slamch
338* ..
339* .. External Subroutines ..
340 EXTERNAL sgemv, slacpy, slag2, slaln2, xerbla
341* ..
342* .. Intrinsic Functions ..
343 INTRINSIC abs, max, min
344* ..
345* .. Executable Statements ..
346*
347* Decode and Test the input parameters
348*
349 IF( lsame( howmny, 'A' ) ) THEN
350 ihwmny = 1
351 ilall = .true.
352 ilback = .false.
353 ELSE IF( lsame( howmny, 'S' ) ) THEN
354 ihwmny = 2
355 ilall = .false.
356 ilback = .false.
357 ELSE IF( lsame( howmny, 'B' ) ) THEN
358 ihwmny = 3
359 ilall = .true.
360 ilback = .true.
361 ELSE
362 ihwmny = -1
363 ilall = .true.
364 END IF
365*
366 IF( lsame( side, 'R' ) ) THEN
367 iside = 1
368 compl = .false.
369 compr = .true.
370 ELSE IF( lsame( side, 'L' ) ) THEN
371 iside = 2
372 compl = .true.
373 compr = .false.
374 ELSE IF( lsame( side, 'B' ) ) THEN
375 iside = 3
376 compl = .true.
377 compr = .true.
378 ELSE
379 iside = -1
380 END IF
381*
382 info = 0
383 IF( iside.LT.0 ) THEN
384 info = -1
385 ELSE IF( ihwmny.LT.0 ) THEN
386 info = -2
387 ELSE IF( n.LT.0 ) THEN
388 info = -4
389 ELSE IF( lds.LT.max( 1, n ) ) THEN
390 info = -6
391 ELSE IF( ldp.LT.max( 1, n ) ) THEN
392 info = -8
393 END IF
394 IF( info.NE.0 ) THEN
395 CALL xerbla( 'STGEVC', -info )
396 RETURN
397 END IF
398*
399* Count the number of eigenvectors to be computed
400*
401 IF( .NOT.ilall ) THEN
402 im = 0
403 ilcplx = .false.
404 DO 10 j = 1, n
405 IF( ilcplx ) THEN
406 ilcplx = .false.
407 GO TO 10
408 END IF
409 IF( j.LT.n ) THEN
410 IF( s( j+1, j ).NE.zero )
411 $ ilcplx = .true.
412 END IF
413 IF( ilcplx ) THEN
414 IF( SELECT( j ) .OR. SELECT( j+1 ) )
415 $ im = im + 2
416 ELSE
417 IF( SELECT( j ) )
418 $ im = im + 1
419 END IF
420 10 CONTINUE
421 ELSE
422 im = n
423 END IF
424*
425* Check 2-by-2 diagonal blocks of A, B
426*
427 ilabad = .false.
428 ilbbad = .false.
429 DO 20 j = 1, n - 1
430 IF( s( j+1, j ).NE.zero ) THEN
431 IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
432 $ p( j, j+1 ).NE.zero )ilbbad = .true.
433 IF( j.LT.n-1 ) THEN
434 IF( s( j+2, j+1 ).NE.zero )
435 $ ilabad = .true.
436 END IF
437 END IF
438 20 CONTINUE
439*
440 IF( ilabad ) THEN
441 info = -5
442 ELSE IF( ilbbad ) THEN
443 info = -7
444 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
445 info = -10
446 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
447 info = -12
448 ELSE IF( mm.LT.im ) THEN
449 info = -13
450 END IF
451 IF( info.NE.0 ) THEN
452 CALL xerbla( 'STGEVC', -info )
453 RETURN
454 END IF
455*
456* Quick return if possible
457*
458 m = im
459 IF( n.EQ.0 )
460 $ RETURN
461*
462* Machine Constants
463*
464 safmin = slamch( 'Safe minimum' )
465 big = one / safmin
466 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
467 small = safmin*n / ulp
468 big = one / small
469 bignum = one / ( safmin*n )
470*
471* Compute the 1-norm of each column of the strictly upper triangular
472* part (i.e., excluding all elements belonging to the diagonal
473* blocks) of A and B to check for possible overflow in the
474* triangular solver.
475*
476 anorm = abs( s( 1, 1 ) )
477 IF( n.GT.1 )
478 $ anorm = anorm + abs( s( 2, 1 ) )
479 bnorm = abs( p( 1, 1 ) )
480 work( 1 ) = zero
481 work( n+1 ) = zero
482*
483 DO 50 j = 2, n
484 temp = zero
485 temp2 = zero
486 IF( s( j, j-1 ).EQ.zero ) THEN
487 iend = j - 1
488 ELSE
489 iend = j - 2
490 END IF
491 DO 30 i = 1, iend
492 temp = temp + abs( s( i, j ) )
493 temp2 = temp2 + abs( p( i, j ) )
494 30 CONTINUE
495 work( j ) = temp
496 work( n+j ) = temp2
497 DO 40 i = iend + 1, min( j+1, n )
498 temp = temp + abs( s( i, j ) )
499 temp2 = temp2 + abs( p( i, j ) )
500 40 CONTINUE
501 anorm = max( anorm, temp )
502 bnorm = max( bnorm, temp2 )
503 50 CONTINUE
504*
505 ascale = one / max( anorm, safmin )
506 bscale = one / max( bnorm, safmin )
507*
508* Left eigenvectors
509*
510 IF( compl ) THEN
511 ieig = 0
512*
513* Main loop over eigenvalues
514*
515 ilcplx = .false.
516 DO 220 je = 1, n
517*
518* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
519* (b) this would be the second of a complex pair.
520* Check for complex eigenvalue, so as to be sure of which
521* entry(-ies) of SELECT to look at.
522*
523 IF( ilcplx ) THEN
524 ilcplx = .false.
525 GO TO 220
526 END IF
527 nw = 1
528 IF( je.LT.n ) THEN
529 IF( s( je+1, je ).NE.zero ) THEN
530 ilcplx = .true.
531 nw = 2
532 END IF
533 END IF
534 IF( ilall ) THEN
535 ilcomp = .true.
536 ELSE IF( ilcplx ) THEN
537 ilcomp = SELECT( je ) .OR. SELECT( je+1 )
538 ELSE
539 ilcomp = SELECT( je )
540 END IF
541 IF( .NOT.ilcomp )
542 $ GO TO 220
543*
544* Decide if (a) singular pencil, (b) real eigenvalue, or
545* (c) complex eigenvalue.
546*
547 IF( .NOT.ilcplx ) THEN
548 IF( abs( s( je, je ) ).LE.safmin .AND.
549 $ abs( p( je, je ) ).LE.safmin ) THEN
550*
551* Singular matrix pencil -- return unit eigenvector
552*
553 ieig = ieig + 1
554 DO 60 jr = 1, n
555 vl( jr, ieig ) = zero
556 60 CONTINUE
557 vl( ieig, ieig ) = one
558 GO TO 220
559 END IF
560 END IF
561*
562* Clear vector
563*
564 DO 70 jr = 1, nw*n
565 work( 2*n+jr ) = zero
566 70 CONTINUE
567* T
568* Compute coefficients in ( a A - b B ) y = 0
569* a is ACOEF
570* b is BCOEFR + i*BCOEFI
571*
572 IF( .NOT.ilcplx ) THEN
573*
574* Real eigenvalue
575*
576 temp = one / max( abs( s( je, je ) )*ascale,
577 $ abs( p( je, je ) )*bscale, safmin )
578 salfar = ( temp*s( je, je ) )*ascale
579 sbeta = ( temp*p( je, je ) )*bscale
580 acoef = sbeta*ascale
581 bcoefr = salfar*bscale
582 bcoefi = zero
583*
584* Scale to avoid underflow
585*
586 scale = one
587 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
588 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
589 $ small
590 IF( lsa )
591 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
592 IF( lsb )
593 $ scale = max( scale, ( small / abs( salfar ) )*
594 $ min( bnorm, big ) )
595 IF( lsa .OR. lsb ) THEN
596 scale = min( scale, one /
597 $ ( safmin*max( one, abs( acoef ),
598 $ abs( bcoefr ) ) ) )
599 IF( lsa ) THEN
600 acoef = ascale*( scale*sbeta )
601 ELSE
602 acoef = scale*acoef
603 END IF
604 IF( lsb ) THEN
605 bcoefr = bscale*( scale*salfar )
606 ELSE
607 bcoefr = scale*bcoefr
608 END IF
609 END IF
610 acoefa = abs( acoef )
611 bcoefa = abs( bcoefr )
612*
613* First component is 1
614*
615 work( 2*n+je ) = one
616 xmax = one
617 ELSE
618*
619* Complex eigenvalue
620*
621 CALL slag2( s( je, je ), lds, p( je, je ), ldp,
622 $ safmin*safety, acoef, temp, bcoefr, temp2,
623 $ bcoefi )
624 bcoefi = -bcoefi
625 IF( bcoefi.EQ.zero ) THEN
626 info = je
627 RETURN
628 END IF
629*
630* Scale to avoid over/underflow
631*
632 acoefa = abs( acoef )
633 bcoefa = abs( bcoefr ) + abs( bcoefi )
634 scale = one
635 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
636 $ scale = ( safmin / ulp ) / acoefa
637 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
638 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
639 IF( safmin*acoefa.GT.ascale )
640 $ scale = ascale / ( safmin*acoefa )
641 IF( safmin*bcoefa.GT.bscale )
642 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
643 IF( scale.NE.one ) THEN
644 acoef = scale*acoef
645 acoefa = abs( acoef )
646 bcoefr = scale*bcoefr
647 bcoefi = scale*bcoefi
648 bcoefa = abs( bcoefr ) + abs( bcoefi )
649 END IF
650*
651* Compute first two components of eigenvector
652*
653 temp = acoef*s( je+1, je )
654 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
655 temp2i = -bcoefi*p( je, je )
656 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
657 work( 2*n+je ) = one
658 work( 3*n+je ) = zero
659 work( 2*n+je+1 ) = -temp2r / temp
660 work( 3*n+je+1 ) = -temp2i / temp
661 ELSE
662 work( 2*n+je+1 ) = one
663 work( 3*n+je+1 ) = zero
664 temp = acoef*s( je, je+1 )
665 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
666 $ s( je+1, je+1 ) ) / temp
667 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
668 END IF
669 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
670 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
671 END IF
672*
673 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
674*
675* T
676* Triangular solve of (a A - b B) y = 0
677*
678* T
679* (rowwise in (a A - b B) , or columnwise in (a A - b B) )
680*
681 il2by2 = .false.
682*
683 DO 160 j = je + nw, n
684 IF( il2by2 ) THEN
685 il2by2 = .false.
686 GO TO 160
687 END IF
688*
689 na = 1
690 bdiag( 1 ) = p( j, j )
691 IF( j.LT.n ) THEN
692 IF( s( j+1, j ).NE.zero ) THEN
693 il2by2 = .true.
694 bdiag( 2 ) = p( j+1, j+1 )
695 na = 2
696 END IF
697 END IF
698*
699* Check whether scaling is necessary for dot products
700*
701 xscale = one / max( one, xmax )
702 temp = max( work( j ), work( n+j ),
703 $ acoefa*work( j )+bcoefa*work( n+j ) )
704 IF( il2by2 )
705 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
706 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
707 IF( temp.GT.bignum*xscale ) THEN
708 DO 90 jw = 0, nw - 1
709 DO 80 jr = je, j - 1
710 work( ( jw+2 )*n+jr ) = xscale*
711 $ work( ( jw+2 )*n+jr )
712 80 CONTINUE
713 90 CONTINUE
714 xmax = xmax*xscale
715 END IF
716*
717* Compute dot products
718*
719* j-1
720* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
721* k=je
722*
723* To reduce the op count, this is done as
724*
725* _ j-1 _ j-1
726* a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
727* k=je k=je
728*
729* which may cause underflow problems if A or B are close
730* to underflow. (E.g., less than SMALL.)
731*
732*
733 DO 120 jw = 1, nw
734 DO 110 ja = 1, na
735 sums( ja, jw ) = zero
736 sump( ja, jw ) = zero
737*
738 DO 100 jr = je, j - 1
739 sums( ja, jw ) = sums( ja, jw ) +
740 $ s( jr, j+ja-1 )*
741 $ work( ( jw+1 )*n+jr )
742 sump( ja, jw ) = sump( ja, jw ) +
743 $ p( jr, j+ja-1 )*
744 $ work( ( jw+1 )*n+jr )
745 100 CONTINUE
746 110 CONTINUE
747 120 CONTINUE
748*
749 DO 130 ja = 1, na
750 IF( ilcplx ) THEN
751 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
752 $ bcoefr*sump( ja, 1 ) -
753 $ bcoefi*sump( ja, 2 )
754 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
755 $ bcoefr*sump( ja, 2 ) +
756 $ bcoefi*sump( ja, 1 )
757 ELSE
758 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
759 $ bcoefr*sump( ja, 1 )
760 END IF
761 130 CONTINUE
762*
763* T
764* Solve ( a A - b B ) y = SUM(,)
765* with scaling and perturbation of the denominator
766*
767 CALL slaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
768 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
769 $ bcoefi, work( 2*n+j ), n, scale, temp,
770 $ iinfo )
771 IF( scale.LT.one ) THEN
772 DO 150 jw = 0, nw - 1
773 DO 140 jr = je, j - 1
774 work( ( jw+2 )*n+jr ) = scale*
775 $ work( ( jw+2 )*n+jr )
776 140 CONTINUE
777 150 CONTINUE
778 xmax = scale*xmax
779 END IF
780 xmax = max( xmax, temp )
781 160 CONTINUE
782*
783* Copy eigenvector to VL, back transforming if
784* HOWMNY='B'.
785*
786 ieig = ieig + 1
787 IF( ilback ) THEN
788 DO 170 jw = 0, nw - 1
789 CALL sgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
790 $ work( ( jw+2 )*n+je ), 1, zero,
791 $ work( ( jw+4 )*n+1 ), 1 )
792 170 CONTINUE
793 CALL slacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
794 $ ldvl )
795 ibeg = 1
796 ELSE
797 CALL slacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
798 $ ldvl )
799 ibeg = je
800 END IF
801*
802* Scale eigenvector
803*
804 xmax = zero
805 IF( ilcplx ) THEN
806 DO 180 j = ibeg, n
807 xmax = max( xmax, abs( vl( j, ieig ) )+
808 $ abs( vl( j, ieig+1 ) ) )
809 180 CONTINUE
810 ELSE
811 DO 190 j = ibeg, n
812 xmax = max( xmax, abs( vl( j, ieig ) ) )
813 190 CONTINUE
814 END IF
815*
816 IF( xmax.GT.safmin ) THEN
817 xscale = one / xmax
818*
819 DO 210 jw = 0, nw - 1
820 DO 200 jr = ibeg, n
821 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
822 200 CONTINUE
823 210 CONTINUE
824 END IF
825 ieig = ieig + nw - 1
826*
827 220 CONTINUE
828 END IF
829*
830* Right eigenvectors
831*
832 IF( compr ) THEN
833 ieig = im + 1
834*
835* Main loop over eigenvalues
836*
837 ilcplx = .false.
838 DO 500 je = n, 1, -1
839*
840* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
841* (b) this would be the second of a complex pair.
842* Check for complex eigenvalue, so as to be sure of which
843* entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
844* or SELECT(JE-1).
845* If this is a complex pair, the 2-by-2 diagonal block
846* corresponding to the eigenvalue is in rows/columns JE-1:JE
847*
848 IF( ilcplx ) THEN
849 ilcplx = .false.
850 GO TO 500
851 END IF
852 nw = 1
853 IF( je.GT.1 ) THEN
854 IF( s( je, je-1 ).NE.zero ) THEN
855 ilcplx = .true.
856 nw = 2
857 END IF
858 END IF
859 IF( ilall ) THEN
860 ilcomp = .true.
861 ELSE IF( ilcplx ) THEN
862 ilcomp = SELECT( je ) .OR. SELECT( je-1 )
863 ELSE
864 ilcomp = SELECT( je )
865 END IF
866 IF( .NOT.ilcomp )
867 $ GO TO 500
868*
869* Decide if (a) singular pencil, (b) real eigenvalue, or
870* (c) complex eigenvalue.
871*
872 IF( .NOT.ilcplx ) THEN
873 IF( abs( s( je, je ) ).LE.safmin .AND.
874 $ abs( p( je, je ) ).LE.safmin ) THEN
875*
876* Singular matrix pencil -- unit eigenvector
877*
878 ieig = ieig - 1
879 DO 230 jr = 1, n
880 vr( jr, ieig ) = zero
881 230 CONTINUE
882 vr( ieig, ieig ) = one
883 GO TO 500
884 END IF
885 END IF
886*
887* Clear vector
888*
889 DO 250 jw = 0, nw - 1
890 DO 240 jr = 1, n
891 work( ( jw+2 )*n+jr ) = zero
892 240 CONTINUE
893 250 CONTINUE
894*
895* Compute coefficients in ( a A - b B ) x = 0
896* a is ACOEF
897* b is BCOEFR + i*BCOEFI
898*
899 IF( .NOT.ilcplx ) THEN
900*
901* Real eigenvalue
902*
903 temp = one / max( abs( s( je, je ) )*ascale,
904 $ abs( p( je, je ) )*bscale, safmin )
905 salfar = ( temp*s( je, je ) )*ascale
906 sbeta = ( temp*p( je, je ) )*bscale
907 acoef = sbeta*ascale
908 bcoefr = salfar*bscale
909 bcoefi = zero
910*
911* Scale to avoid underflow
912*
913 scale = one
914 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
915 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
916 $ small
917 IF( lsa )
918 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
919 IF( lsb )
920 $ scale = max( scale, ( small / abs( salfar ) )*
921 $ min( bnorm, big ) )
922 IF( lsa .OR. lsb ) THEN
923 scale = min( scale, one /
924 $ ( safmin*max( one, abs( acoef ),
925 $ abs( bcoefr ) ) ) )
926 IF( lsa ) THEN
927 acoef = ascale*( scale*sbeta )
928 ELSE
929 acoef = scale*acoef
930 END IF
931 IF( lsb ) THEN
932 bcoefr = bscale*( scale*salfar )
933 ELSE
934 bcoefr = scale*bcoefr
935 END IF
936 END IF
937 acoefa = abs( acoef )
938 bcoefa = abs( bcoefr )
939*
940* First component is 1
941*
942 work( 2*n+je ) = one
943 xmax = one
944*
945* Compute contribution from column JE of A and B to sum
946* (See "Further Details", above.)
947*
948 DO 260 jr = 1, je - 1
949 work( 2*n+jr ) = bcoefr*p( jr, je ) -
950 $ acoef*s( jr, je )
951 260 CONTINUE
952 ELSE
953*
954* Complex eigenvalue
955*
956 CALL slag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
957 $ safmin*safety, acoef, temp, bcoefr, temp2,
958 $ bcoefi )
959 IF( bcoefi.EQ.zero ) THEN
960 info = je - 1
961 RETURN
962 END IF
963*
964* Scale to avoid over/underflow
965*
966 acoefa = abs( acoef )
967 bcoefa = abs( bcoefr ) + abs( bcoefi )
968 scale = one
969 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
970 $ scale = ( safmin / ulp ) / acoefa
971 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
972 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
973 IF( safmin*acoefa.GT.ascale )
974 $ scale = ascale / ( safmin*acoefa )
975 IF( safmin*bcoefa.GT.bscale )
976 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
977 IF( scale.NE.one ) THEN
978 acoef = scale*acoef
979 acoefa = abs( acoef )
980 bcoefr = scale*bcoefr
981 bcoefi = scale*bcoefi
982 bcoefa = abs( bcoefr ) + abs( bcoefi )
983 END IF
984*
985* Compute first two components of eigenvector
986* and contribution to sums
987*
988 temp = acoef*s( je, je-1 )
989 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
990 temp2i = -bcoefi*p( je, je )
991 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
992 work( 2*n+je ) = one
993 work( 3*n+je ) = zero
994 work( 2*n+je-1 ) = -temp2r / temp
995 work( 3*n+je-1 ) = -temp2i / temp
996 ELSE
997 work( 2*n+je-1 ) = one
998 work( 3*n+je-1 ) = zero
999 temp = acoef*s( je-1, je )
1000 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1001 $ s( je-1, je-1 ) ) / temp
1002 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1003 END IF
1004*
1005 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1006 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1007*
1008* Compute contribution from columns JE and JE-1
1009* of A and B to the sums.
1010*
1011 creala = acoef*work( 2*n+je-1 )
1012 cimaga = acoef*work( 3*n+je-1 )
1013 crealb = bcoefr*work( 2*n+je-1 ) -
1014 $ bcoefi*work( 3*n+je-1 )
1015 cimagb = bcoefi*work( 2*n+je-1 ) +
1016 $ bcoefr*work( 3*n+je-1 )
1017 cre2a = acoef*work( 2*n+je )
1018 cim2a = acoef*work( 3*n+je )
1019 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1020 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1021 DO 270 jr = 1, je - 2
1022 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1023 $ crealb*p( jr, je-1 ) -
1024 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1025 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1026 $ cimagb*p( jr, je-1 ) -
1027 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1028 270 CONTINUE
1029 END IF
1030*
1031 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1032*
1033* Columnwise triangular solve of (a A - b B) x = 0
1034*
1035 il2by2 = .false.
1036 DO 370 j = je - nw, 1, -1
1037*
1038* If a 2-by-2 block, is in position j-1:j, wait until
1039* next iteration to process it (when it will be j:j+1)
1040*
1041 IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1042 IF( s( j, j-1 ).NE.zero ) THEN
1043 il2by2 = .true.
1044 GO TO 370
1045 END IF
1046 END IF
1047 bdiag( 1 ) = p( j, j )
1048 IF( il2by2 ) THEN
1049 na = 2
1050 bdiag( 2 ) = p( j+1, j+1 )
1051 ELSE
1052 na = 1
1053 END IF
1054*
1055* Compute x(j) (and x(j+1), if 2-by-2 block)
1056*
1057 CALL slaln2( .false., na, nw, dmin, acoef, s( j, j ),
1058 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1059 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1060 $ iinfo )
1061 IF( scale.LT.one ) THEN
1062*
1063 DO 290 jw = 0, nw - 1
1064 DO 280 jr = 1, je
1065 work( ( jw+2 )*n+jr ) = scale*
1066 $ work( ( jw+2 )*n+jr )
1067 280 CONTINUE
1068 290 CONTINUE
1069 END IF
1070 xmax = max( scale*xmax, temp )
1071*
1072 DO 310 jw = 1, nw
1073 DO 300 ja = 1, na
1074 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1075 300 CONTINUE
1076 310 CONTINUE
1077*
1078* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1079*
1080 IF( j.GT.1 ) THEN
1081*
1082* Check whether scaling is necessary for sum.
1083*
1084 xscale = one / max( one, xmax )
1085 temp = acoefa*work( j ) + bcoefa*work( n+j )
1086 IF( il2by2 )
1087 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1088 $ work( n+j+1 ) )
1089 temp = max( temp, acoefa, bcoefa )
1090 IF( temp.GT.bignum*xscale ) THEN
1091*
1092 DO 330 jw = 0, nw - 1
1093 DO 320 jr = 1, je
1094 work( ( jw+2 )*n+jr ) = xscale*
1095 $ work( ( jw+2 )*n+jr )
1096 320 CONTINUE
1097 330 CONTINUE
1098 xmax = xmax*xscale
1099 END IF
1100*
1101* Compute the contributions of the off-diagonals of
1102* column j (and j+1, if 2-by-2 block) of A and B to the
1103* sums.
1104*
1105*
1106 DO 360 ja = 1, na
1107 IF( ilcplx ) THEN
1108 creala = acoef*work( 2*n+j+ja-1 )
1109 cimaga = acoef*work( 3*n+j+ja-1 )
1110 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1111 $ bcoefi*work( 3*n+j+ja-1 )
1112 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1113 $ bcoefr*work( 3*n+j+ja-1 )
1114 DO 340 jr = 1, j - 1
1115 work( 2*n+jr ) = work( 2*n+jr ) -
1116 $ creala*s( jr, j+ja-1 ) +
1117 $ crealb*p( jr, j+ja-1 )
1118 work( 3*n+jr ) = work( 3*n+jr ) -
1119 $ cimaga*s( jr, j+ja-1 ) +
1120 $ cimagb*p( jr, j+ja-1 )
1121 340 CONTINUE
1122 ELSE
1123 creala = acoef*work( 2*n+j+ja-1 )
1124 crealb = bcoefr*work( 2*n+j+ja-1 )
1125 DO 350 jr = 1, j - 1
1126 work( 2*n+jr ) = work( 2*n+jr ) -
1127 $ creala*s( jr, j+ja-1 ) +
1128 $ crealb*p( jr, j+ja-1 )
1129 350 CONTINUE
1130 END IF
1131 360 CONTINUE
1132 END IF
1133*
1134 il2by2 = .false.
1135 370 CONTINUE
1136*
1137* Copy eigenvector to VR, back transforming if
1138* HOWMNY='B'.
1139*
1140 ieig = ieig - nw
1141 IF( ilback ) THEN
1142*
1143 DO 410 jw = 0, nw - 1
1144 DO 380 jr = 1, n
1145 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1146 $ vr( jr, 1 )
1147 380 CONTINUE
1148*
1149* A series of compiler directives to defeat
1150* vectorization for the next loop
1151*
1152*
1153 DO 400 jc = 2, je
1154 DO 390 jr = 1, n
1155 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1156 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1157 390 CONTINUE
1158 400 CONTINUE
1159 410 CONTINUE
1160*
1161 DO 430 jw = 0, nw - 1
1162 DO 420 jr = 1, n
1163 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1164 420 CONTINUE
1165 430 CONTINUE
1166*
1167 iend = n
1168 ELSE
1169 DO 450 jw = 0, nw - 1
1170 DO 440 jr = 1, n
1171 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1172 440 CONTINUE
1173 450 CONTINUE
1174*
1175 iend = je
1176 END IF
1177*
1178* Scale eigenvector
1179*
1180 xmax = zero
1181 IF( ilcplx ) THEN
1182 DO 460 j = 1, iend
1183 xmax = max( xmax, abs( vr( j, ieig ) )+
1184 $ abs( vr( j, ieig+1 ) ) )
1185 460 CONTINUE
1186 ELSE
1187 DO 470 j = 1, iend
1188 xmax = max( xmax, abs( vr( j, ieig ) ) )
1189 470 CONTINUE
1190 END IF
1191*
1192 IF( xmax.GT.safmin ) THEN
1193 xscale = one / xmax
1194 DO 490 jw = 0, nw - 1
1195 DO 480 jr = 1, iend
1196 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1197 480 CONTINUE
1198 490 CONTINUE
1199 END IF
1200 500 CONTINUE
1201 END IF
1202*
1203 RETURN
1204*
1205* End of STGEVC
1206*
1207 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
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 slag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition slag2.f:156
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition slaln2.f:218
subroutine stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC
Definition stgevc.f:295