LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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 November 2015
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.6.0) --
225 * -- LAPACK is a software package provided by Univ. of Tennessee, --
226 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * November 2015
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  CALL dlaset( 'Full', m, m, zero, zero, work, m )
599  work( 1 ) = one
600  t( 1, 1 ) = one
601  idum = lwork - m*m - 2
602  IF( n2.GT.1 ) THEN
603  CALL dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
604  $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
605  work( m+1 ) = -work( 2 )
606  work( m+2 ) = work( 1 )
607  t( n2, n2 ) = t( 1, 1 )
608  t( 1, 2 ) = -t( 2, 1 )
609  END IF
610  work( m*m ) = one
611  t( m, m ) = one
612 *
613  IF( n1.GT.1 ) THEN
614  CALL dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
615  $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
616  $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
617  $ t( m, m-1 ) )
618  work( m*m ) = work( n2*m+n2+1 )
619  work( m*m-1 ) = -work( n2*m+n2+2 )
620  t( m, m ) = t( n2+1, n2+1 )
621  t( m-1, m ) = -t( m, m-1 )
622  END IF
623  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624  $ lda, zero, work( m*m+1 ), n2 )
625  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
626  $ lda )
627  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628  $ ldb, zero, work( m*m+1 ), n2 )
629  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
630  $ ldb )
631  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
632  $ work( m*m+1 ), m )
633  CALL dlacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
634  CALL dgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636  CALL dlacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637  CALL dgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639  CALL dlacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640  CALL dgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
641  $ work, m )
642  CALL dlacpy( 'Full', m, m, work, m, ir, ldst )
643 *
644 * Accumulate transformations into Q and Z if requested.
645 *
646  IF( wantq ) THEN
647  CALL dgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
648  $ ldst, zero, work, n )
649  CALL dlacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
650 *
651  END IF
652 *
653  IF( wantz ) THEN
654  CALL dgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655  $ ldst, zero, work, n )
656  CALL dlacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
657 *
658  END IF
659 *
660 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
661 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
662 *
663  i = j1 + m
664  IF( i.LE.n ) THEN
665  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
666  $ a( j1, i ), lda, zero, work, m )
667  CALL dlacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
668  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
669  $ b( j1, i ), ldb, zero, work, m )
670  CALL dlacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
671  END IF
672  i = j1 - 1
673  IF( i.GT.0 ) THEN
674  CALL dgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
675  $ ldst, zero, work, i )
676  CALL dlacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
677  CALL dgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678  $ ldst, zero, work, i )
679  CALL dlacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
680  END IF
681 *
682 * Exit with INFO = 0 if swap was successfully performed.
683 *
684  RETURN
685 *
686  END IF
687 *
688 * Exit with INFO = 1 if swap was rejected.
689 *
690  70 CONTINUE
691 *
692  info = 1
693  RETURN
694 *
695 * End of DTGEX2
696 *
697  END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
DORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition: dorg2r.f:116
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: dtgsy2.f:276
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgeqr2.f:123
subroutine dlagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
Definition: dlagv2.f:159
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:161
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgerq2.f:125
subroutine dormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition: dormr2.f:161
subroutine dtgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition: dtgex2.f:223
subroutine dorgr2(M, N, K, A, LDA, TAU, WORK, INFO)
DORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
Definition: dorgr2.f:116