LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ stgevc()

subroutine stgevc ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
real, dimension( lds, * )  S,
integer  LDS,
real, dimension( ldp, * )  P,
integer  LDP,
real, dimension( ldvl, * )  VL,
integer  LDVL,
real, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
real, dimension( * )  WORK,
integer  INFO 
)

STGEVC

Download STGEVC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STGEVC computes some or all of the right and/or left eigenvectors of
 a pair of real matrices (S,P), where S is a quasi-triangular matrix
 and P is upper triangular.  Matrix pairs of this type are produced by
 the generalized Schur factorization of a matrix pair (A,B):

    A = Q*S*Z**T,  B = Q*P*Z**T

 as computed by SGGHRD + SHGEQZ.

 The right eigenvector x and the left eigenvector y of (S,P)
 corresponding to an eigenvalue w are defined by:

    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,

 where y**H denotes the conjugate tranpose of y.
 The eigenvalues are not input to this routine, but are computed
 directly from the diagonal blocks of S and P.

 This routine returns the matrices X and/or Y of right and left
 eigenvectors of (S,P), or the products Z*X and/or Q*Y,
 where Z and Q are input matrices.
 If Q and Z are the orthogonal factors from the generalized Schur
 factorization of a matrix pair (A,B), then Z*X and Q*Y
 are the matrices of right and left eigenvectors of (A,B).
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'R': compute right eigenvectors only;
          = 'L': compute left eigenvectors only;
          = 'B': compute both right and left eigenvectors.
[in]HOWMNY
          HOWMNY is CHARACTER*1
          = 'A': compute all right and/or left eigenvectors;
          = 'B': compute all right and/or left eigenvectors,
                 backtransformed by the matrices in VR and/or VL;
          = 'S': compute selected right and/or left eigenvectors,
                 specified by the logical array SELECT.
[in]SELECT
          SELECT is LOGICAL array, dimension (N)
          If HOWMNY='S', SELECT specifies the eigenvectors to be
          computed.  If w(j) is a real eigenvalue, the corresponding
          real eigenvector is computed if SELECT(j) is .TRUE..
          If w(j) and w(j+1) are the real and imaginary parts of a
          complex eigenvalue, the corresponding complex eigenvector
          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
          set to .FALSE..
          Not referenced if HOWMNY = 'A' or 'B'.
[in]N
          N is INTEGER
          The order of the matrices S and P.  N >= 0.
[in]S
          S is REAL array, dimension (LDS,N)
          The upper quasi-triangular matrix S from a generalized Schur
          factorization, as computed by SHGEQZ.
[in]LDS
          LDS is INTEGER
          The leading dimension of array S.  LDS >= max(1,N).
[in]P
          P is REAL array, dimension (LDP,N)
          The upper triangular matrix P from a generalized Schur
          factorization, as computed by SHGEQZ.
          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
          of S must be in positive diagonal form.
[in]LDP
          LDP is INTEGER
          The leading dimension of array P.  LDP >= max(1,N).
[in,out]VL
          VL is REAL array, dimension (LDVL,MM)
          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
          contain an N-by-N matrix Q (usually the orthogonal matrix Q
          of left Schur vectors returned by SHGEQZ).
          On exit, if SIDE = 'L' or 'B', VL contains:
          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
          if HOWMNY = 'B', the matrix Q*Y;
          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
                      SELECT, stored consecutively in the columns of
                      VL, in the same order as their eigenvalues.

          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part, and the second the imaginary part.

          Not referenced if SIDE = 'R'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of array VL.  LDVL >= 1, and if
          SIDE = 'L' or 'B', LDVL >= N.
[in,out]VR
          VR is REAL array, dimension (LDVR,MM)
          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
          contain an N-by-N matrix Z (usually the orthogonal matrix Z
          of right Schur vectors returned by SHGEQZ).

          On exit, if SIDE = 'R' or 'B', VR contains:
          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
          if HOWMNY = 'B' or 'b', the matrix Z*X;
          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
                      specified by SELECT, stored consecutively in the
                      columns of VR, in the same order as their
                      eigenvalues.

          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part and the second the imaginary part.

          Not referenced if SIDE = 'L'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1, and if
          SIDE = 'R' or 'B', LDVR >= N.
[in]MM
          MM is INTEGER
          The number of columns in the arrays VL and/or VR. MM >= M.
[out]M
          M is INTEGER
          The number of columns in the arrays VL and/or VR actually
          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
          is set to N.  Each selected real eigenvector occupies one
          column and each selected complex eigenvector occupies two
          columns.
[out]WORK
          WORK is REAL array, dimension (6*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
                eigenvalue.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  Allocation of workspace:
  ---------- -- ---------

     WORK( j ) = 1-norm of j-th column of A, above the diagonal
     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
     WORK( 2*N+1:3*N ) = real part of eigenvector
     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector

  Rowwise vs. columnwise solution methods:
  ------- --  ---------- -------- -------

  Finding a generalized eigenvector consists basically of solving the
  singular triangular system

   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)

  Consider finding the i-th right eigenvector (assume all eigenvalues
  are real). The equation to be solved is:
       n                   i
  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
      k=j                 k=j

  where  C = (A - w B)  (The components v(i+1:n) are 0.)

  The "rowwise" method is:

  (1)  v(i) := 1
  for j = i-1,. . .,1:
                          i
      (2) compute  s = - sum C(j,k) v(k)   and
                        k=j+1

      (3) v(j) := s / C(j,j)

  Step 2 is sometimes called the "dot product" step, since it is an
  inner product between the j-th row and the portion of the eigenvector
  that has been computed so far.

  The "columnwise" method consists basically in doing the sums
  for all the rows in parallel.  As each v(j) is computed, the
  contribution of v(j) times the j-th column of C is added to the
  partial sums.  Since FORTRAN arrays are stored columnwise, this has
  the advantage that at each step, the elements of C that are accessed
  are adjacent to one another, whereas with the rowwise method, the
  elements accessed at a step are spaced LDS (and LDP) words apart.

  When finding left eigenvectors, the matrix in question is the
  transpose of the one in storage, so the rowwise method then
  actually accesses columns of A and B at each step, and so is the
  preferred method.

Definition at line 293 of file stgevc.f.

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*
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
slacpy
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
slaln2
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
slag2
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
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: