LAPACK 3.12.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 transpose 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, 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
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
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
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: