stgevc.f

famous linear algebra library (LAPACK) ports to windows

      SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
     $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          HOWMNY, SIDE
      INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      REAL               P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
     $                   VR( LDVR, * ), WORK( * )
*     ..
*
*
*  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).
* 
*  Arguments
*  =========
*
*  SIDE    (input) CHARACTER*1
*          = 'R': compute right eigenvectors only;
*          = 'L': compute left eigenvectors only;
*          = 'B': compute both right and left eigenvectors.
*
*  HOWMNY  (input) 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.
*
*  SELECT  (input) 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'.
*
*  N       (input) INTEGER
*          The order of the matrices S and P.  N >= 0.
*
*  S       (input) REAL array, dimension (LDS,N)
*          The upper quasi-triangular matrix S from a generalized Schur
*          factorization, as computed by SHGEQZ.
*
*  LDS     (input) INTEGER
*          The leading dimension of array S.  LDS >= max(1,N).
*
*  P       (input) 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.
*
*  LDP     (input) INTEGER
*          The leading dimension of array P.  LDP >= max(1,N).
*
*  VL      (input/output) 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'.
*
*  LDVL    (input) INTEGER
*          The leading dimension of array VL.  LDVL >= 1, and if
*          SIDE = 'L' or 'B', LDVL >= N.
*
*  VR      (input/output) 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'.
*
*  LDVR    (input) INTEGER
*          The leading dimension of the array VR.  LDVR >= 1, and if
*          SIDE = 'R' or 'B', LDVR >= N.
*
*  MM      (input) INTEGER
*          The number of columns in the arrays VL and/or VR. MM >= M.
*
*  M       (output) 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.
*
*  WORK    (workspace) REAL array, dimension (6*N)
*
*  INFO    (output) 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.
*
*  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.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, SAFETY
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0,
     $                   SAFETY = 1.0E+2 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
     $                   ILBBAD, ILCOMP, ILCPLX, LSA, LSB
      INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
     $                   J, JA, JC, JE, JR, JW, NA, NW
      REAL               ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
     $                   BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
     $                   CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
     $                   CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
     $                   SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
     $                   XSCALE
*     ..
*     .. Local Arrays ..
      REAL               BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
     $                   SUMP( 2, 2 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH
      EXTERNAL           LSAME, SLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMV, SLABAD, SLACPY, SLAG2, SLALN2, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Decode and Test the input parameters
*
      IF( LSAME( HOWMNY, 'A' ) ) THEN
         IHWMNY = 1
         ILALL = .TRUE.
         ILBACK = .FALSE.
      ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
         IHWMNY = 2
         ILALL = .FALSE.
         ILBACK = .FALSE.
      ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
         IHWMNY = 3
         ILALL = .TRUE.
         ILBACK = .TRUE.
      ELSE
         IHWMNY = -1
         ILALL = .TRUE.
      END IF
*
      IF( LSAME( SIDE, 'R' ) ) THEN
         ISIDE = 1
         COMPL = .FALSE.
         COMPR = .TRUE.
      ELSE IF( LSAME( SIDE, 'L' ) ) THEN
         ISIDE = 2
         COMPL = .TRUE.
         COMPR = .FALSE.
      ELSE IF( LSAME( SIDE, 'B' ) ) THEN
         ISIDE = 3
         COMPL = .TRUE.
         COMPR = .TRUE.
      ELSE
         ISIDE = -1
      END IF
*
      INFO = 0
      IF( ISIDE.LT.0 ) THEN
         INFO = -1
      ELSE IF( IHWMNY.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'STGEVC', -INFO )
         RETURN
      END IF
*
*     Count the number of eigenvectors to be computed
*
      IF( .NOT.ILALL ) THEN
         IM = 0
         ILCPLX = .FALSE.
         DO 10 J = 1, N
            IF( ILCPLX ) THEN
               ILCPLX = .FALSE.
               GO TO 10
            END IF
            IF( J.LT.N ) THEN
               IF( S( J+1, J ).NE.ZERO )
     $            ILCPLX = .TRUE.
            END IF
            IF( ILCPLX ) THEN
               IF( SELECT( J ) .OR. SELECT( J+1 ) )
     $            IM = IM + 2
            ELSE
               IF( SELECT( J ) )
     $            IM = IM + 1
            END IF
   10    CONTINUE
      ELSE
         IM = N
      END IF
*
*     Check 2-by-2 diagonal blocks of A, B
*
      ILABAD = .FALSE.
      ILBBAD = .FALSE.
      DO 20 J = 1, N - 1
         IF( S( J+1, J ).NE.ZERO ) THEN
            IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
     $          P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
            IF( J.LT.N-1 ) THEN
               IF( S( J+2, J+1 ).NE.ZERO )
     $            ILABAD = .TRUE.
            END IF
         END IF
   20 CONTINUE
*
      IF( ILABAD ) THEN
         INFO = -5
      ELSE IF( ILBBAD ) THEN
         INFO = -7
      ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
         INFO = -10
      ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
         INFO = -12
      ELSE IF( MM.LT.IM ) THEN
         INFO = -13
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'STGEVC', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      M = IM
      IF( N.EQ.0 )
     $   RETURN
*
*     Machine Constants
*
      SAFMIN = SLAMCH( 'Safe minimum' )
      BIG = ONE / SAFMIN
      CALL SLABAD( SAFMIN, BIG )
      ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
      SMALL = SAFMIN*N / ULP
      BIG = ONE / SMALL
      BIGNUM = ONE / ( SAFMIN*N )
*
*     Compute the 1-norm of each column of the strictly upper triangular
*     part (i.e., excluding all elements belonging to the diagonal
*     blocks) of A and B to check for possible overflow in the
*     triangular solver.
*
      ANORM = ABS( S( 1, 1 ) )
      IF( N.GT.1 )
     $   ANORM = ANORM + ABS( S( 2, 1 ) )
      BNORM = ABS( P( 1, 1 ) )
      WORK( 1 ) = ZERO
      WORK( N+1 ) = ZERO
*
      DO 50 J = 2, N