DTGSNA  - 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 DTGSNA( 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( * )

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


PURPOSE

       DTGSNA 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 DGGES),
       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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DTGEVC.  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) DOUBLE PRECISION 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 DTGEVC.  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) DOUBLE PRECISION 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) DOUBLE PRECISION 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)  DOUBLE  PRECISION  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



FURTHER DETAILS

       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


          Note that if the default method for computing DIF(i) is
       wanted
          (see  DLATDF),  then  the  parameter DIFDRI (see below)
       should be
          changed from 3 to 4 (routine DLATDF(IJOB =  2  will  be
       used)).
          See DTGSYL 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


          and d2 is an upper  bound  on  Difl((S11,  T11),  (S22,
       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
          DLATDF),  then  the parameter DIFDRI (see below) should
       be changed
          from 3 to 4 (routine DLATDF(IJOB = 2  will  be  used)).
       See DTGSYL
          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


       [3] B. Kagstrom and P.  Poromaa,  LAPACK-Style  Algorithms
       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