LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
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: