1 SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
2 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER HOWMNY, JOB
12 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
13 * ..
14 * .. Array Arguments ..
15 LOGICAL SELECT( * )
16 INTEGER IWORK( * )
17 DOUBLE PRECISION DIF( * ), S( * )
18 COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
19 $ VR( LDVR, * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZTGSNA estimates reciprocal condition numbers for specified
26 * eigenvalues and/or eigenvectors of a matrix pair (A, B).
27 *
28 * (A, B) must be in generalized Schur canonical form, that is, A and
29 * B are both upper triangular.
30 *
31 * Arguments
32 * =========
33 *
34 * JOB (input) CHARACTER*1
35 * Specifies whether condition numbers are required for
36 * eigenvalues (S) or eigenvectors (DIF):
37 * = 'E': for eigenvalues only (S);
38 * = 'V': for eigenvectors only (DIF);
39 * = 'B': for both eigenvalues and eigenvectors (S and DIF).
40 *
41 * HOWMNY (input) CHARACTER*1
42 * = 'A': compute condition numbers for all eigenpairs;
43 * = 'S': compute condition numbers for selected eigenpairs
44 * specified by the array SELECT.
45 *
46 * SELECT (input) LOGICAL array, dimension (N)
47 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
48 * condition numbers are required. To select condition numbers
49 * for the corresponding j-th eigenvalue and/or eigenvector,
50 * SELECT(j) must be set to .TRUE..
51 * If HOWMNY = 'A', SELECT is not referenced.
52 *
53 * N (input) INTEGER
54 * The order of the square matrix pair (A, B). N >= 0.
55 *
56 * A (input) COMPLEX*16 array, dimension (LDA,N)
57 * The upper triangular matrix A in the pair (A,B).
58 *
59 * LDA (input) INTEGER
60 * The leading dimension of the array A. LDA >= max(1,N).
61 *
62 * B (input) COMPLEX*16 array, dimension (LDB,N)
63 * The upper triangular matrix B in the pair (A, B).
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * VL (input) COMPLEX*16 array, dimension (LDVL,M)
69 * IF JOB = 'E' or 'B', VL must contain left eigenvectors of
70 * (A, B), corresponding to the eigenpairs specified by HOWMNY
71 * and SELECT. The eigenvectors must be stored in consecutive
72 * columns of VL, as returned by ZTGEVC.
73 * If JOB = 'V', VL is not referenced.
74 *
75 * LDVL (input) INTEGER
76 * The leading dimension of the array VL. LDVL >= 1; and
77 * If JOB = 'E' or 'B', LDVL >= N.
78 *
79 * VR (input) COMPLEX*16 array, dimension (LDVR,M)
80 * IF JOB = 'E' or 'B', VR must contain right eigenvectors of
81 * (A, B), corresponding to the eigenpairs specified by HOWMNY
82 * and SELECT. The eigenvectors must be stored in consecutive
83 * columns of VR, as returned by ZTGEVC.
84 * If JOB = 'V', VR is not referenced.
85 *
86 * LDVR (input) INTEGER
87 * The leading dimension of the array VR. LDVR >= 1;
88 * If JOB = 'E' or 'B', LDVR >= N.
89 *
90 * S (output) DOUBLE PRECISION array, dimension (MM)
91 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
92 * selected eigenvalues, stored in consecutive elements of the
93 * array.
94 * If JOB = 'V', S is not referenced.
95 *
96 * DIF (output) DOUBLE PRECISION array, dimension (MM)
97 * If JOB = 'V' or 'B', the estimated reciprocal condition
98 * numbers of the selected eigenvectors, stored in consecutive
99 * elements of the array.
100 * If the eigenvalues cannot be reordered to compute DIF(j),
101 * DIF(j) is set to 0; this can only occur when the true value
102 * would be very small anyway.
103 * For each eigenvalue/vector specified by SELECT, DIF stores
104 * a Frobenius norm-based estimate of Difl.
105 * If JOB = 'E', DIF is not referenced.
106 *
107 * MM (input) INTEGER
108 * The number of elements in the arrays S and DIF. MM >= M.
109 *
110 * M (output) INTEGER
111 * The number of elements of the arrays S and DIF used to store
112 * the specified condition numbers; for each selected eigenvalue
113 * one element is used. If HOWMNY = 'A', M is set to N.
114 *
115 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
116 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
117 *
118 * LWORK (input) INTEGER
119 * The dimension of the array WORK. LWORK >= max(1,N).
120 * If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
121 *
122 * IWORK (workspace) INTEGER array, dimension (N+2)
123 * If JOB = 'E', IWORK is not referenced.
124 *
125 * INFO (output) INTEGER
126 * = 0: Successful exit
127 * < 0: If INFO = -i, the i-th argument had an illegal value
128 *
129 * Further Details
130 * ===============
131 *
132 * The reciprocal of the condition number of the i-th generalized
133 * eigenvalue w = (a, b) is defined as
134 *
135 * S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
136 *
137 * where u and v are the right and left eigenvectors of (A, B)
138 * corresponding to w; |z| denotes the absolute value of the complex
139 * number, and norm(u) denotes the 2-norm of the vector u. The pair
140 * (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
141 * matrix pair (A, B). If both a and b equal zero, then (A,B) is
142 * singular and S(I) = -1 is returned.
143 *
144 * An approximate error bound on the chordal distance between the i-th
145 * computed generalized eigenvalue w and the corresponding exact
146 * eigenvalue lambda is
147 *
148 * chord(w, lambda) <= EPS * norm(A, B) / S(I),
149 *
150 * where EPS is the machine precision.
151 *
152 * The reciprocal of the condition number of the right eigenvector u
153 * and left eigenvector v corresponding to the generalized eigenvalue w
154 * is defined as follows. Suppose
155 *
156 * (A, B) = ( a * ) ( b * ) 1
157 * ( 0 A22 ),( 0 B22 ) n-1
158 * 1 n-1 1 n-1
159 *
160 * Then the reciprocal condition number DIF(I) is
161 *
162 * Difl[(a, b), (A22, B22)] = sigma-min( Zl )
163 *
164 * where sigma-min(Zl) denotes the smallest singular value of
165 *
166 * Zl = [ kron(a, In-1) -kron(1, A22) ]
167 * [ kron(b, In-1) -kron(1, B22) ].
168 *
169 * Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
170 * transpose of X. kron(X, Y) is the Kronecker product between the
171 * matrices X and Y.
172 *
173 * We approximate the smallest singular value of Zl with an upper
174 * bound. This is done by ZLATDF.
175 *
176 * An approximate error bound for a computed eigenvector VL(i) or
177 * VR(i) is given by
178 *
179 * EPS * norm(A, B) / DIF(i).
180 *
181 * See ref. [2-3] for more details and further references.
182 *
183 * Based on contributions by
184 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
185 * Umea University, S-901 87 Umea, Sweden.
186 *
187 * References
188 * ==========
189 *
190 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
191 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
192 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
193 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
194 *
195 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
196 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
197 * Estimation: Theory, Algorithms and Software, Report
198 * UMINF - 94.04, Department of Computing Science, Umea University,
199 * S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
200 * To appear in Numerical Algorithms, 1996.
201 *
202 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
203 * for Solving the Generalized Sylvester Equation and Estimating the
204 * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
205 * Department of Computing Science, Umea University, S-901 87 Umea,
206 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
207 * Note 75.
208 * To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
209 *
210 * =====================================================================
211 *
212 * .. Parameters ..
213 DOUBLE PRECISION ZERO, ONE
214 INTEGER IDIFJB
215 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
216 * ..
217 * .. Local Scalars ..
218 LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS
219 INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
220 DOUBLE PRECISION BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
221 COMPLEX*16 YHAX, YHBX
222 * ..
223 * .. Local Arrays ..
224 COMPLEX*16 DUMMY( 1 ), DUMMY1( 1 )
225 * ..
226 * .. External Functions ..
227 LOGICAL LSAME
228 DOUBLE PRECISION DLAMCH, DLAPY2, DZNRM2
229 COMPLEX*16 ZDOTC
230 EXTERNAL LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
231 * ..
232 * .. External Subroutines ..
233 EXTERNAL DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
234 * ..
235 * .. Intrinsic Functions ..
236 INTRINSIC ABS, DCMPLX, MAX
237 * ..
238 * .. Executable Statements ..
239 *
240 * Decode and test the input parameters
241 *
242 WANTBH = LSAME( JOB, 'B' )
243 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
244 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
245 *
246 SOMCON = LSAME( HOWMNY, 'S' )
247 *
248 INFO = 0
249 LQUERY = ( LWORK.EQ.-1 )
250 *
251 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
252 INFO = -1
253 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
254 INFO = -2
255 ELSE IF( N.LT.0 ) THEN
256 INFO = -4
257 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
258 INFO = -6
259 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
260 INFO = -8
261 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
262 INFO = -10
263 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
264 INFO = -12
265 ELSE
266 *
267 * Set M to the number of eigenpairs for which condition numbers
268 * are required, and test MM.
269 *
270 IF( SOMCON ) THEN
271 M = 0
272 DO 10 K = 1, N
273 IF( SELECT( K ) )
274 $ M = M + 1
275 10 CONTINUE
276 ELSE
277 M = N
278 END IF
279 *
280 IF( N.EQ.0 ) THEN
281 LWMIN = 1
282 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
283 LWMIN = 2*N*N
284 ELSE
285 LWMIN = N
286 END IF
287 WORK( 1 ) = LWMIN
288 *
289 IF( MM.LT.M ) THEN
290 INFO = -15
291 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
292 INFO = -18
293 END IF
294 END IF
295 *
296 IF( INFO.NE.0 ) THEN
297 CALL XERBLA( 'ZTGSNA', -INFO )
298 RETURN
299 ELSE IF( LQUERY ) THEN
300 RETURN
301 END IF
302 *
303 * Quick return if possible
304 *
305 IF( N.EQ.0 )
306 $ RETURN
307 *
308 * Get machine constants
309 *
310 EPS = DLAMCH( 'P' )
311 SMLNUM = DLAMCH( 'S' ) / EPS
312 BIGNUM = ONE / SMLNUM
313 CALL DLABAD( SMLNUM, BIGNUM )
314 KS = 0
315 DO 20 K = 1, N
316 *
317 * Determine whether condition numbers are required for the k-th
318 * eigenpair.
319 *
320 IF( SOMCON ) THEN
321 IF( .NOT.SELECT( K ) )
322 $ GO TO 20
323 END IF
324 *
325 KS = KS + 1
326 *
327 IF( WANTS ) THEN
328 *
329 * Compute the reciprocal condition number of the k-th
330 * eigenvalue.
331 *
332 RNRM = DZNRM2( N, VR( 1, KS ), 1 )
333 LNRM = DZNRM2( N, VL( 1, KS ), 1 )
334 CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
335 $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
336 YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
337 CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
338 $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
339 YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
340 COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
341 IF( COND.EQ.ZERO ) THEN
342 S( KS ) = -ONE
343 ELSE
344 S( KS ) = COND / ( RNRM*LNRM )
345 END IF
346 END IF
347 *
348 IF( WANTDF ) THEN
349 IF( N.EQ.1 ) THEN
350 DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
351 ELSE
352 *
353 * Estimate the reciprocal condition number of the k-th
354 * eigenvectors.
355 *
356 * Copy the matrix (A, B) to the array WORK and move the
357 * (k,k)th pair to the (1,1) position.
358 *
359 CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
360 CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
361 IFST = K
362 ILST = 1
363 *
364 CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
365 $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
366 *
367 IF( IERR.GT.0 ) THEN
368 *
369 * Ill-conditioned problem - swap rejected.
370 *
371 DIF( KS ) = ZERO
372 ELSE
373 *
374 * Reordering successful, solve generalized Sylvester
375 * equation for R and L,
376 * A22 * R - L * A11 = A12
377 * B22 * R - L * B11 = B12,
378 * and compute estimate of Difl[(A11,B11), (A22, B22)].
379 *
380 N1 = 1
381 N2 = N - N1
382 I = N*N + 1
383 CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
384 $ N, WORK, N, WORK( N1+1 ), N,
385 $ WORK( N*N1+N1+I ), N, WORK( I ), N,
386 $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
387 $ 1, IWORK, IERR )
388 END IF
389 END IF
390 END IF
391 *
392 20 CONTINUE
393 WORK( 1 ) = LWMIN
394 RETURN
395 *
396 * End of ZTGSNA
397 *
398 END
2 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER HOWMNY, JOB
12 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
13 * ..
14 * .. Array Arguments ..
15 LOGICAL SELECT( * )
16 INTEGER IWORK( * )
17 DOUBLE PRECISION DIF( * ), S( * )
18 COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
19 $ VR( LDVR, * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZTGSNA estimates reciprocal condition numbers for specified
26 * eigenvalues and/or eigenvectors of a matrix pair (A, B).
27 *
28 * (A, B) must be in generalized Schur canonical form, that is, A and
29 * B are both upper triangular.
30 *
31 * Arguments
32 * =========
33 *
34 * JOB (input) CHARACTER*1
35 * Specifies whether condition numbers are required for
36 * eigenvalues (S) or eigenvectors (DIF):
37 * = 'E': for eigenvalues only (S);
38 * = 'V': for eigenvectors only (DIF);
39 * = 'B': for both eigenvalues and eigenvectors (S and DIF).
40 *
41 * HOWMNY (input) CHARACTER*1
42 * = 'A': compute condition numbers for all eigenpairs;
43 * = 'S': compute condition numbers for selected eigenpairs
44 * specified by the array SELECT.
45 *
46 * SELECT (input) LOGICAL array, dimension (N)
47 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
48 * condition numbers are required. To select condition numbers
49 * for the corresponding j-th eigenvalue and/or eigenvector,
50 * SELECT(j) must be set to .TRUE..
51 * If HOWMNY = 'A', SELECT is not referenced.
52 *
53 * N (input) INTEGER
54 * The order of the square matrix pair (A, B). N >= 0.
55 *
56 * A (input) COMPLEX*16 array, dimension (LDA,N)
57 * The upper triangular matrix A in the pair (A,B).
58 *
59 * LDA (input) INTEGER
60 * The leading dimension of the array A. LDA >= max(1,N).
61 *
62 * B (input) COMPLEX*16 array, dimension (LDB,N)
63 * The upper triangular matrix B in the pair (A, B).
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * VL (input) COMPLEX*16 array, dimension (LDVL,M)
69 * IF JOB = 'E' or 'B', VL must contain left eigenvectors of
70 * (A, B), corresponding to the eigenpairs specified by HOWMNY
71 * and SELECT. The eigenvectors must be stored in consecutive
72 * columns of VL, as returned by ZTGEVC.
73 * If JOB = 'V', VL is not referenced.
74 *
75 * LDVL (input) INTEGER
76 * The leading dimension of the array VL. LDVL >= 1; and
77 * If JOB = 'E' or 'B', LDVL >= N.
78 *
79 * VR (input) COMPLEX*16 array, dimension (LDVR,M)
80 * IF JOB = 'E' or 'B', VR must contain right eigenvectors of
81 * (A, B), corresponding to the eigenpairs specified by HOWMNY
82 * and SELECT. The eigenvectors must be stored in consecutive
83 * columns of VR, as returned by ZTGEVC.
84 * If JOB = 'V', VR is not referenced.
85 *
86 * LDVR (input) INTEGER
87 * The leading dimension of the array VR. LDVR >= 1;
88 * If JOB = 'E' or 'B', LDVR >= N.
89 *
90 * S (output) DOUBLE PRECISION array, dimension (MM)
91 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
92 * selected eigenvalues, stored in consecutive elements of the
93 * array.
94 * If JOB = 'V', S is not referenced.
95 *
96 * DIF (output) DOUBLE PRECISION array, dimension (MM)
97 * If JOB = 'V' or 'B', the estimated reciprocal condition
98 * numbers of the selected eigenvectors, stored in consecutive
99 * elements of the array.
100 * If the eigenvalues cannot be reordered to compute DIF(j),
101 * DIF(j) is set to 0; this can only occur when the true value
102 * would be very small anyway.
103 * For each eigenvalue/vector specified by SELECT, DIF stores
104 * a Frobenius norm-based estimate of Difl.
105 * If JOB = 'E', DIF is not referenced.
106 *
107 * MM (input) INTEGER
108 * The number of elements in the arrays S and DIF. MM >= M.
109 *
110 * M (output) INTEGER
111 * The number of elements of the arrays S and DIF used to store
112 * the specified condition numbers; for each selected eigenvalue
113 * one element is used. If HOWMNY = 'A', M is set to N.
114 *
115 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
116 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
117 *
118 * LWORK (input) INTEGER
119 * The dimension of the array WORK. LWORK >= max(1,N).
120 * If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
121 *
122 * IWORK (workspace) INTEGER array, dimension (N+2)
123 * If JOB = 'E', IWORK is not referenced.
124 *
125 * INFO (output) INTEGER
126 * = 0: Successful exit
127 * < 0: If INFO = -i, the i-th argument had an illegal value
128 *
129 * Further Details
130 * ===============
131 *
132 * The reciprocal of the condition number of the i-th generalized
133 * eigenvalue w = (a, b) is defined as
134 *
135 * S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
136 *
137 * where u and v are the right and left eigenvectors of (A, B)
138 * corresponding to w; |z| denotes the absolute value of the complex
139 * number, and norm(u) denotes the 2-norm of the vector u. The pair
140 * (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
141 * matrix pair (A, B). If both a and b equal zero, then (A,B) is
142 * singular and S(I) = -1 is returned.
143 *
144 * An approximate error bound on the chordal distance between the i-th
145 * computed generalized eigenvalue w and the corresponding exact
146 * eigenvalue lambda is
147 *
148 * chord(w, lambda) <= EPS * norm(A, B) / S(I),
149 *
150 * where EPS is the machine precision.
151 *
152 * The reciprocal of the condition number of the right eigenvector u
153 * and left eigenvector v corresponding to the generalized eigenvalue w
154 * is defined as follows. Suppose
155 *
156 * (A, B) = ( a * ) ( b * ) 1
157 * ( 0 A22 ),( 0 B22 ) n-1
158 * 1 n-1 1 n-1
159 *
160 * Then the reciprocal condition number DIF(I) is
161 *
162 * Difl[(a, b), (A22, B22)] = sigma-min( Zl )
163 *
164 * where sigma-min(Zl) denotes the smallest singular value of
165 *
166 * Zl = [ kron(a, In-1) -kron(1, A22) ]
167 * [ kron(b, In-1) -kron(1, B22) ].
168 *
169 * Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
170 * transpose of X. kron(X, Y) is the Kronecker product between the
171 * matrices X and Y.
172 *
173 * We approximate the smallest singular value of Zl with an upper
174 * bound. This is done by ZLATDF.
175 *
176 * An approximate error bound for a computed eigenvector VL(i) or
177 * VR(i) is given by
178 *
179 * EPS * norm(A, B) / DIF(i).
180 *
181 * See ref. [2-3] for more details and further references.
182 *
183 * Based on contributions by
184 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
185 * Umea University, S-901 87 Umea, Sweden.
186 *
187 * References
188 * ==========
189 *
190 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
191 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
192 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
193 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
194 *
195 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
196 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
197 * Estimation: Theory, Algorithms and Software, Report
198 * UMINF - 94.04, Department of Computing Science, Umea University,
199 * S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
200 * To appear in Numerical Algorithms, 1996.
201 *
202 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
203 * for Solving the Generalized Sylvester Equation and Estimating the
204 * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
205 * Department of Computing Science, Umea University, S-901 87 Umea,
206 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
207 * Note 75.
208 * To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
209 *
210 * =====================================================================
211 *
212 * .. Parameters ..
213 DOUBLE PRECISION ZERO, ONE
214 INTEGER IDIFJB
215 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
216 * ..
217 * .. Local Scalars ..
218 LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS
219 INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
220 DOUBLE PRECISION BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
221 COMPLEX*16 YHAX, YHBX
222 * ..
223 * .. Local Arrays ..
224 COMPLEX*16 DUMMY( 1 ), DUMMY1( 1 )
225 * ..
226 * .. External Functions ..
227 LOGICAL LSAME
228 DOUBLE PRECISION DLAMCH, DLAPY2, DZNRM2
229 COMPLEX*16 ZDOTC
230 EXTERNAL LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
231 * ..
232 * .. External Subroutines ..
233 EXTERNAL DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
234 * ..
235 * .. Intrinsic Functions ..
236 INTRINSIC ABS, DCMPLX, MAX
237 * ..
238 * .. Executable Statements ..
239 *
240 * Decode and test the input parameters
241 *
242 WANTBH = LSAME( JOB, 'B' )
243 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
244 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
245 *
246 SOMCON = LSAME( HOWMNY, 'S' )
247 *
248 INFO = 0
249 LQUERY = ( LWORK.EQ.-1 )
250 *
251 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
252 INFO = -1
253 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
254 INFO = -2
255 ELSE IF( N.LT.0 ) THEN
256 INFO = -4
257 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
258 INFO = -6
259 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
260 INFO = -8
261 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
262 INFO = -10
263 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
264 INFO = -12
265 ELSE
266 *
267 * Set M to the number of eigenpairs for which condition numbers
268 * are required, and test MM.
269 *
270 IF( SOMCON ) THEN
271 M = 0
272 DO 10 K = 1, N
273 IF( SELECT( K ) )
274 $ M = M + 1
275 10 CONTINUE
276 ELSE
277 M = N
278 END IF
279 *
280 IF( N.EQ.0 ) THEN
281 LWMIN = 1
282 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
283 LWMIN = 2*N*N
284 ELSE
285 LWMIN = N
286 END IF
287 WORK( 1 ) = LWMIN
288 *
289 IF( MM.LT.M ) THEN
290 INFO = -15
291 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
292 INFO = -18
293 END IF
294 END IF
295 *
296 IF( INFO.NE.0 ) THEN
297 CALL XERBLA( 'ZTGSNA', -INFO )
298 RETURN
299 ELSE IF( LQUERY ) THEN
300 RETURN
301 END IF
302 *
303 * Quick return if possible
304 *
305 IF( N.EQ.0 )
306 $ RETURN
307 *
308 * Get machine constants
309 *
310 EPS = DLAMCH( 'P' )
311 SMLNUM = DLAMCH( 'S' ) / EPS
312 BIGNUM = ONE / SMLNUM
313 CALL DLABAD( SMLNUM, BIGNUM )
314 KS = 0
315 DO 20 K = 1, N
316 *
317 * Determine whether condition numbers are required for the k-th
318 * eigenpair.
319 *
320 IF( SOMCON ) THEN
321 IF( .NOT.SELECT( K ) )
322 $ GO TO 20
323 END IF
324 *
325 KS = KS + 1
326 *
327 IF( WANTS ) THEN
328 *
329 * Compute the reciprocal condition number of the k-th
330 * eigenvalue.
331 *
332 RNRM = DZNRM2( N, VR( 1, KS ), 1 )
333 LNRM = DZNRM2( N, VL( 1, KS ), 1 )
334 CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
335 $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
336 YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
337 CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
338 $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
339 YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
340 COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
341 IF( COND.EQ.ZERO ) THEN
342 S( KS ) = -ONE
343 ELSE
344 S( KS ) = COND / ( RNRM*LNRM )
345 END IF
346 END IF
347 *
348 IF( WANTDF ) THEN
349 IF( N.EQ.1 ) THEN
350 DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
351 ELSE
352 *
353 * Estimate the reciprocal condition number of the k-th
354 * eigenvectors.
355 *
356 * Copy the matrix (A, B) to the array WORK and move the
357 * (k,k)th pair to the (1,1) position.
358 *
359 CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
360 CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
361 IFST = K
362 ILST = 1
363 *
364 CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
365 $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
366 *
367 IF( IERR.GT.0 ) THEN
368 *
369 * Ill-conditioned problem - swap rejected.
370 *
371 DIF( KS ) = ZERO
372 ELSE
373 *
374 * Reordering successful, solve generalized Sylvester
375 * equation for R and L,
376 * A22 * R - L * A11 = A12
377 * B22 * R - L * B11 = B12,
378 * and compute estimate of Difl[(A11,B11), (A22, B22)].
379 *
380 N1 = 1
381 N2 = N - N1
382 I = N*N + 1
383 CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
384 $ N, WORK, N, WORK( N1+1 ), N,
385 $ WORK( N*N1+N1+I ), N, WORK( I ), N,
386 $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
387 $ 1, IWORK, IERR )
388 END IF
389 END IF
390 END IF
391 *
392 20 CONTINUE
393 WORK( 1 ) = LWMIN
394 RETURN
395 *
396 * End of ZTGSNA
397 *
398 END