STGSNA  - estimate reciprocal condition numbers for speci­
       fied eigenvalues and/or eigenvectors of a matrix pair  (A,
       B)  in  generalized  real  Schur canonical form (or of any
       matrix pair (Q*A*Z', Q*B*Z') with  orthogonal  matrices  Q
       and Z, where Z' denotes the transpose of Z


SYNOPSIS

       SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB,
                          VL, LDVL, VR,  LDVR,  S,  DIF,  MM,  M,
                          WORK, LWORK, IWORK, INFO )

           CHARACTER      HOWMNY, JOB

           INTEGER        INFO,  LDA,  LDB, LDVL, LDVR, LWORK, M,
                          MM, N

           LOGICAL        SELECT( * )

           INTEGER        IWORK( * )

           REAL           A( LDA, * ), B( LDB, * ), DIF( * ),  S(
                          *  ),  VL(  LDVL,  *  ), VR( LDVR, * ),
                          WORK( * )


PURPOSE

       STGSNA estimates reciprocal condition numbers  for  speci­
       fied  eigenvalues and/or eigenvectors of a matrix pair (A,
       B) in generalized real Schur canonical  form  (or  of  any
       matrix  pair  (Q*A*Z',  Q*B*Z') with orthogonal matrices Q
       and Z, where Z' denotes the transpose of Z.  (A,  B)  must
       be  in generalized real Schur form (as returned by SGGES),
       i.e. A is block upper triangular with  1-by-1  and  2-by-2
       diagonal blocks. B is upper triangular.


ARGUMENTS

       JOB     (input) CHARACTER*1
               Specifies  whether  condition numbers are required
               for eigenvalues (S) or eigenvectors (DIF):
               = 'E': for eigenvalues only (S);
               = 'V': for eigenvectors only (DIF);
               = 'B': for both eigenvalues  and  eigenvectors  (S
               and DIF).

       HOWMNY  (input) CHARACTER*1
               =  'A':  compute  condition numbers for all eigen­
               pairs;
               = 'S':  compute  condition  numbers  for  selected
               eigenpairs specified by the array SELECT.

       SELECT  (input) LOGICAL array, dimension (N)
               If  HOWMNY  = 'S', SELECT specifies the eigenpairs

               select  condition numbers for the eigenpair corre­
               sponding to a real eigenvalue w(j), SELECT(j) must
               be set to .TRUE.. To select condition numbers cor­
               responding to a complex conjugate pair  of  eigen­
               values   w(j)  and  w(j+1),  either  SELECT(j)  or
               SELECT(j+1) or both, must be set  to  .TRUE..   If
               HOWMNY = 'A', SELECT is not referenced.

       N       (input) INTEGER
               The  order  of the square matrix pair (A, B). N >=
               0.

       A       (input) REAL array, dimension (LDA,N)
               The upper quasi-triangular matrix A  in  the  pair
               (A,B).

       LDA     (input) INTEGER
               The  leading  dimension  of  the  array  A. LDA >=
               max(1,N).

       B       (input) REAL array, dimension (LDB,N)
               The upper triangular matrix B in the pair (A,B).

       LDB     (input) INTEGER
               The leading dimension  of  the  array  B.  LDB  >=
               max(1,N).

       VL      (input) REAL array, dimension (LDVL,M)
               If  JOB  = 'E' or 'B', VL must contain left eigen­
               vectors of (A, B), corresponding to the eigenpairs
               specified  by  HOWMNY and SELECT. The eigenvectors
               must be stored in consecutive columns  of  VL,  as
               returned  by STGEVC.  If JOB = 'V', VL is not ref­
               erenced.

       LDVL    (input) INTEGER
               The leading dimension of the array VL. LDVL >=  1.
               If JOB = 'E' or 'B', LDVL >= N.

       VR      (input) REAL array, dimension (LDVR,M)
               If  JOB = 'E' or 'B', VR must contain right eigen­
               vectors of (A, B), corresponding to the eigenpairs
               specified  by  HOWMNY and SELECT. The eigenvectors
               must be stored in consecutive columns  ov  VR,  as
               returned  by STGEVC.  If JOB = 'V', VR is not ref­
               erenced.

       LDVR    (input) INTEGER
               The leading dimension of the array VR. LDVR >=  1.
               If JOB = 'E' or 'B', LDVR >= N.

       S       (output) REAL array, dimension (MM)
               If  JOB  =  'E'  or  'B', the reciprocal condition

               consecutive  elements  of the array. For a complex
               conjugate pair of eigenvalues two consecutive ele­
               ments  of  S are set to the same value. Thus S(j),
               DIF(j), and the j-th columns of VL and VR all cor­
               respond  to the same eigenpair (but not in general
               the j-th  eigenpair,  unless  all  eigenpairs  are
               selected).  If JOB = 'V', S is not referenced.

       DIF     (output) REAL array, dimension (MM)
               If JOB = 'V' or 'B', the estimated reciprocal con­
               dition  numbers  of  the  selected   eigenvectors,
               stored in consecutive elements of the array. For a
               complex eigenvector two  consecutive  elements  of
               DIF  are set to the same value. If the eigenvalues
               cannot be reordered to compute DIF(j),  DIF(j)  is
               set  to 0; this can only occur when the true value
               would be very small anyway.  If JOB = 'E', DIF  is
               not referenced.

       MM      (input) INTEGER
               The number of elements in the arrays S and DIF. MM
               >= M.

       M       (output) INTEGER
               The number of elements of the  arrays  S  and  DIF
               used to store the specified condition numbers; for
               each selected real eigenvalue one element is used,
               and  for  each  selected complex conjugate pair of
               eigenvalues, two elements are used.  If  HOWMNY  =
               'A', M is set to N.

       WORK    (workspace/output) REAL array, dimension (LWORK)
               If  JOB = 'E', WORK is not referenced.  Otherwise,
               on exit, if INFO = 0, WORK(1) returns the  optimal
               LWORK.

       LWORK   (input) INTEGER
               The  dimension  of the array WORK. LWORK >= N.  If
               JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.

               If LWORK = -1, then a workspace query is  assumed;
               the  routine  only  calculates the optimal size of
               the WORK array, returns this value  as  the  first
               entry  of  the  WORK  array,  and no error message
               related to LWORK is issued by XERBLA.

       IWORK   (workspace) INTEGER array, dimension (N + 6)
               If JOB = 'E', IWORK is not referenced.

       INFO    (output) INTEGER
               =0: Successful exit
               <0: If INFO = -i, the i-th argument had an illegal
               value

       The  reciprocal  of  the condition number of a generalized
       eigenvalue w = (a, b) is defined as

            S(w)   =    (|u'Av|**2    +    |u'Bv|**2)**(1/2)    /
       (norm(u)*norm(v))

       where  u  and v are the left and right eigenvectors of (A,
       B) corresponding to w; |z| denotes the absolute  value  of
       the  complex number, and norm(u) denotes the 2-norm of the
       vector u.
       The pair (a, b) corresponds to an eigenvalue w  =  a/b  (=
       u'Av/u'Bv)  of  the  matrix  pair  (A, B). If both a and b
       equal zero, then (A B)  is  singular  and  S(I)  =  -1  is
       returned.

       An approximate error bound on the chordal distance between
       the i-th computed generalized eigenvalue w and the  corre­
       sponding exact eigenvalue lambda is

            chord(w, lambda) <= EPS * norm(A, B) / S(I)

       where EPS is the machine precision.

       The  reciprocal  of  the  condition number DIF(i) of right
       eigenvector u and left eigenvector v corresponding to  the
       generalized eigenvalue w is defined as follows:

       a) If the i-th eigenvalue w = (a,b) is real

          Suppose  U  and  V  are orthogonal transformations such
       that

                     U'*(A, B)*V  = (S, T) = ( a   *  ) ( b  *  )
       1
                                             ( 0  S22 ),( 0 T22 )
       n-1
                                               1  n-1     1 n-1

          Then the reciprocal condition number DIF(i) is

                     Difl((a, b), (S22, T22)) = sigma-min( Zl ),

          where sigma-min(Zl) denotes the smallest singular value
       of the
          2(n-1)-by-2(n-1) matrix

              Zl = [ kron(a, In-1)  -kron(1, S22) ]
                   [ kron(b, In-1)  -kron(1, T22) ] .

          Here  In-1  is the identity matrix of size n-1. kron(X,
       Y) is the
          Kronecker product between the matrices X and Y.

       wanted
          (see  SLATDF),  then  the  parameter DIFDRI (see below)
       should be
          changed from 3 to 4 (routine SLATDF(IJOB =  2  will  be
       used)).
          See STGSYL for more details.

       b) If the i-th and (i+1)-th eigenvalues are complex conju­
       gate pair,

          Suppose U and V  are  orthogonal  transformations  such
       that

                     U'*(A, B)*V = (S, T) = ( S11  *   ) ( T11  *
       )  2
                                            (  0     S22  ),(   0
       T22) n-2
                                              2       n-2       2
       n-2

          and (S11, T11) corresponds  to  the  complex  conjugate
       eigenvalue
          pair (w, conjg(w)). There exist unitary matrices U1 and
       V1 such
          that

              U1'*S11*V1 = ( s11 s12 )   and U1'*T11*V1 =  (  t11
       t12 )
                           (   0   s22  )                    (  0
       t22 )

          where the generalized eigenvalues w = s11/t11 and
          conjg(w) = s22/t22.

          Then the reciprocal condition number DIF(i) is  bounded
       by

              min( d1, max( 1, |real(s11)/real(s22)| )*d2 )

          where,  d1  =  Difl((s11,  t11),  (s22,  t22)) = sigma-
       min(Z1), where
          Z1 is the complex 2-by-2 matrix

                   Z1 =  [ s11  -s22 ]
                         [ t11  -t22 ],

          This is done by computing (using real arithmetic) the
          roots of the characteristical polynomial det(Z1' * Z1 -
       lambda I),
          where  Z1'  denotes  the  conjugate transpose of Z1 and
       det(X) denotes
          the determinant of X.

       T22)), i.e. an
          upper    bound    on   sigma-min(Z2),   where   Z2   is
       (2n-2)-by-(2n-2)

                   Z2 = [ kron(S11', In-2)  -kron(I2, S22) ]
                        [ kron(T11', In-2)  -kron(I2, T22) ]

          Note that if the default method for  computing  DIF  is
       wanted (see
          SLATDF),  then  the parameter DIFDRI (see below) should
       be changed
          from 3 to 4 (routine SLATDF(IJOB = 2  will  be  used)).
       See STGSYL
          for more details.

       For each eigenvalue/vector specified by SELECT, DIF stores
       a Frobenius norm-based estimate of Difl.

       An approximate error bound for the i-th computed eigenvec­
       tor VL(i) or VR(i) is given by

                  EPS * norm(A, B) / DIF(i).

       See ref. [2-3] for more details and further references.

       Based on contributions by
          Bo  Kagstrom and Peter Poromaa, Department of Computing
       Science,
          Umea University, S-901 87 Umea, Sweden.

       References
       ==========

       [1] B. Kagstrom; A Direct Method for Reordering  Eigenval­
       ues in the
           Generalized  Real  Schur Form of a Regular Matrix Pair
       (A, B), in
           M.S. Moonen et al  (eds),  Linear  Algebra  for  Large
       Scale and
           Real-Time Applications, Kluwer Academic Publ. 1993, pp
       195-218.

       [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with
       Specified
           Eigenvalues of a Regular Matrix Pair (A, B) and Condi­
       tion
           Estimation: Theory, Algorithms and Software,
           Report UMINF - 94.04, Department of Computing Science,
       Umea
           University,  S-901  87  Umea,  Sweden,  1994.  Also as
       LAPACK Working
           Note 87. To appear in Numerical Algorithms, 1996.

       and Software
           for  Solving  the  Generalized  Sylvester Equation and
       Estimating the
           Separation between Regular Matrix Pairs, Report  UMINF
       - 93.23,
           Department  of  Computing  Science,  Umea  University,
       S-901 87 Umea,
           Sweden, December 1993, Revised  April  1994,  Also  as
       LAPACK Working
           Note  75.   To appear in ACM Trans. on Math. Software,
       Vol 22,
           No 1, 1996.


Man(1) output converted with man2html