LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dtgex2.f
Go to the documentation of this file.
1 *> \brief \b DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTGEX2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
22 * LDZ, J1, N1, N2, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * LOGICAL WANTQ, WANTZ
26 * INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30 * $ WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
40 *> of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
41 *> (A, B) by an orthogonal equivalence transformation.
42 *>
43 *> (A, B) must be in generalized real Schur canonical form (as returned
44 *> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
45 *> diagonal blocks. B is upper triangular.
46 *>
47 *> Optionally, the matrices Q and Z of generalized Schur vectors are
48 *> updated.
49 *>
50 *> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
51 *> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
52 *>
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] WANTQ
59 *> \verbatim
60 *> WANTQ is LOGICAL
61 *> .TRUE. : update the left transformation matrix Q;
62 *> .FALSE.: do not update Q.
63 *> \endverbatim
64 *>
65 *> \param[in] WANTZ
66 *> \verbatim
67 *> WANTZ is LOGICAL
68 *> .TRUE. : update the right transformation matrix Z;
69 *> .FALSE.: do not update Z.
70 *> \endverbatim
71 *>
72 *> \param[in] N
73 *> \verbatim
74 *> N is INTEGER
75 *> The order of the matrices A and B. N >= 0.
76 *> \endverbatim
77 *>
78 *> \param[in,out] A
79 *> \verbatim
80 *> A is DOUBLE PRECISION array, dimensions (LDA,N)
81 *> On entry, the matrix A in the pair (A, B).
82 *> On exit, the updated matrix A.
83 *> \endverbatim
84 *>
85 *> \param[in] LDA
86 *> \verbatim
87 *> LDA is INTEGER
88 *> The leading dimension of the array A. LDA >= max(1,N).
89 *> \endverbatim
90 *>
91 *> \param[in,out] B
92 *> \verbatim
93 *> B is DOUBLE PRECISION array, dimensions (LDB,N)
94 *> On entry, the matrix B in the pair (A, B).
95 *> On exit, the updated matrix B.
96 *> \endverbatim
97 *>
98 *> \param[in] LDB
99 *> \verbatim
100 *> LDB is INTEGER
101 *> The leading dimension of the array B. LDB >= max(1,N).
102 *> \endverbatim
103 *>
104 *> \param[in,out] Q
105 *> \verbatim
106 *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
107 *> On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
108 *> On exit, the updated matrix Q.
109 *> Not referenced if WANTQ = .FALSE..
110 *> \endverbatim
111 *>
112 *> \param[in] LDQ
113 *> \verbatim
114 *> LDQ is INTEGER
115 *> The leading dimension of the array Q. LDQ >= 1.
116 *> If WANTQ = .TRUE., LDQ >= N.
117 *> \endverbatim
118 *>
119 *> \param[in,out] Z
120 *> \verbatim
121 *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
122 *> On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
123 *> On exit, the updated matrix Z.
124 *> Not referenced if WANTZ = .FALSE..
125 *> \endverbatim
126 *>
127 *> \param[in] LDZ
128 *> \verbatim
129 *> LDZ is INTEGER
130 *> The leading dimension of the array Z. LDZ >= 1.
131 *> If WANTZ = .TRUE., LDZ >= N.
132 *> \endverbatim
133 *>
134 *> \param[in] J1
135 *> \verbatim
136 *> J1 is INTEGER
137 *> The index to the first block (A11, B11). 1 <= J1 <= N.
138 *> \endverbatim
139 *>
140 *> \param[in] N1
141 *> \verbatim
142 *> N1 is INTEGER
143 *> The order of the first block (A11, B11). N1 = 0, 1 or 2.
144 *> \endverbatim
145 *>
146 *> \param[in] N2
147 *> \verbatim
148 *> N2 is INTEGER
149 *> The order of the second block (A22, B22). N2 = 0, 1 or 2.
150 *> \endverbatim
151 *>
152 *> \param[out] WORK
153 *> \verbatim
154 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
155 *> \endverbatim
156 *>
157 *> \param[in] LWORK
158 *> \verbatim
159 *> LWORK is INTEGER
160 *> The dimension of the array WORK.
161 *> LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
162 *> \endverbatim
163 *>
164 *> \param[out] INFO
165 *> \verbatim
166 *> INFO is INTEGER
167 *> =0: Successful exit
168 *> >0: If INFO = 1, the transformed matrix (A, B) would be
169 *> too far from generalized Schur form; the blocks are
170 *> not swapped and (A, B) and (Q, Z) are unchanged.
171 *> The problem of swapping is too ill-conditioned.
172 *> <0: If INFO = -16: LWORK is too small. Appropriate value
173 *> for LWORK is returned in WORK(1).
174 *> \endverbatim
175 *
176 * Authors:
177 * ========
178 *
179 *> \author Univ. of Tennessee
180 *> \author Univ. of California Berkeley
181 *> \author Univ. of Colorado Denver
182 *> \author NAG Ltd.
183 *
184 *> \date September 2012
185 *
186 *> \ingroup doubleGEauxiliary
187 *
188 *> \par Further Details:
189 * =====================
190 *>
191 *> In the current code both weak and strong stability tests are
192 *> performed. The user can omit the strong stability test by changing
193 *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
194 *> details.
195 *
196 *> \par Contributors:
197 * ==================
198 *>
199 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
200 *> Umea University, S-901 87 Umea, Sweden.
201 *
202 *> \par References:
203 * ================
204 *>
205 *> \verbatim
206 *>
207 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
208 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
209 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
210 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
211 *>
212 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
213 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
214 *> Estimation: Theory, Algorithms and Software,
215 *> Report UMINF - 94.04, Department of Computing Science, Umea
216 *> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
217 *> Note 87. To appear in Numerical Algorithms, 1996.
218 *> \endverbatim
219 *>
220 * =====================================================================
221  SUBROUTINE dtgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
222  $ ldz, j1, n1, n2, work, lwork, info )
223 *
224 * -- LAPACK auxiliary routine (version 3.4.2) --
225 * -- LAPACK is a software package provided by Univ. of Tennessee, --
226 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * September 2012
228 *
229 * .. Scalar Arguments ..
230  LOGICAL wantq, wantz
231  INTEGER info, j1, lda, ldb, ldq, ldz, lwork, n, n1, n2
232 * ..
233 * .. Array Arguments ..
234  DOUBLE PRECISION a( lda, * ), b( ldb, * ), q( ldq, * ),
235  $ work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 * Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
240 * loops. Sven Hammarling, 1/5/02.
241 *
242 * .. Parameters ..
243  DOUBLE PRECISION zero, one
244  parameter( zero = 0.0d+0, one = 1.0d+0 )
245  DOUBLE PRECISION twenty
246  parameter( twenty = 2.0d+01 )
247  INTEGER ldst
248  parameter( ldst = 4 )
249  LOGICAL wands
250  parameter( wands = .true. )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL dtrong, weak
254  INTEGER i, idum, linfo, m
255  DOUBLE PRECISION bqra21, brqa21, ddum, dnorm, dscale, dsum, eps,
256  $ f, g, sa, sb, scale, smlnum, ss, thresh, ws
257 * ..
258 * .. Local Arrays ..
259  INTEGER iwork( ldst )
260  DOUBLE PRECISION ai( 2 ), ar( 2 ), be( 2 ), ir( ldst, ldst ),
261  $ ircop( ldst, ldst ), li( ldst, ldst ),
262  $ licop( ldst, ldst ), s( ldst, ldst ),
263  $ scpy( ldst, ldst ), t( ldst, ldst ),
264  $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
265 * ..
266 * .. External Functions ..
267  DOUBLE PRECISION dlamch
268  EXTERNAL dlamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL dgemm, dgeqr2, dgerq2, dlacpy, dlagv2, dlartg,
273  $ drot, dscal, dtgsy2
274 * ..
275 * .. Intrinsic Functions ..
276  INTRINSIC abs, max, sqrt
277 * ..
278 * .. Executable Statements ..
279 *
280  info = 0
281 *
282 * Quick return if possible
283 *
284  IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
285  $ return
286  IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
287  $ return
288  m = n1 + n2
289  IF( lwork.LT.max( 1, n*m, m*m*2 ) ) THEN
290  info = -16
291  work( 1 ) = max( 1, n*m, m*m*2 )
292  return
293  END IF
294 *
295  weak = .false.
296  dtrong = .false.
297 *
298 * Make a local copy of selected block
299 *
300  CALL dlaset( 'Full', ldst, ldst, zero, zero, li, ldst )
301  CALL dlaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
302  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
303  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
304 *
305 * Compute threshold for testing acceptance of swapping.
306 *
307  eps = dlamch( 'P' )
308  smlnum = dlamch( 'S' ) / eps
309  dscale = zero
310  dsum = one
311  CALL dlacpy( 'Full', m, m, s, ldst, work, m )
312  CALL dlassq( m*m, work, 1, dscale, dsum )
313  CALL dlacpy( 'Full', m, m, t, ldst, work, m )
314  CALL dlassq( m*m, work, 1, dscale, dsum )
315  dnorm = dscale*sqrt( dsum )
316 *
317 * THRES has been changed from
318 * THRESH = MAX( TEN*EPS*SA, SMLNUM )
319 * to
320 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
321 * on 04/01/10.
322 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
323 * Jim Demmel and Guillaume Revy. See forum post 1783.
324 *
325  thresh = max( twenty*eps*dnorm, smlnum )
326 *
327  IF( m.EQ.2 ) THEN
328 *
329 * CASE 1: Swap 1-by-1 and 1-by-1 blocks.
330 *
331 * Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
332 * using Givens rotations and perform the swap tentatively.
333 *
334  f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
335  g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
336  sb = abs( t( 2, 2 ) )
337  sa = abs( s( 2, 2 ) )
338  CALL dlartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
339  ir( 2, 1 ) = -ir( 1, 2 )
340  ir( 2, 2 ) = ir( 1, 1 )
341  CALL drot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
342  $ ir( 2, 1 ) )
343  CALL drot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
344  $ ir( 2, 1 ) )
345  IF( sa.GE.sb ) THEN
346  CALL dlartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
347  $ ddum )
348  ELSE
349  CALL dlartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  END IF
352  CALL drot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
353  $ li( 2, 1 ) )
354  CALL drot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
355  $ li( 2, 1 ) )
356  li( 2, 2 ) = li( 1, 1 )
357  li( 1, 2 ) = -li( 2, 1 )
358 *
359 * Weak stability test:
360 * |S21| + |T21| <= O(EPS * F-norm((S, T)))
361 *
362  ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
363  weak = ws.LE.thresh
364  IF( .NOT.weak )
365  $ go to 70
366 *
367  IF( wands ) THEN
368 *
369 * Strong stability test:
370 * F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))
371 *
372  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
373  $ m )
374  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
375  $ work, m )
376  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
377  $ work( m*m+1 ), m )
378  dscale = zero
379  dsum = one
380  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
381 *
382  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
383  $ m )
384  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
385  $ work, m )
386  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
387  $ work( m*m+1 ), m )
388  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389  ss = dscale*sqrt( dsum )
390  dtrong = ss.LE.thresh
391  IF( .NOT.dtrong )
392  $ go to 70
393  END IF
394 *
395 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
396 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
397 *
398  CALL drot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
399  $ ir( 2, 1 ) )
400  CALL drot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
401  $ ir( 2, 1 ) )
402  CALL drot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403  $ li( 1, 1 ), li( 2, 1 ) )
404  CALL drot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
405  $ li( 1, 1 ), li( 2, 1 ) )
406 *
407 * Set N1-by-N2 (2,1) - blocks to ZERO.
408 *
409  a( j1+1, j1 ) = zero
410  b( j1+1, j1 ) = zero
411 *
412 * Accumulate transformations into Q and Z if requested.
413 *
414  IF( wantz )
415  $ CALL drot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
416  $ ir( 2, 1 ) )
417  IF( wantq )
418  $ CALL drot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
419  $ li( 2, 1 ) )
420 *
421 * Exit with INFO = 0 if swap was successfully performed.
422 *
423  return
424 *
425  ELSE
426 *
427 * CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
428 * and 2-by-2 blocks.
429 *
430 * Solve the generalized Sylvester equation
431 * S11 * R - L * S22 = SCALE * S12
432 * T11 * R - L * T22 = SCALE * T12
433 * for R and L. Solutions in LI and IR.
434 *
435  CALL dlacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436  CALL dlacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
437  $ ir( n2+1, n1+1 ), ldst )
438  CALL dtgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
439  $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
440  $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
441  $ linfo )
442 *
443 * Compute orthogonal matrix QL:
444 *
445 * QL**T * LI = [ TL ]
446 * [ 0 ]
447 * where
448 * LI = [ -L ]
449 * [ SCALE * identity(N2) ]
450 *
451  DO 10 i = 1, n2
452  CALL dscal( n1, -one, li( 1, i ), 1 )
453  li( n1+i, i ) = scale
454  10 continue
455  CALL dgeqr2( m, n2, li, ldst, taul, work, linfo )
456  IF( linfo.NE.0 )
457  $ go to 70
458  CALL dorg2r( m, m, n2, li, ldst, taul, work, linfo )
459  IF( linfo.NE.0 )
460  $ go to 70
461 *
462 * Compute orthogonal matrix RQ:
463 *
464 * IR * RQ**T = [ 0 TR],
465 *
466 * where IR = [ SCALE * identity(N1), R ]
467 *
468  DO 20 i = 1, n1
469  ir( n2+i, i ) = scale
470  20 continue
471  CALL dgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
472  IF( linfo.NE.0 )
473  $ go to 70
474  CALL dorgr2( m, m, n1, ir, ldst, taur, work, linfo )
475  IF( linfo.NE.0 )
476  $ go to 70
477 *
478 * Perform the swapping tentatively:
479 *
480  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
481  $ work, m )
482  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
483  $ ldst )
484  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
485  $ work, m )
486  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
487  $ ldst )
488  CALL dlacpy( 'F', m, m, s, ldst, scpy, ldst )
489  CALL dlacpy( 'F', m, m, t, ldst, tcpy, ldst )
490  CALL dlacpy( 'F', m, m, ir, ldst, ircop, ldst )
491  CALL dlacpy( 'F', m, m, li, ldst, licop, ldst )
492 *
493 * Triangularize the B-part by an RQ factorization.
494 * Apply transformation (from left) to A-part, giving S.
495 *
496  CALL dgerq2( m, m, t, ldst, taur, work, linfo )
497  IF( linfo.NE.0 )
498  $ go to 70
499  CALL dormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
500  $ linfo )
501  IF( linfo.NE.0 )
502  $ go to 70
503  CALL dormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
504  $ linfo )
505  IF( linfo.NE.0 )
506  $ go to 70
507 *
508 * Compute F-norm(S21) in BRQA21. (T21 is 0.)
509 *
510  dscale = zero
511  dsum = one
512  DO 30 i = 1, n2
513  CALL dlassq( n1, s( n2+1, i ), 1, dscale, dsum )
514  30 continue
515  brqa21 = dscale*sqrt( dsum )
516 *
517 * Triangularize the B-part by a QR factorization.
518 * Apply transformation (from right) to A-part, giving S.
519 *
520  CALL dgeqr2( m, m, tcpy, ldst, taul, work, linfo )
521  IF( linfo.NE.0 )
522  $ go to 70
523  CALL dorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
524  $ work, info )
525  CALL dorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
526  $ work, info )
527  IF( linfo.NE.0 )
528  $ go to 70
529 *
530 * Compute F-norm(S21) in BQRA21. (T21 is 0.)
531 *
532  dscale = zero
533  dsum = one
534  DO 40 i = 1, n2
535  CALL dlassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
536  40 continue
537  bqra21 = dscale*sqrt( dsum )
538 *
539 * Decide which method to use.
540 * Weak stability test:
541 * F-norm(S21) <= O(EPS * F-norm((S, T)))
542 *
543  IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresh ) THEN
544  CALL dlacpy( 'F', m, m, scpy, ldst, s, ldst )
545  CALL dlacpy( 'F', m, m, tcpy, ldst, t, ldst )
546  CALL dlacpy( 'F', m, m, ircop, ldst, ir, ldst )
547  CALL dlacpy( 'F', m, m, licop, ldst, li, ldst )
548  ELSE IF( brqa21.GE.thresh ) THEN
549  go to 70
550  END IF
551 *
552 * Set lower triangle of B-part to zero
553 *
554  CALL dlaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
555 *
556  IF( wands ) THEN
557 *
558 * Strong stability test:
559 * F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
560 *
561  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
562  $ m )
563  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
564  $ work, m )
565  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
566  $ work( m*m+1 ), m )
567  dscale = zero
568  dsum = one
569  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
570 *
571  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
572  $ m )
573  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
574  $ work, m )
575  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
576  $ work( m*m+1 ), m )
577  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578  ss = dscale*sqrt( dsum )
579  dtrong = ( ss.LE.thresh )
580  IF( .NOT.dtrong )
581  $ go to 70
582 *
583  END IF
584 *
585 * If the swap is accepted ("weakly" and "strongly"), apply the
586 * transformations and set N1-by-N2 (2,1)-block to zero.
587 *
588  CALL dlaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
589 *
590 * copy back M-by-M diagonal block starting at index J1 of (A, B)
591 *
592  CALL dlacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
593  CALL dlacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
594  CALL dlaset( 'Full', ldst, ldst, zero, zero, t, ldst )
595 *
596 * Standardize existing 2-by-2 blocks.
597 *
598  DO 50 i = 1, m*m
599  work(i) = zero
600  50 continue
601  work( 1 ) = one
602  t( 1, 1 ) = one
603  idum = lwork - m*m - 2
604  IF( n2.GT.1 ) THEN
605  CALL dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
606  $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
607  work( m+1 ) = -work( 2 )
608  work( m+2 ) = work( 1 )
609  t( n2, n2 ) = t( 1, 1 )
610  t( 1, 2 ) = -t( 2, 1 )
611  END IF
612  work( m*m ) = one
613  t( m, m ) = one
614 *
615  IF( n1.GT.1 ) THEN
616  CALL dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
617  $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
618  $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
619  $ t( m, m-1 ) )
620  work( m*m ) = work( n2*m+n2+1 )
621  work( m*m-1 ) = -work( n2*m+n2+2 )
622  t( m, m ) = t( n2+1, n2+1 )
623  t( m-1, m ) = -t( m, m-1 )
624  END IF
625  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
626  $ lda, zero, work( m*m+1 ), n2 )
627  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
628  $ lda )
629  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
630  $ ldb, zero, work( m*m+1 ), n2 )
631  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
632  $ ldb )
633  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
634  $ work( m*m+1 ), m )
635  CALL dlacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
636  CALL dgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
637  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
638  CALL dlacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
639  CALL dgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
640  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
641  CALL dlacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
642  CALL dgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
643  $ work, m )
644  CALL dlacpy( 'Full', m, m, work, m, ir, ldst )
645 *
646 * Accumulate transformations into Q and Z if requested.
647 *
648  IF( wantq ) THEN
649  CALL dgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
650  $ ldst, zero, work, n )
651  CALL dlacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
652 *
653  END IF
654 *
655  IF( wantz ) THEN
656  CALL dgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
657  $ ldst, zero, work, n )
658  CALL dlacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
659 *
660  END IF
661 *
662 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
663 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
664 *
665  i = j1 + m
666  IF( i.LE.n ) THEN
667  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
668  $ a( j1, i ), lda, zero, work, m )
669  CALL dlacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
670  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
671  $ b( j1, i ), lda, zero, work, m )
672  CALL dlacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
673  END IF
674  i = j1 - 1
675  IF( i.GT.0 ) THEN
676  CALL dgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
677  $ ldst, zero, work, i )
678  CALL dlacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
679  CALL dgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
680  $ ldst, zero, work, i )
681  CALL dlacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
682  END IF
683 *
684 * Exit with INFO = 0 if swap was successfully performed.
685 *
686  return
687 *
688  END IF
689 *
690 * Exit with INFO = 1 if swap was rejected.
691 *
692  70 continue
693 *
694  info = 1
695  return
696 *
697 * End of DTGEX2
698 *
699  END