ctgsna - estimate reciprocal condition numbers for specified eigenvalues and/or eigenvectors of a matrix pair (A, B)
SUBROUTINE CTGSNA( JOB, HOWMNT, SELECT, N, A, LDA, B, LDB, VL, LDVL, * VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO) CHARACTER * 1 JOB, HOWMNT COMPLEX A(LDA,*), B(LDB,*), VL(LDVL,*), VR(LDVR,*), WORK(*) INTEGER N, LDA, LDB, LDVL, LDVR, MM, M, LWORK, INFO INTEGER IWORK(*) LOGICAL SELECT(*) REAL S(*), DIF(*)
SUBROUTINE CTGSNA_64( JOB, HOWMNT, SELECT, N, A, LDA, B, LDB, VL, * LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO) CHARACTER * 1 JOB, HOWMNT COMPLEX A(LDA,*), B(LDB,*), VL(LDVL,*), VR(LDVR,*), WORK(*) INTEGER*8 N, LDA, LDB, LDVL, LDVR, MM, M, LWORK, INFO INTEGER*8 IWORK(*) LOGICAL*8 SELECT(*) REAL S(*), DIF(*)
SUBROUTINE TGSNA( JOB, HOWMNT, SELECT, [N], A, [LDA], B, [LDB], VL, * [LDVL], VR, [LDVR], S, DIF, MM, M, [WORK], [LWORK], [IWORK], * [INFO]) CHARACTER(LEN=1) :: JOB, HOWMNT COMPLEX, DIMENSION(:) :: WORK COMPLEX, DIMENSION(:,:) :: A, B, VL, VR INTEGER :: N, LDA, LDB, LDVL, LDVR, MM, M, LWORK, INFO INTEGER, DIMENSION(:) :: IWORK LOGICAL, DIMENSION(:) :: SELECT REAL, DIMENSION(:) :: S, DIF
SUBROUTINE TGSNA_64( JOB, HOWMNT, SELECT, [N], A, [LDA], B, [LDB], * VL, [LDVL], VR, [LDVR], S, DIF, MM, M, [WORK], [LWORK], [IWORK], * [INFO]) CHARACTER(LEN=1) :: JOB, HOWMNT COMPLEX, DIMENSION(:) :: WORK COMPLEX, DIMENSION(:,:) :: A, B, VL, VR INTEGER(8) :: N, LDA, LDB, LDVL, LDVR, MM, M, LWORK, INFO INTEGER(8), DIMENSION(:) :: IWORK LOGICAL(8), DIMENSION(:) :: SELECT REAL, DIMENSION(:) :: S, DIF
#include <sunperf.h>
void ctgsna(char job, char howmnt, logical *select, int n, complex *a, int lda, complex *b, int ldb, complex *vl, int ldvl, complex *vr, int ldvr, float *s, float *dif, int mm, int *m, int *info);
void ctgsna_64(char job, char howmnt, logical *select, long n, complex *a, long lda, complex *b, long ldb, complex *vl, long ldvl, complex *vr, long ldvr, float *s, float *dif, long mm, long *m, long *info);
ctgsna estimates reciprocal condition numbers for specified eigenvalues and/or eigenvectors of a matrix pair (A, B).
(A, B) must be in generalized Schur canonical form, that is, A and B are both upper triangular.
= 'E': for eigenvalues only (S);
= 'V': for eigenvectors only (DIF);
= 'B': for both eigenvalues and eigenvectors (S and DIF).
= 'A': compute condition numbers for all eigenpairs;
= 'S': compute condition numbers for selected eigenpairs specified by the array SELECT.
SELECT(j) must be set to .TRUE..
If HOWMNT = 'A', SELECT is not referenced.
DIF(j) is set to 0; this can only occur when the true value
would be very small anyway.
For each eigenvalue/vector specified by SELECT, DIF stores
a Frobenius norm-based estimate of Difl.
If JOB = 'E', DIF is not referenced.
WORK(1) returns the optimal LWORK.
dimension(N+2)
If JOB = 'E', IWORK is not referenced.
= 0: Successful exit
< 0: If INFO = -i, the i-th argument had an illegal value
The reciprocal of the condition number of the i-th generalized eigenvalue w = (a, b) is defined as
S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
where u and v are the right and left 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 ( = v'Au/v'Bu) 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 of the right eigenvector u and left eigenvector v corresponding to the generalized eigenvalue w is defined as follows. Suppose
(A, B) = ( a * ) ( b * ) 1
( 0 A22 ),( 0 B22 ) n-1
1 n-1 1 n-1
Then the reciprocal condition number DIF(I) is
Difl[(a, b), (A22, B22)] = sigma-min( Zl )
where sigma-min(Zl) denotes the smallest singular value of
Zl = [ kron(a, In-1) -kron(1, A22) ]
[ kron(b, In-1) -kron(1, B22) ].
Here In-1 is the identity matrix of size n-1 and X' is the conjugate transpose of X. kron(X, Y) is the Kronecker product between the matrices X and Y.
We approximate the smallest singular value of Zl with an upper bound. This is done by CLATDF.
An approximate error bound for a 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.