LAPACK 3.11.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 tranpose 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 realGEcomputational
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, slabad, 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 CALL slabad( safmin, big )
467 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
468 small = safmin*n / ulp
469 big = one / small
470 bignum = one / ( safmin*n )
471*
472* Compute the 1-norm of each column of the strictly upper triangular
473* part (i.e., excluding all elements belonging to the diagonal
474* blocks) of A and B to check for possible overflow in the
475* triangular solver.
476*
477 anorm = abs( s( 1, 1 ) )
478 IF( n.GT.1 )
479 $ anorm = anorm + abs( s( 2, 1 ) )
480 bnorm = abs( p( 1, 1 ) )
481 work( 1 ) = zero
482 work( n+1 ) = zero
483*
484 DO 50 j = 2, n
485 temp = zero
486 temp2 = zero
487 IF( s( j, j-1 ).EQ.zero ) THEN
488 iend = j - 1
489 ELSE
490 iend = j - 2
491 END IF
492 DO 30 i = 1, iend
493 temp = temp + abs( s( i, j ) )
494 temp2 = temp2 + abs( p( i, j ) )
495 30 CONTINUE
496 work( j ) = temp
497 work( n+j ) = temp2
498 DO 40 i = iend + 1, min( j+1, n )
499 temp = temp + abs( s( i, j ) )
500 temp2 = temp2 + abs( p( i, j ) )
501 40 CONTINUE
502 anorm = max( anorm, temp )
503 bnorm = max( bnorm, temp2 )
504 50 CONTINUE
505*
506 ascale = one / max( anorm, safmin )
507 bscale = one / max( bnorm, safmin )
508*
509* Left eigenvectors
510*
511 IF( compl ) THEN
512 ieig = 0
513*
514* Main loop over eigenvalues
515*
516 ilcplx = .false.
517 DO 220 je = 1, n
518*
519* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
520* (b) this would be the second of a complex pair.
521* Check for complex eigenvalue, so as to be sure of which
522* entry(-ies) of SELECT to look at.
523*
524 IF( ilcplx ) THEN
525 ilcplx = .false.
526 GO TO 220
527 END IF
528 nw = 1
529 IF( je.LT.n ) THEN
530 IF( s( je+1, je ).NE.zero ) THEN
531 ilcplx = .true.
532 nw = 2
533 END IF
534 END IF
535 IF( ilall ) THEN
536 ilcomp = .true.
537 ELSE IF( ilcplx ) THEN
538 ilcomp = SELECT( je ) .OR. SELECT( je+1 )
539 ELSE
540 ilcomp = SELECT( je )
541 END IF
542 IF( .NOT.ilcomp )
543 $ GO TO 220
544*
545* Decide if (a) singular pencil, (b) real eigenvalue, or
546* (c) complex eigenvalue.
547*
548 IF( .NOT.ilcplx ) THEN
549 IF( abs( s( je, je ) ).LE.safmin .AND.
550 $ abs( p( je, je ) ).LE.safmin ) THEN
551*
552* Singular matrix pencil -- return unit eigenvector
553*
554 ieig = ieig + 1
555 DO 60 jr = 1, n
556 vl( jr, ieig ) = zero
557 60 CONTINUE
558 vl( ieig, ieig ) = one
559 GO TO 220
560 END IF
561 END IF
562*
563* Clear vector
564*
565 DO 70 jr = 1, nw*n
566 work( 2*n+jr ) = zero
567 70 CONTINUE
568* T
569* Compute coefficients in ( a A - b B ) y = 0
570* a is ACOEF
571* b is BCOEFR + i*BCOEFI
572*
573 IF( .NOT.ilcplx ) THEN
574*
575* Real eigenvalue
576*
577 temp = one / max( abs( s( je, je ) )*ascale,
578 $ abs( p( je, je ) )*bscale, safmin )
579 salfar = ( temp*s( je, je ) )*ascale
580 sbeta = ( temp*p( je, je ) )*bscale
581 acoef = sbeta*ascale
582 bcoefr = salfar*bscale
583 bcoefi = zero
584*
585* Scale to avoid underflow
586*
587 scale = one
588 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
589 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
590 $ small
591 IF( lsa )
592 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
593 IF( lsb )
594 $ scale = max( scale, ( small / abs( salfar ) )*
595 $ min( bnorm, big ) )
596 IF( lsa .OR. lsb ) THEN
597 scale = min( scale, one /
598 $ ( safmin*max( one, abs( acoef ),
599 $ abs( bcoefr ) ) ) )
600 IF( lsa ) THEN
601 acoef = ascale*( scale*sbeta )
602 ELSE
603 acoef = scale*acoef
604 END IF
605 IF( lsb ) THEN
606 bcoefr = bscale*( scale*salfar )
607 ELSE
608 bcoefr = scale*bcoefr
609 END IF
610 END IF
611 acoefa = abs( acoef )
612 bcoefa = abs( bcoefr )
613*
614* First component is 1
615*
616 work( 2*n+je ) = one
617 xmax = one
618 ELSE
619*
620* Complex eigenvalue
621*
622 CALL slag2( s( je, je ), lds, p( je, je ), ldp,
623 $ safmin*safety, acoef, temp, bcoefr, temp2,
624 $ bcoefi )
625 bcoefi = -bcoefi
626 IF( bcoefi.EQ.zero ) THEN
627 info = je
628 RETURN
629 END IF
630*
631* Scale to avoid over/underflow
632*
633 acoefa = abs( acoef )
634 bcoefa = abs( bcoefr ) + abs( bcoefi )
635 scale = one
636 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
637 $ scale = ( safmin / ulp ) / acoefa
638 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
639 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
640 IF( safmin*acoefa.GT.ascale )
641 $ scale = ascale / ( safmin*acoefa )
642 IF( safmin*bcoefa.GT.bscale )
643 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
644 IF( scale.NE.one ) THEN
645 acoef = scale*acoef
646 acoefa = abs( acoef )
647 bcoefr = scale*bcoefr
648 bcoefi = scale*bcoefi
649 bcoefa = abs( bcoefr ) + abs( bcoefi )
650 END IF
651*
652* Compute first two components of eigenvector
653*
654 temp = acoef*s( je+1, je )
655 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
656 temp2i = -bcoefi*p( je, je )
657 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
658 work( 2*n+je ) = one
659 work( 3*n+je ) = zero
660 work( 2*n+je+1 ) = -temp2r / temp
661 work( 3*n+je+1 ) = -temp2i / temp
662 ELSE
663 work( 2*n+je+1 ) = one
664 work( 3*n+je+1 ) = zero
665 temp = acoef*s( je, je+1 )
666 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
667 $ s( je+1, je+1 ) ) / temp
668 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
669 END IF
670 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
671 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
672 END IF
673*
674 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
675*
676* T
677* Triangular solve of (a A - b B) y = 0
678*
679* T
680* (rowwise in (a A - b B) , or columnwise in (a A - b B) )
681*
682 il2by2 = .false.
683*
684 DO 160 j = je + nw, n
685 IF( il2by2 ) THEN
686 il2by2 = .false.
687 GO TO 160
688 END IF
689*
690 na = 1
691 bdiag( 1 ) = p( j, j )
692 IF( j.LT.n ) THEN
693 IF( s( j+1, j ).NE.zero ) THEN
694 il2by2 = .true.
695 bdiag( 2 ) = p( j+1, j+1 )
696 na = 2
697 END IF
698 END IF
699*
700* Check whether scaling is necessary for dot products
701*
702 xscale = one / max( one, xmax )
703 temp = max( work( j ), work( n+j ),
704 $ acoefa*work( j )+bcoefa*work( n+j ) )
705 IF( il2by2 )
706 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
707 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
708 IF( temp.GT.bignum*xscale ) THEN
709 DO 90 jw = 0, nw - 1
710 DO 80 jr = je, j - 1
711 work( ( jw+2 )*n+jr ) = xscale*
712 $ work( ( jw+2 )*n+jr )
713 80 CONTINUE
714 90 CONTINUE
715 xmax = xmax*xscale
716 END IF
717*
718* Compute dot products
719*
720* j-1
721* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
722* k=je
723*
724* To reduce the op count, this is done as
725*
726* _ j-1 _ j-1
727* a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
728* k=je k=je
729*
730* which may cause underflow problems if A or B are close
731* to underflow. (E.g., less than SMALL.)
732*
733*
734 DO 120 jw = 1, nw
735 DO 110 ja = 1, na
736 sums( ja, jw ) = zero
737 sump( ja, jw ) = zero
738*
739 DO 100 jr = je, j - 1
740 sums( ja, jw ) = sums( ja, jw ) +
741 $ s( jr, j+ja-1 )*
742 $ work( ( jw+1 )*n+jr )
743 sump( ja, jw ) = sump( ja, jw ) +
744 $ p( jr, j+ja-1 )*
745 $ work( ( jw+1 )*n+jr )
746 100 CONTINUE
747 110 CONTINUE
748 120 CONTINUE
749*
750 DO 130 ja = 1, na
751 IF( ilcplx ) THEN
752 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
753 $ bcoefr*sump( ja, 1 ) -
754 $ bcoefi*sump( ja, 2 )
755 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
756 $ bcoefr*sump( ja, 2 ) +
757 $ bcoefi*sump( ja, 1 )
758 ELSE
759 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
760 $ bcoefr*sump( ja, 1 )
761 END IF
762 130 CONTINUE
763*
764* T
765* Solve ( a A - b B ) y = SUM(,)
766* with scaling and perturbation of the denominator
767*
768 CALL slaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
769 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
770 $ bcoefi, work( 2*n+j ), n, scale, temp,
771 $ iinfo )
772 IF( scale.LT.one ) THEN
773 DO 150 jw = 0, nw - 1
774 DO 140 jr = je, j - 1
775 work( ( jw+2 )*n+jr ) = scale*
776 $ work( ( jw+2 )*n+jr )
777 140 CONTINUE
778 150 CONTINUE
779 xmax = scale*xmax
780 END IF
781 xmax = max( xmax, temp )
782 160 CONTINUE
783*
784* Copy eigenvector to VL, back transforming if
785* HOWMNY='B'.
786*
787 ieig = ieig + 1
788 IF( ilback ) THEN
789 DO 170 jw = 0, nw - 1
790 CALL sgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
791 $ work( ( jw+2 )*n+je ), 1, zero,
792 $ work( ( jw+4 )*n+1 ), 1 )
793 170 CONTINUE
794 CALL slacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
795 $ ldvl )
796 ibeg = 1
797 ELSE
798 CALL slacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
799 $ ldvl )
800 ibeg = je
801 END IF
802*
803* Scale eigenvector
804*
805 xmax = zero
806 IF( ilcplx ) THEN
807 DO 180 j = ibeg, n
808 xmax = max( xmax, abs( vl( j, ieig ) )+
809 $ abs( vl( j, ieig+1 ) ) )
810 180 CONTINUE
811 ELSE
812 DO 190 j = ibeg, n
813 xmax = max( xmax, abs( vl( j, ieig ) ) )
814 190 CONTINUE
815 END IF
816*
817 IF( xmax.GT.safmin ) THEN
818 xscale = one / xmax
819*
820 DO 210 jw = 0, nw - 1
821 DO 200 jr = ibeg, n
822 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
823 200 CONTINUE
824 210 CONTINUE
825 END IF
826 ieig = ieig + nw - 1
827*
828 220 CONTINUE
829 END IF
830*
831* Right eigenvectors
832*
833 IF( compr ) THEN
834 ieig = im + 1
835*
836* Main loop over eigenvalues
837*
838 ilcplx = .false.
839 DO 500 je = n, 1, -1
840*
841* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
842* (b) this would be the second of a complex pair.
843* Check for complex eigenvalue, so as to be sure of which
844* entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
845* or SELECT(JE-1).
846* If this is a complex pair, the 2-by-2 diagonal block
847* corresponding to the eigenvalue is in rows/columns JE-1:JE
848*
849 IF( ilcplx ) THEN
850 ilcplx = .false.
851 GO TO 500
852 END IF
853 nw = 1
854 IF( je.GT.1 ) THEN
855 IF( s( je, je-1 ).NE.zero ) THEN
856 ilcplx = .true.
857 nw = 2
858 END IF
859 END IF
860 IF( ilall ) THEN
861 ilcomp = .true.
862 ELSE IF( ilcplx ) THEN
863 ilcomp = SELECT( je ) .OR. SELECT( je-1 )
864 ELSE
865 ilcomp = SELECT( je )
866 END IF
867 IF( .NOT.ilcomp )
868 $ GO TO 500
869*
870* Decide if (a) singular pencil, (b) real eigenvalue, or
871* (c) complex eigenvalue.
872*
873 IF( .NOT.ilcplx ) THEN
874 IF( abs( s( je, je ) ).LE.safmin .AND.
875 $ abs( p( je, je ) ).LE.safmin ) THEN
876*
877* Singular matrix pencil -- unit eigenvector
878*
879 ieig = ieig - 1
880 DO 230 jr = 1, n
881 vr( jr, ieig ) = zero
882 230 CONTINUE
883 vr( ieig, ieig ) = one
884 GO TO 500
885 END IF
886 END IF
887*
888* Clear vector
889*
890 DO 250 jw = 0, nw - 1
891 DO 240 jr = 1, n
892 work( ( jw+2 )*n+jr ) = zero
893 240 CONTINUE
894 250 CONTINUE
895*
896* Compute coefficients in ( a A - b B ) x = 0
897* a is ACOEF
898* b is BCOEFR + i*BCOEFI
899*
900 IF( .NOT.ilcplx ) THEN
901*
902* Real eigenvalue
903*
904 temp = one / max( abs( s( je, je ) )*ascale,
905 $ abs( p( je, je ) )*bscale, safmin )
906 salfar = ( temp*s( je, je ) )*ascale
907 sbeta = ( temp*p( je, je ) )*bscale
908 acoef = sbeta*ascale
909 bcoefr = salfar*bscale
910 bcoefi = zero
911*
912* Scale to avoid underflow
913*
914 scale = one
915 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
916 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
917 $ small
918 IF( lsa )
919 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
920 IF( lsb )
921 $ scale = max( scale, ( small / abs( salfar ) )*
922 $ min( bnorm, big ) )
923 IF( lsa .OR. lsb ) THEN
924 scale = min( scale, one /
925 $ ( safmin*max( one, abs( acoef ),
926 $ abs( bcoefr ) ) ) )
927 IF( lsa ) THEN
928 acoef = ascale*( scale*sbeta )
929 ELSE
930 acoef = scale*acoef
931 END IF
932 IF( lsb ) THEN
933 bcoefr = bscale*( scale*salfar )
934 ELSE
935 bcoefr = scale*bcoefr
936 END IF
937 END IF
938 acoefa = abs( acoef )
939 bcoefa = abs( bcoefr )
940*
941* First component is 1
942*
943 work( 2*n+je ) = one
944 xmax = one
945*
946* Compute contribution from column JE of A and B to sum
947* (See "Further Details", above.)
948*
949 DO 260 jr = 1, je - 1
950 work( 2*n+jr ) = bcoefr*p( jr, je ) -
951 $ acoef*s( jr, je )
952 260 CONTINUE
953 ELSE
954*
955* Complex eigenvalue
956*
957 CALL slag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
958 $ safmin*safety, acoef, temp, bcoefr, temp2,
959 $ bcoefi )
960 IF( bcoefi.EQ.zero ) THEN
961 info = je - 1
962 RETURN
963 END IF
964*
965* Scale to avoid over/underflow
966*
967 acoefa = abs( acoef )
968 bcoefa = abs( bcoefr ) + abs( bcoefi )
969 scale = one
970 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
971 $ scale = ( safmin / ulp ) / acoefa
972 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
973 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
974 IF( safmin*acoefa.GT.ascale )
975 $ scale = ascale / ( safmin*acoefa )
976 IF( safmin*bcoefa.GT.bscale )
977 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
978 IF( scale.NE.one ) THEN
979 acoef = scale*acoef
980 acoefa = abs( acoef )
981 bcoefr = scale*bcoefr
982 bcoefi = scale*bcoefi
983 bcoefa = abs( bcoefr ) + abs( bcoefi )
984 END IF
985*
986* Compute first two components of eigenvector
987* and contribution to sums
988*
989 temp = acoef*s( je, je-1 )
990 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
991 temp2i = -bcoefi*p( je, je )
992 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
993 work( 2*n+je ) = one
994 work( 3*n+je ) = zero
995 work( 2*n+je-1 ) = -temp2r / temp
996 work( 3*n+je-1 ) = -temp2i / temp
997 ELSE
998 work( 2*n+je-1 ) = one
999 work( 3*n+je-1 ) = zero
1000 temp = acoef*s( je-1, je )
1001 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1002 $ s( je-1, je-1 ) ) / temp
1003 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1004 END IF
1005*
1006 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1007 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1008*
1009* Compute contribution from columns JE and JE-1
1010* of A and B to the sums.
1011*
1012 creala = acoef*work( 2*n+je-1 )
1013 cimaga = acoef*work( 3*n+je-1 )
1014 crealb = bcoefr*work( 2*n+je-1 ) -
1015 $ bcoefi*work( 3*n+je-1 )
1016 cimagb = bcoefi*work( 2*n+je-1 ) +
1017 $ bcoefr*work( 3*n+je-1 )
1018 cre2a = acoef*work( 2*n+je )
1019 cim2a = acoef*work( 3*n+je )
1020 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1021 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1022 DO 270 jr = 1, je - 2
1023 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1024 $ crealb*p( jr, je-1 ) -
1025 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1026 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1027 $ cimagb*p( jr, je-1 ) -
1028 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1029 270 CONTINUE
1030 END IF
1031*
1032 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1033*
1034* Columnwise triangular solve of (a A - b B) x = 0
1035*
1036 il2by2 = .false.
1037 DO 370 j = je - nw, 1, -1
1038*
1039* If a 2-by-2 block, is in position j-1:j, wait until
1040* next iteration to process it (when it will be j:j+1)
1041*
1042 IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1043 IF( s( j, j-1 ).NE.zero ) THEN
1044 il2by2 = .true.
1045 GO TO 370
1046 END IF
1047 END IF
1048 bdiag( 1 ) = p( j, j )
1049 IF( il2by2 ) THEN
1050 na = 2
1051 bdiag( 2 ) = p( j+1, j+1 )
1052 ELSE
1053 na = 1
1054 END IF
1055*
1056* Compute x(j) (and x(j+1), if 2-by-2 block)
1057*
1058 CALL slaln2( .false., na, nw, dmin, acoef, s( j, j ),
1059 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1060 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1061 $ iinfo )
1062 IF( scale.LT.one ) THEN
1063*
1064 DO 290 jw = 0, nw - 1
1065 DO 280 jr = 1, je
1066 work( ( jw+2 )*n+jr ) = scale*
1067 $ work( ( jw+2 )*n+jr )
1068 280 CONTINUE
1069 290 CONTINUE
1070 END IF
1071 xmax = max( scale*xmax, temp )
1072*
1073 DO 310 jw = 1, nw
1074 DO 300 ja = 1, na
1075 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1076 300 CONTINUE
1077 310 CONTINUE
1078*
1079* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1080*
1081 IF( j.GT.1 ) THEN
1082*
1083* Check whether scaling is necessary for sum.
1084*
1085 xscale = one / max( one, xmax )
1086 temp = acoefa*work( j ) + bcoefa*work( n+j )
1087 IF( il2by2 )
1088 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1089 $ work( n+j+1 ) )
1090 temp = max( temp, acoefa, bcoefa )
1091 IF( temp.GT.bignum*xscale ) THEN
1092*
1093 DO 330 jw = 0, nw - 1
1094 DO 320 jr = 1, je
1095 work( ( jw+2 )*n+jr ) = xscale*
1096 $ work( ( jw+2 )*n+jr )
1097 320 CONTINUE
1098 330 CONTINUE
1099 xmax = xmax*xscale
1100 END IF
1101*
1102* Compute the contributions of the off-diagonals of
1103* column j (and j+1, if 2-by-2 block) of A and B to the
1104* sums.
1105*
1106*
1107 DO 360 ja = 1, na
1108 IF( ilcplx ) THEN
1109 creala = acoef*work( 2*n+j+ja-1 )
1110 cimaga = acoef*work( 3*n+j+ja-1 )
1111 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1112 $ bcoefi*work( 3*n+j+ja-1 )
1113 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1114 $ bcoefr*work( 3*n+j+ja-1 )
1115 DO 340 jr = 1, j - 1
1116 work( 2*n+jr ) = work( 2*n+jr ) -
1117 $ creala*s( jr, j+ja-1 ) +
1118 $ crealb*p( jr, j+ja-1 )
1119 work( 3*n+jr ) = work( 3*n+jr ) -
1120 $ cimaga*s( jr, j+ja-1 ) +
1121 $ cimagb*p( jr, j+ja-1 )
1122 340 CONTINUE
1123 ELSE
1124 creala = acoef*work( 2*n+j+ja-1 )
1125 crealb = bcoefr*work( 2*n+j+ja-1 )
1126 DO 350 jr = 1, j - 1
1127 work( 2*n+jr ) = work( 2*n+jr ) -
1128 $ creala*s( jr, j+ja-1 ) +
1129 $ crealb*p( jr, j+ja-1 )
1130 350 CONTINUE
1131 END IF
1132 360 CONTINUE
1133 END IF
1134*
1135 il2by2 = .false.
1136 370 CONTINUE
1137*
1138* Copy eigenvector to VR, back transforming if
1139* HOWMNY='B'.
1140*
1141 ieig = ieig - nw
1142 IF( ilback ) THEN
1143*
1144 DO 410 jw = 0, nw - 1
1145 DO 380 jr = 1, n
1146 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1147 $ vr( jr, 1 )
1148 380 CONTINUE
1149*
1150* A series of compiler directives to defeat
1151* vectorization for the next loop
1152*
1153*
1154 DO 400 jc = 2, je
1155 DO 390 jr = 1, n
1156 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1157 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1158 390 CONTINUE
1159 400 CONTINUE
1160 410 CONTINUE
1161*
1162 DO 430 jw = 0, nw - 1
1163 DO 420 jr = 1, n
1164 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1165 420 CONTINUE
1166 430 CONTINUE
1167*
1168 iend = n
1169 ELSE
1170 DO 450 jw = 0, nw - 1
1171 DO 440 jr = 1, n
1172 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1173 440 CONTINUE
1174 450 CONTINUE
1175*
1176 iend = je
1177 END IF
1178*
1179* Scale eigenvector
1180*
1181 xmax = zero
1182 IF( ilcplx ) THEN
1183 DO 460 j = 1, iend
1184 xmax = max( xmax, abs( vr( j, ieig ) )+
1185 $ abs( vr( j, ieig+1 ) ) )
1186 460 CONTINUE
1187 ELSE
1188 DO 470 j = 1, iend
1189 xmax = max( xmax, abs( vr( j, ieig ) ) )
1190 470 CONTINUE
1191 END IF
1192*
1193 IF( xmax.GT.safmin ) THEN
1194 xscale = one / xmax
1195 DO 490 jw = 0, nw - 1
1196 DO 480 jr = 1, iend
1197 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1198 480 CONTINUE
1199 490 CONTINUE
1200 END IF
1201 500 CONTINUE
1202 END IF
1203*
1204 RETURN
1205*
1206* End of STGEVC
1207*
1208 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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 stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:295
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 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 sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156