stgsna(3S)

STGSNA - estimate reciprocal condition numbers for specified 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

As shipped in IRIX 6.5.19. Added in IRIX 6.5.15.

NAME
     STGSNA - estimate reciprocal condition numbers for specified 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( * )

IMPLEMENTATION
     These routines are part of the SCSL Scientific Library and can be loaded
     using either the -lscs or the -lscs_mp option.  The -lscs_mp option
     directs the linker to use the multi-processor version of the library.

     When linking to SCSL with -lscs or -lscs_mp, the default integer size is
     4 bytes (32 bits). Another version of SCSL is available in which integers
     are 8 bytes (64 bits).  This version allows the user access to larger
     memory sizes and helps when porting legacy Cray codes.  It can be loaded
     by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
     only one of the two versions; 4-byte integer and 8-byte integer library
     calls cannot be mixed.

PURPOSE
     STGSNA estimates reciprocal condition numbers for specified 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 eigenpairs;
             = '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 for which
             condition numbers are required. To select condition numbers for
             the eigenpair corresponding to a real eigenvalue w(j), SELECT(j)
             must be set to .TRUE.. To select condition numbers corresponding
             to a complex conjugate pair of eigenvalues 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 eigenvectors 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 referenced.

     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 eigenvectors 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 referenced.

     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 numbers of the
             selected eigenvalues, stored in consecutive elements of the
             array. For a complex conjugate pair of eigenvalues two
             consecutive elements of S are set to the same value. Thus S(j),
             DIF(j), and the j-th columns of VL and VR all correspond 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 condition 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

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 corresponding 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.

        Note that if the default method for computing DIF(i) 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.

     b) If the i-th and (i+1)-th eigenvalues are complex conjugate 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.

        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
        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 eigenvector 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 Eigenvalues 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 Condition
         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.

     [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.


SEE ALSO
     INTRO_LAPACK(3S), INTRO_SCSL(3S)

     This man page is available only online.