1 SUBROUTINE DTGSNA( 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 A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
18 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DTGSNA estimates reciprocal condition numbers for specified
25 * eigenvalues and/or eigenvectors of a matrix pair (A, B) in
26 * generalized real Schur canonical form (or of any matrix pair
27 * (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
28 * Z**T denotes the transpose of Z.
29 *
30 * (A, B) must be in generalized real Schur form (as returned by DGGES),
31 * i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
32 * blocks. B is upper triangular.
33 *
34 *
35 * Arguments
36 * =========
37 *
38 * JOB (input) CHARACTER*1
39 * Specifies whether condition numbers are required for
40 * eigenvalues (S) or eigenvectors (DIF):
41 * = 'E': for eigenvalues only (S);
42 * = 'V': for eigenvectors only (DIF);
43 * = 'B': for both eigenvalues and eigenvectors (S and DIF).
44 *
45 * HOWMNY (input) CHARACTER*1
46 * = 'A': compute condition numbers for all eigenpairs;
47 * = 'S': compute condition numbers for selected eigenpairs
48 * specified by the array SELECT.
49 *
50 * SELECT (input) LOGICAL array, dimension (N)
51 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
52 * condition numbers are required. To select condition numbers
53 * for the eigenpair corresponding to a real eigenvalue w(j),
54 * SELECT(j) must be set to .TRUE.. To select condition numbers
55 * corresponding to a complex conjugate pair of eigenvalues w(j)
56 * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
57 * set to .TRUE..
58 * If HOWMNY = 'A', SELECT is not referenced.
59 *
60 * N (input) INTEGER
61 * The order of the square matrix pair (A, B). N >= 0.
62 *
63 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
64 * The upper quasi-triangular matrix A in the pair (A,B).
65 *
66 * LDA (input) INTEGER
67 * The leading dimension of the array A. LDA >= max(1,N).
68 *
69 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
70 * The upper triangular matrix B in the pair (A,B).
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the array B. LDB >= max(1,N).
74 *
75 * VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
76 * If JOB = 'E' or 'B', VL must contain left eigenvectors of
77 * (A, B), corresponding to the eigenpairs specified by HOWMNY
78 * and SELECT. The eigenvectors must be stored in consecutive
79 * columns of VL, as returned by DTGEVC.
80 * If JOB = 'V', VL is not referenced.
81 *
82 * LDVL (input) INTEGER
83 * The leading dimension of the array VL. LDVL >= 1.
84 * If JOB = 'E' or 'B', LDVL >= N.
85 *
86 * VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
87 * If JOB = 'E' or 'B', VR must contain right eigenvectors of
88 * (A, B), corresponding to the eigenpairs specified by HOWMNY
89 * and SELECT. The eigenvectors must be stored in consecutive
90 * columns ov VR, as returned by DTGEVC.
91 * If JOB = 'V', VR is not referenced.
92 *
93 * LDVR (input) INTEGER
94 * The leading dimension of the array VR. LDVR >= 1.
95 * If JOB = 'E' or 'B', LDVR >= N.
96 *
97 * S (output) DOUBLE PRECISION array, dimension (MM)
98 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
99 * selected eigenvalues, stored in consecutive elements of the
100 * array. For a complex conjugate pair of eigenvalues two
101 * consecutive elements of S are set to the same value. Thus
102 * S(j), DIF(j), and the j-th columns of VL and VR all
103 * correspond to the same eigenpair (but not in general the
104 * j-th eigenpair, unless all eigenpairs are selected).
105 * If JOB = 'V', S is not referenced.
106 *
107 * DIF (output) DOUBLE PRECISION array, dimension (MM)
108 * If JOB = 'V' or 'B', the estimated reciprocal condition
109 * numbers of the selected eigenvectors, stored in consecutive
110 * elements of the array. For a complex eigenvector two
111 * consecutive elements of DIF are set to the same value. If
112 * the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
113 * is set to 0; this can only occur when the true value would be
114 * very small anyway.
115 * If JOB = 'E', DIF is not referenced.
116 *
117 * MM (input) INTEGER
118 * The number of elements in the arrays S and DIF. MM >= M.
119 *
120 * M (output) INTEGER
121 * The number of elements of the arrays S and DIF used to store
122 * the specified condition numbers; for each selected real
123 * eigenvalue one element is used, and for each selected complex
124 * conjugate pair of eigenvalues, two elements are used.
125 * If HOWMNY = 'A', M is set to N.
126 *
127 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
128 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129 *
130 * LWORK (input) INTEGER
131 * The dimension of the array WORK. LWORK >= max(1,N).
132 * If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
133 *
134 * If LWORK = -1, then a workspace query is assumed; the routine
135 * only calculates the optimal size of the WORK array, returns
136 * this value as the first entry of the WORK array, and no error
137 * message related to LWORK is issued by XERBLA.
138 *
139 * IWORK (workspace) INTEGER array, dimension (N + 6)
140 * If JOB = 'E', IWORK is not referenced.
141 *
142 * INFO (output) INTEGER
143 * =0: Successful exit
144 * <0: If INFO = -i, the i-th argument had an illegal value
145 *
146 *
147 * Further Details
148 * ===============
149 *
150 * The reciprocal of the condition number of a generalized eigenvalue
151 * w = (a, b) is defined as
152 *
153 * S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
154 *
155 * where u and v are the left and right eigenvectors of (A, B)
156 * corresponding to w; |z| denotes the absolute value of the complex
157 * number, and norm(u) denotes the 2-norm of the vector u.
158 * The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
159 * of the matrix pair (A, B). If both a and b equal zero, then (A B) is
160 * singular and S(I) = -1 is returned.
161 *
162 * An approximate error bound on the chordal distance between the i-th
163 * computed generalized eigenvalue w and the corresponding exact
164 * eigenvalue lambda is
165 *
166 * chord(w, lambda) <= EPS * norm(A, B) / S(I)
167 *
168 * where EPS is the machine precision.
169 *
170 * The reciprocal of the condition number DIF(i) of right eigenvector u
171 * and left eigenvector v corresponding to the generalized eigenvalue w
172 * is defined as follows:
173 *
174 * a) If the i-th eigenvalue w = (a,b) is real
175 *
176 * Suppose U and V are orthogonal transformations such that
177 *
178 * U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
179 * ( 0 S22 ),( 0 T22 ) n-1
180 * 1 n-1 1 n-1
181 *
182 * Then the reciprocal condition number DIF(i) is
183 *
184 * Difl((a, b), (S22, T22)) = sigma-min( Zl ),
185 *
186 * where sigma-min(Zl) denotes the smallest singular value of the
187 * 2(n-1)-by-2(n-1) matrix
188 *
189 * Zl = [ kron(a, In-1) -kron(1, S22) ]
190 * [ kron(b, In-1) -kron(1, T22) ] .
191 *
192 * Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
193 * Kronecker product between the matrices X and Y.
194 *
195 * Note that if the default method for computing DIF(i) is wanted
196 * (see DLATDF), then the parameter DIFDRI (see below) should be
197 * changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
198 * See DTGSYL for more details.
199 *
200 * b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
201 *
202 * Suppose U and V are orthogonal transformations such that
203 *
204 * U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
205 * ( 0 S22 ),( 0 T22) n-2
206 * 2 n-2 2 n-2
207 *
208 * and (S11, T11) corresponds to the complex conjugate eigenvalue
209 * pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
210 * that
211 *
212 * U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
213 * ( 0 s22 ) ( 0 t22 )
214 *
215 * where the generalized eigenvalues w = s11/t11 and
216 * conjg(w) = s22/t22.
217 *
218 * Then the reciprocal condition number DIF(i) is bounded by
219 *
220 * min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
221 *
222 * where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
223 * Z1 is the complex 2-by-2 matrix
224 *
225 * Z1 = [ s11 -s22 ]
226 * [ t11 -t22 ],
227 *
228 * This is done by computing (using real arithmetic) the
229 * roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
230 * where Z1**T denotes the transpose of Z1 and det(X) denotes
231 * the determinant of X.
232 *
233 * and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
234 * upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
235 *
236 * Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ]
237 * [ kron(T11**T, In-2) -kron(I2, T22) ]
238 *
239 * Note that if the default method for computing DIF is wanted (see
240 * DLATDF), then the parameter DIFDRI (see below) should be changed
241 * from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
242 * for more details.
243 *
244 * For each eigenvalue/vector specified by SELECT, DIF stores a
245 * Frobenius norm-based estimate of Difl.
246 *
247 * An approximate error bound for the i-th computed eigenvector VL(i) or
248 * VR(i) is given by
249 *
250 * EPS * norm(A, B) / DIF(i).
251 *
252 * See ref. [2-3] for more details and further references.
253 *
254 * Based on contributions by
255 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
256 * Umea University, S-901 87 Umea, Sweden.
257 *
258 * References
259 * ==========
260 *
261 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
262 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
263 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
264 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
265 *
266 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
267 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
268 * Estimation: Theory, Algorithms and Software,
269 * Report UMINF - 94.04, Department of Computing Science, Umea
270 * University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
271 * Note 87. To appear in Numerical Algorithms, 1996.
272 *
273 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
274 * for Solving the Generalized Sylvester Equation and Estimating the
275 * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
276 * Department of Computing Science, Umea University, S-901 87 Umea,
277 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
278 * Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
279 * No 1, 1996.
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284 INTEGER DIFDRI
285 PARAMETER ( DIFDRI = 3 )
286 DOUBLE PRECISION ZERO, ONE, TWO, FOUR
287 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
288 $ FOUR = 4.0D+0 )
289 * ..
290 * .. Local Scalars ..
291 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
292 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
293 DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
294 $ EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
295 $ TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
296 $ UHBVI
297 * ..
298 * .. Local Arrays ..
299 DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 )
300 * ..
301 * .. External Functions ..
302 LOGICAL LSAME
303 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
304 EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
305 * ..
306 * .. External Subroutines ..
307 EXTERNAL DGEMV, DLACPY, DLAG2, DTGEXC, DTGSYL, XERBLA
308 * ..
309 * .. Intrinsic Functions ..
310 INTRINSIC MAX, MIN, SQRT
311 * ..
312 * .. Executable Statements ..
313 *
314 * Decode and test the input parameters
315 *
316 WANTBH = LSAME( JOB, 'B' )
317 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
318 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
319 *
320 SOMCON = LSAME( HOWMNY, 'S' )
321 *
322 INFO = 0
323 LQUERY = ( LWORK.EQ.-1 )
324 *
325 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
326 INFO = -1
327 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
328 INFO = -2
329 ELSE IF( N.LT.0 ) THEN
330 INFO = -4
331 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
332 INFO = -6
333 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
334 INFO = -8
335 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
336 INFO = -10
337 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
338 INFO = -12
339 ELSE
340 *
341 * Set M to the number of eigenpairs for which condition numbers
342 * are required, and test MM.
343 *
344 IF( SOMCON ) THEN
345 M = 0
346 PAIR = .FALSE.
347 DO 10 K = 1, N
348 IF( PAIR ) THEN
349 PAIR = .FALSE.
350 ELSE
351 IF( K.LT.N ) THEN
352 IF( A( K+1, K ).EQ.ZERO ) THEN
353 IF( SELECT( K ) )
354 $ M = M + 1
355 ELSE
356 PAIR = .TRUE.
357 IF( SELECT( K ) .OR. SELECT( K+1 ) )
358 $ M = M + 2
359 END IF
360 ELSE
361 IF( SELECT( N ) )
362 $ M = M + 1
363 END IF
364 END IF
365 10 CONTINUE
366 ELSE
367 M = N
368 END IF
369 *
370 IF( N.EQ.0 ) THEN
371 LWMIN = 1
372 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
373 LWMIN = 2*N*( N + 2 ) + 16
374 ELSE
375 LWMIN = N
376 END IF
377 WORK( 1 ) = LWMIN
378 *
379 IF( MM.LT.M ) THEN
380 INFO = -15
381 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
382 INFO = -18
383 END IF
384 END IF
385 *
386 IF( INFO.NE.0 ) THEN
387 CALL XERBLA( 'DTGSNA', -INFO )
388 RETURN
389 ELSE IF( LQUERY ) THEN
390 RETURN
391 END IF
392 *
393 * Quick return if possible
394 *
395 IF( N.EQ.0 )
396 $ RETURN
397 *
398 * Get machine constants
399 *
400 EPS = DLAMCH( 'P' )
401 SMLNUM = DLAMCH( 'S' ) / EPS
402 KS = 0
403 PAIR = .FALSE.
404 *
405 DO 20 K = 1, N
406 *
407 * Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
408 *
409 IF( PAIR ) THEN
410 PAIR = .FALSE.
411 GO TO 20
412 ELSE
413 IF( K.LT.N )
414 $ PAIR = A( K+1, K ).NE.ZERO
415 END IF
416 *
417 * Determine whether condition numbers are required for the k-th
418 * eigenpair.
419 *
420 IF( SOMCON ) THEN
421 IF( PAIR ) THEN
422 IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
423 $ GO TO 20
424 ELSE
425 IF( .NOT.SELECT( K ) )
426 $ GO TO 20
427 END IF
428 END IF
429 *
430 KS = KS + 1
431 *
432 IF( WANTS ) THEN
433 *
434 * Compute the reciprocal condition number of the k-th
435 * eigenvalue.
436 *
437 IF( PAIR ) THEN
438 *
439 * Complex eigenvalue pair.
440 *
441 RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
442 $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
443 LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
444 $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
445 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
446 $ WORK, 1 )
447 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
448 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
449 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
450 $ ZERO, WORK, 1 )
451 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
452 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
453 UHAV = TMPRR + TMPII
454 UHAVI = TMPIR - TMPRI
455 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
456 $ WORK, 1 )
457 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
458 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
459 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
460 $ ZERO, WORK, 1 )
461 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
462 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
463 UHBV = TMPRR + TMPII
464 UHBVI = TMPIR - TMPRI
465 UHAV = DLAPY2( UHAV, UHAVI )
466 UHBV = DLAPY2( UHBV, UHBVI )
467 COND = DLAPY2( UHAV, UHBV )
468 S( KS ) = COND / ( RNRM*LNRM )
469 S( KS+1 ) = S( KS )
470 *
471 ELSE
472 *
473 * Real eigenvalue.
474 *
475 RNRM = DNRM2( N, VR( 1, KS ), 1 )
476 LNRM = DNRM2( N, VL( 1, KS ), 1 )
477 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
478 $ WORK, 1 )
479 UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
480 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
481 $ WORK, 1 )
482 UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
483 COND = DLAPY2( UHAV, UHBV )
484 IF( COND.EQ.ZERO ) THEN
485 S( KS ) = -ONE
486 ELSE
487 S( KS ) = COND / ( RNRM*LNRM )
488 END IF
489 END IF
490 END IF
491 *
492 IF( WANTDF ) THEN
493 IF( N.EQ.1 ) THEN
494 DIF( KS ) = DLAPY2( A( 1, 1 ), B( 1, 1 ) )
495 GO TO 20
496 END IF
497 *
498 * Estimate the reciprocal condition number of the k-th
499 * eigenvectors.
500 IF( PAIR ) THEN
501 *
502 * Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)).
503 * Compute the eigenvalue(s) at position K.
504 *
505 WORK( 1 ) = A( K, K )
506 WORK( 2 ) = A( K+1, K )
507 WORK( 3 ) = A( K, K+1 )
508 WORK( 4 ) = A( K+1, K+1 )
509 WORK( 5 ) = B( K, K )
510 WORK( 6 ) = B( K+1, K )
511 WORK( 7 ) = B( K, K+1 )
512 WORK( 8 ) = B( K+1, K+1 )
513 CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
514 $ DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
515 ALPRQT = ONE
516 C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
517 C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
518 ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 )
519 ROOT2 = C2 / ROOT1
520 ROOT1 = ROOT1 / TWO
521 COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
522 END IF
523 *
524 * Copy the matrix (A, B) to the array WORK and swap the
525 * diagonal block beginning at A(k,k) to the (1,1) position.
526 *
527 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
528 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
529 IFST = K
530 ILST = 1
531 *
532 CALL DTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
533 $ DUMMY, 1, DUMMY1, 1, IFST, ILST,
534 $ WORK( N*N*2+1 ), LWORK-2*N*N, IERR )
535 *
536 IF( IERR.GT.0 ) THEN
537 *
538 * Ill-conditioned problem - swap rejected.
539 *
540 DIF( KS ) = ZERO
541 ELSE
542 *
543 * Reordering successful, solve generalized Sylvester
544 * equation for R and L,
545 * A22 * R - L * A11 = A12
546 * B22 * R - L * B11 = B12,
547 * and compute estimate of Difl((A11,B11), (A22, B22)).
548 *
549 N1 = 1
550 IF( WORK( 2 ).NE.ZERO )
551 $ N1 = 2
552 N2 = N - N1
553 IF( N2.EQ.0 ) THEN
554 DIF( KS ) = COND
555 ELSE
556 I = N*N + 1
557 IZ = 2*N*N + 1
558 CALL DTGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ),
559 $ N, WORK, N, WORK( N1+1 ), N,
560 $ WORK( N*N1+N1+I ), N, WORK( I ), N,
561 $ WORK( N1+I ), N, SCALE, DIF( KS ),
562 $ WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
563 *
564 IF( PAIR )
565 $ DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
566 $ COND )
567 END IF
568 END IF
569 IF( PAIR )
570 $ DIF( KS+1 ) = DIF( KS )
571 END IF
572 IF( PAIR )
573 $ KS = KS + 1
574 *
575 20 CONTINUE
576 WORK( 1 ) = LWMIN
577 RETURN
578 *
579 * End of DTGSNA
580 *
581 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 A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
18 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DTGSNA estimates reciprocal condition numbers for specified
25 * eigenvalues and/or eigenvectors of a matrix pair (A, B) in
26 * generalized real Schur canonical form (or of any matrix pair
27 * (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
28 * Z**T denotes the transpose of Z.
29 *
30 * (A, B) must be in generalized real Schur form (as returned by DGGES),
31 * i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
32 * blocks. B is upper triangular.
33 *
34 *
35 * Arguments
36 * =========
37 *
38 * JOB (input) CHARACTER*1
39 * Specifies whether condition numbers are required for
40 * eigenvalues (S) or eigenvectors (DIF):
41 * = 'E': for eigenvalues only (S);
42 * = 'V': for eigenvectors only (DIF);
43 * = 'B': for both eigenvalues and eigenvectors (S and DIF).
44 *
45 * HOWMNY (input) CHARACTER*1
46 * = 'A': compute condition numbers for all eigenpairs;
47 * = 'S': compute condition numbers for selected eigenpairs
48 * specified by the array SELECT.
49 *
50 * SELECT (input) LOGICAL array, dimension (N)
51 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
52 * condition numbers are required. To select condition numbers
53 * for the eigenpair corresponding to a real eigenvalue w(j),
54 * SELECT(j) must be set to .TRUE.. To select condition numbers
55 * corresponding to a complex conjugate pair of eigenvalues w(j)
56 * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
57 * set to .TRUE..
58 * If HOWMNY = 'A', SELECT is not referenced.
59 *
60 * N (input) INTEGER
61 * The order of the square matrix pair (A, B). N >= 0.
62 *
63 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
64 * The upper quasi-triangular matrix A in the pair (A,B).
65 *
66 * LDA (input) INTEGER
67 * The leading dimension of the array A. LDA >= max(1,N).
68 *
69 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
70 * The upper triangular matrix B in the pair (A,B).
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the array B. LDB >= max(1,N).
74 *
75 * VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
76 * If JOB = 'E' or 'B', VL must contain left eigenvectors of
77 * (A, B), corresponding to the eigenpairs specified by HOWMNY
78 * and SELECT. The eigenvectors must be stored in consecutive
79 * columns of VL, as returned by DTGEVC.
80 * If JOB = 'V', VL is not referenced.
81 *
82 * LDVL (input) INTEGER
83 * The leading dimension of the array VL. LDVL >= 1.
84 * If JOB = 'E' or 'B', LDVL >= N.
85 *
86 * VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
87 * If JOB = 'E' or 'B', VR must contain right eigenvectors of
88 * (A, B), corresponding to the eigenpairs specified by HOWMNY
89 * and SELECT. The eigenvectors must be stored in consecutive
90 * columns ov VR, as returned by DTGEVC.
91 * If JOB = 'V', VR is not referenced.
92 *
93 * LDVR (input) INTEGER
94 * The leading dimension of the array VR. LDVR >= 1.
95 * If JOB = 'E' or 'B', LDVR >= N.
96 *
97 * S (output) DOUBLE PRECISION array, dimension (MM)
98 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
99 * selected eigenvalues, stored in consecutive elements of the
100 * array. For a complex conjugate pair of eigenvalues two
101 * consecutive elements of S are set to the same value. Thus
102 * S(j), DIF(j), and the j-th columns of VL and VR all
103 * correspond to the same eigenpair (but not in general the
104 * j-th eigenpair, unless all eigenpairs are selected).
105 * If JOB = 'V', S is not referenced.
106 *
107 * DIF (output) DOUBLE PRECISION array, dimension (MM)
108 * If JOB = 'V' or 'B', the estimated reciprocal condition
109 * numbers of the selected eigenvectors, stored in consecutive
110 * elements of the array. For a complex eigenvector two
111 * consecutive elements of DIF are set to the same value. If
112 * the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
113 * is set to 0; this can only occur when the true value would be
114 * very small anyway.
115 * If JOB = 'E', DIF is not referenced.
116 *
117 * MM (input) INTEGER
118 * The number of elements in the arrays S and DIF. MM >= M.
119 *
120 * M (output) INTEGER
121 * The number of elements of the arrays S and DIF used to store
122 * the specified condition numbers; for each selected real
123 * eigenvalue one element is used, and for each selected complex
124 * conjugate pair of eigenvalues, two elements are used.
125 * If HOWMNY = 'A', M is set to N.
126 *
127 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
128 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129 *
130 * LWORK (input) INTEGER
131 * The dimension of the array WORK. LWORK >= max(1,N).
132 * If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
133 *
134 * If LWORK = -1, then a workspace query is assumed; the routine
135 * only calculates the optimal size of the WORK array, returns
136 * this value as the first entry of the WORK array, and no error
137 * message related to LWORK is issued by XERBLA.
138 *
139 * IWORK (workspace) INTEGER array, dimension (N + 6)
140 * If JOB = 'E', IWORK is not referenced.
141 *
142 * INFO (output) INTEGER
143 * =0: Successful exit
144 * <0: If INFO = -i, the i-th argument had an illegal value
145 *
146 *
147 * Further Details
148 * ===============
149 *
150 * The reciprocal of the condition number of a generalized eigenvalue
151 * w = (a, b) is defined as
152 *
153 * S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
154 *
155 * where u and v are the left and right eigenvectors of (A, B)
156 * corresponding to w; |z| denotes the absolute value of the complex
157 * number, and norm(u) denotes the 2-norm of the vector u.
158 * The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
159 * of the matrix pair (A, B). If both a and b equal zero, then (A B) is
160 * singular and S(I) = -1 is returned.
161 *
162 * An approximate error bound on the chordal distance between the i-th
163 * computed generalized eigenvalue w and the corresponding exact
164 * eigenvalue lambda is
165 *
166 * chord(w, lambda) <= EPS * norm(A, B) / S(I)
167 *
168 * where EPS is the machine precision.
169 *
170 * The reciprocal of the condition number DIF(i) of right eigenvector u
171 * and left eigenvector v corresponding to the generalized eigenvalue w
172 * is defined as follows:
173 *
174 * a) If the i-th eigenvalue w = (a,b) is real
175 *
176 * Suppose U and V are orthogonal transformations such that
177 *
178 * U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
179 * ( 0 S22 ),( 0 T22 ) n-1
180 * 1 n-1 1 n-1
181 *
182 * Then the reciprocal condition number DIF(i) is
183 *
184 * Difl((a, b), (S22, T22)) = sigma-min( Zl ),
185 *
186 * where sigma-min(Zl) denotes the smallest singular value of the
187 * 2(n-1)-by-2(n-1) matrix
188 *
189 * Zl = [ kron(a, In-1) -kron(1, S22) ]
190 * [ kron(b, In-1) -kron(1, T22) ] .
191 *
192 * Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
193 * Kronecker product between the matrices X and Y.
194 *
195 * Note that if the default method for computing DIF(i) is wanted
196 * (see DLATDF), then the parameter DIFDRI (see below) should be
197 * changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
198 * See DTGSYL for more details.
199 *
200 * b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
201 *
202 * Suppose U and V are orthogonal transformations such that
203 *
204 * U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
205 * ( 0 S22 ),( 0 T22) n-2
206 * 2 n-2 2 n-2
207 *
208 * and (S11, T11) corresponds to the complex conjugate eigenvalue
209 * pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
210 * that
211 *
212 * U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
213 * ( 0 s22 ) ( 0 t22 )
214 *
215 * where the generalized eigenvalues w = s11/t11 and
216 * conjg(w) = s22/t22.
217 *
218 * Then the reciprocal condition number DIF(i) is bounded by
219 *
220 * min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
221 *
222 * where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
223 * Z1 is the complex 2-by-2 matrix
224 *
225 * Z1 = [ s11 -s22 ]
226 * [ t11 -t22 ],
227 *
228 * This is done by computing (using real arithmetic) the
229 * roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
230 * where Z1**T denotes the transpose of Z1 and det(X) denotes
231 * the determinant of X.
232 *
233 * and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
234 * upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
235 *
236 * Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ]
237 * [ kron(T11**T, In-2) -kron(I2, T22) ]
238 *
239 * Note that if the default method for computing DIF is wanted (see
240 * DLATDF), then the parameter DIFDRI (see below) should be changed
241 * from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
242 * for more details.
243 *
244 * For each eigenvalue/vector specified by SELECT, DIF stores a
245 * Frobenius norm-based estimate of Difl.
246 *
247 * An approximate error bound for the i-th computed eigenvector VL(i) or
248 * VR(i) is given by
249 *
250 * EPS * norm(A, B) / DIF(i).
251 *
252 * See ref. [2-3] for more details and further references.
253 *
254 * Based on contributions by
255 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
256 * Umea University, S-901 87 Umea, Sweden.
257 *
258 * References
259 * ==========
260 *
261 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
262 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
263 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
264 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
265 *
266 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
267 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
268 * Estimation: Theory, Algorithms and Software,
269 * Report UMINF - 94.04, Department of Computing Science, Umea
270 * University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
271 * Note 87. To appear in Numerical Algorithms, 1996.
272 *
273 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
274 * for Solving the Generalized Sylvester Equation and Estimating the
275 * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
276 * Department of Computing Science, Umea University, S-901 87 Umea,
277 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
278 * Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
279 * No 1, 1996.
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284 INTEGER DIFDRI
285 PARAMETER ( DIFDRI = 3 )
286 DOUBLE PRECISION ZERO, ONE, TWO, FOUR
287 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
288 $ FOUR = 4.0D+0 )
289 * ..
290 * .. Local Scalars ..
291 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
292 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
293 DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
294 $ EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
295 $ TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
296 $ UHBVI
297 * ..
298 * .. Local Arrays ..
299 DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 )
300 * ..
301 * .. External Functions ..
302 LOGICAL LSAME
303 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
304 EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
305 * ..
306 * .. External Subroutines ..
307 EXTERNAL DGEMV, DLACPY, DLAG2, DTGEXC, DTGSYL, XERBLA
308 * ..
309 * .. Intrinsic Functions ..
310 INTRINSIC MAX, MIN, SQRT
311 * ..
312 * .. Executable Statements ..
313 *
314 * Decode and test the input parameters
315 *
316 WANTBH = LSAME( JOB, 'B' )
317 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
318 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
319 *
320 SOMCON = LSAME( HOWMNY, 'S' )
321 *
322 INFO = 0
323 LQUERY = ( LWORK.EQ.-1 )
324 *
325 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
326 INFO = -1
327 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
328 INFO = -2
329 ELSE IF( N.LT.0 ) THEN
330 INFO = -4
331 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
332 INFO = -6
333 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
334 INFO = -8
335 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
336 INFO = -10
337 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
338 INFO = -12
339 ELSE
340 *
341 * Set M to the number of eigenpairs for which condition numbers
342 * are required, and test MM.
343 *
344 IF( SOMCON ) THEN
345 M = 0
346 PAIR = .FALSE.
347 DO 10 K = 1, N
348 IF( PAIR ) THEN
349 PAIR = .FALSE.
350 ELSE
351 IF( K.LT.N ) THEN
352 IF( A( K+1, K ).EQ.ZERO ) THEN
353 IF( SELECT( K ) )
354 $ M = M + 1
355 ELSE
356 PAIR = .TRUE.
357 IF( SELECT( K ) .OR. SELECT( K+1 ) )
358 $ M = M + 2
359 END IF
360 ELSE
361 IF( SELECT( N ) )
362 $ M = M + 1
363 END IF
364 END IF
365 10 CONTINUE
366 ELSE
367 M = N
368 END IF
369 *
370 IF( N.EQ.0 ) THEN
371 LWMIN = 1
372 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
373 LWMIN = 2*N*( N + 2 ) + 16
374 ELSE
375 LWMIN = N
376 END IF
377 WORK( 1 ) = LWMIN
378 *
379 IF( MM.LT.M ) THEN
380 INFO = -15
381 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
382 INFO = -18
383 END IF
384 END IF
385 *
386 IF( INFO.NE.0 ) THEN
387 CALL XERBLA( 'DTGSNA', -INFO )
388 RETURN
389 ELSE IF( LQUERY ) THEN
390 RETURN
391 END IF
392 *
393 * Quick return if possible
394 *
395 IF( N.EQ.0 )
396 $ RETURN
397 *
398 * Get machine constants
399 *
400 EPS = DLAMCH( 'P' )
401 SMLNUM = DLAMCH( 'S' ) / EPS
402 KS = 0
403 PAIR = .FALSE.
404 *
405 DO 20 K = 1, N
406 *
407 * Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
408 *
409 IF( PAIR ) THEN
410 PAIR = .FALSE.
411 GO TO 20
412 ELSE
413 IF( K.LT.N )
414 $ PAIR = A( K+1, K ).NE.ZERO
415 END IF
416 *
417 * Determine whether condition numbers are required for the k-th
418 * eigenpair.
419 *
420 IF( SOMCON ) THEN
421 IF( PAIR ) THEN
422 IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
423 $ GO TO 20
424 ELSE
425 IF( .NOT.SELECT( K ) )
426 $ GO TO 20
427 END IF
428 END IF
429 *
430 KS = KS + 1
431 *
432 IF( WANTS ) THEN
433 *
434 * Compute the reciprocal condition number of the k-th
435 * eigenvalue.
436 *
437 IF( PAIR ) THEN
438 *
439 * Complex eigenvalue pair.
440 *
441 RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
442 $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
443 LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
444 $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
445 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
446 $ WORK, 1 )
447 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
448 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
449 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
450 $ ZERO, WORK, 1 )
451 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
452 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
453 UHAV = TMPRR + TMPII
454 UHAVI = TMPIR - TMPRI
455 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
456 $ WORK, 1 )
457 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
458 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
459 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
460 $ ZERO, WORK, 1 )
461 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
462 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
463 UHBV = TMPRR + TMPII
464 UHBVI = TMPIR - TMPRI
465 UHAV = DLAPY2( UHAV, UHAVI )
466 UHBV = DLAPY2( UHBV, UHBVI )
467 COND = DLAPY2( UHAV, UHBV )
468 S( KS ) = COND / ( RNRM*LNRM )
469 S( KS+1 ) = S( KS )
470 *
471 ELSE
472 *
473 * Real eigenvalue.
474 *
475 RNRM = DNRM2( N, VR( 1, KS ), 1 )
476 LNRM = DNRM2( N, VL( 1, KS ), 1 )
477 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
478 $ WORK, 1 )
479 UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
480 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
481 $ WORK, 1 )
482 UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
483 COND = DLAPY2( UHAV, UHBV )
484 IF( COND.EQ.ZERO ) THEN
485 S( KS ) = -ONE
486 ELSE
487 S( KS ) = COND / ( RNRM*LNRM )
488 END IF
489 END IF
490 END IF
491 *
492 IF( WANTDF ) THEN
493 IF( N.EQ.1 ) THEN
494 DIF( KS ) = DLAPY2( A( 1, 1 ), B( 1, 1 ) )
495 GO TO 20
496 END IF
497 *
498 * Estimate the reciprocal condition number of the k-th
499 * eigenvectors.
500 IF( PAIR ) THEN
501 *
502 * Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)).
503 * Compute the eigenvalue(s) at position K.
504 *
505 WORK( 1 ) = A( K, K )
506 WORK( 2 ) = A( K+1, K )
507 WORK( 3 ) = A( K, K+1 )
508 WORK( 4 ) = A( K+1, K+1 )
509 WORK( 5 ) = B( K, K )
510 WORK( 6 ) = B( K+1, K )
511 WORK( 7 ) = B( K, K+1 )
512 WORK( 8 ) = B( K+1, K+1 )
513 CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
514 $ DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
515 ALPRQT = ONE
516 C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
517 C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
518 ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 )
519 ROOT2 = C2 / ROOT1
520 ROOT1 = ROOT1 / TWO
521 COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
522 END IF
523 *
524 * Copy the matrix (A, B) to the array WORK and swap the
525 * diagonal block beginning at A(k,k) to the (1,1) position.
526 *
527 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
528 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
529 IFST = K
530 ILST = 1
531 *
532 CALL DTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
533 $ DUMMY, 1, DUMMY1, 1, IFST, ILST,
534 $ WORK( N*N*2+1 ), LWORK-2*N*N, IERR )
535 *
536 IF( IERR.GT.0 ) THEN
537 *
538 * Ill-conditioned problem - swap rejected.
539 *
540 DIF( KS ) = ZERO
541 ELSE
542 *
543 * Reordering successful, solve generalized Sylvester
544 * equation for R and L,
545 * A22 * R - L * A11 = A12
546 * B22 * R - L * B11 = B12,
547 * and compute estimate of Difl((A11,B11), (A22, B22)).
548 *
549 N1 = 1
550 IF( WORK( 2 ).NE.ZERO )
551 $ N1 = 2
552 N2 = N - N1
553 IF( N2.EQ.0 ) THEN
554 DIF( KS ) = COND
555 ELSE
556 I = N*N + 1
557 IZ = 2*N*N + 1
558 CALL DTGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ),
559 $ N, WORK, N, WORK( N1+1 ), N,
560 $ WORK( N*N1+N1+I ), N, WORK( I ), N,
561 $ WORK( N1+I ), N, SCALE, DIF( KS ),
562 $ WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
563 *
564 IF( PAIR )
565 $ DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
566 $ COND )
567 END IF
568 END IF
569 IF( PAIR )
570 $ DIF( KS+1 ) = DIF( KS )
571 END IF
572 IF( PAIR )
573 $ KS = KS + 1
574 *
575 20 CONTINUE
576 WORK( 1 ) = LWMIN
577 RETURN
578 *
579 * End of DTGSNA
580 *
581 END