Contents
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(*)
F95 INTERFACE
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
C INTERFACE
#include <sunperf.h>
void ctgsna(char job, char howmnt, int *select, int n, com-
plex *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, long *select, long n,
complex *a, long lda, complex *b, long ldb, com-
plex *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.
JOB (input)
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).
HOWMNT (input)
= 'A': compute condition numbers for all eigen-
pairs;
= 'S': compute condition numbers for selected
eigenpairs specified by the array SELECT.
SELECT (input)
If HOWMNT = 'S', SELECT specifies the eigenpairs
for which condition numbers are required. To
select condition numbers for the corresponding j-
th eigenvalue and/or eigenvector, SELECT(j) must
be set to .TRUE.. If HOWMNT = 'A', SELECT is not
referenced.
N (input) The order of the square matrix pair (A, B). N >=
0.
A (input) The upper triangular matrix A in the pair (A,B).
LDA (input)
The leading dimension of the array A. LDA >=
max(1,N).
B (input) The upper triangular matrix B in the pair (A, B).
LDB (input)
The leading dimension of the array B. LDB >=
max(1,N).
VL (input)
If JOB = 'E' or 'B', VL must contain left eigen-
vectors of (A, B), corresponding to the eigenpairs
specified by HOWMNT and SELECT. The eigenvectors
must be stored in consecutive columns of VL, as
returned by CTGEVC. If JOB = 'V', VL is not
referenced.
LDVL (input)
The leading dimension of the array VL. LDVL >= 1;
and If JOB = 'E' or 'B', LDVL >= N.
VR (input)
If JOB = 'E' or 'B', VR must contain right eigen-
vectors of (A, B), corresponding to the eigenpairs
specified by HOWMNT and SELECT. The eigenvectors
must be stored in consecutive columns of VR, as
returned by CTGEVC. If JOB = 'V', VR is not
referenced.
LDVR (input)
The leading dimension of the array VR. LDVR >= 1;
If JOB = 'E' or 'B', LDVR >= N.
S (output)
If JOB = 'E' or 'B', the reciprocal condition
numbers of the selected eigenvalues, stored in
consecutive elements of the array. If JOB = 'V',
S is not referenced.
DIF (output)
If JOB = 'V' or 'B', the estimated reciprocal con-
dition numbers of the selected eigenvectors,
stored in consecutive elements of the array. 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.
For each eigenvalue/vector specified by SELECT,
DIF stores a Frobenius norm-based estimate of
Difl. If JOB = 'E', DIF is not referenced.
MM (input)
The number of elements in the arrays S and DIF. MM
>= M.
M (output)
The number of elements of the arrays S and DIF
used to store the specified condition numbers; for
each selected eigenvalue one element is used. If
HOWMNT = 'A', M is set to N.
WORK (workspace)
If JOB = 'E', WORK is not referenced. Otherwise,
on exit, if INFO = 0, WORK(1) returns the optimal
LWORK.
LWORK (input)
The dimension of the array WORK. LWORK >= 1. If
JOB = 'V' or 'B', LWORK >= 2*N*N.
IWORK (workspace)
dimension(N+2) If JOB = 'E', IWORK is not refer-
enced.
INFO (output)
= 0: Successful exit
< 0: If INFO = -i, the i-th argument had an ille-
gal value
The reciprocal of the condition number of the i-th general-
ized 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 eigen-
vector u and left eigenvector v corresponding to the gen-
eralized 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 pro-
duct 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 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.
[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.