LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
stgex2.f
Go to the documentation of this file.
1 *> \brief \b STGEX2 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 STGEX2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgex2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgex2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgex2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STGEX2( 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 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30 * $ WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> STGEX2 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 SGGES), 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 REAL arrays, 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 REAL arrays, 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 REAL array, dimension (LDZ,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 REAL 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 REAL 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( 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 realGEauxiliary
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 stgex2( 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  REAL A( lda, * ), B( ldb, * ), Q( ldq, * ),
235  $ work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 * Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
240 * loops. Sven Hammarling, 1/5/02.
241 *
242 * .. Parameters ..
243  REAL ZERO, ONE
244  parameter( zero = 0.0e+0, one = 1.0e+0 )
245  REAL TWENTY
246  parameter( twenty = 2.0e+01 )
247  INTEGER LDST
248  parameter( ldst = 4 )
249  LOGICAL WANDS
250  parameter( wands = .true. )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL STRONG, WEAK
254  INTEGER I, IDUM, LINFO, M
255  REAL 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  REAL 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  REAL SLAMCH
268  EXTERNAL slamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL sgemm, sgeqr2, sgerq2, slacpy, slagv2, slartg,
273  $ srot, sscal, stgsy2
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( n*m, m*m*2 ) ) THEN
290  info = -16
291  work( 1 ) = max( n*m, m*m*2 )
292  RETURN
293  END IF
294 *
295  weak = .false.
296  strong = .false.
297 *
298 * Make a local copy of selected block
299 *
300  CALL slaset( 'Full', ldst, ldst, zero, zero, li, ldst )
301  CALL slaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
302  CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
303  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
304 *
305 * Compute threshold for testing acceptance of swapping.
306 *
307  eps = slamch( 'P' )
308  smlnum = slamch( 'S' ) / eps
309  dscale = zero
310  dsum = one
311  CALL slacpy( 'Full', m, m, s, ldst, work, m )
312  CALL slassq( m*m, work, 1, dscale, dsum )
313  CALL slacpy( 'Full', m, m, t, ldst, work, m )
314  CALL slassq( 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 slartg( 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 srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
342  $ ir( 2, 1 ) )
343  CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
344  $ ir( 2, 1 ) )
345  IF( sa.GE.sb ) THEN
346  CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
347  $ ddum )
348  ELSE
349  CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  END IF
352  CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
353  $ li( 2, 1 ) )
354  CALL srot( 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 slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
373  $ m )
374  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
375  $ work, m )
376  CALL sgemm( '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 slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
381 *
382  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
383  $ m )
384  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
385  $ work, m )
386  CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
387  $ work( m*m+1 ), m )
388  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389  ss = dscale*sqrt( dsum )
390  strong = ss.LE.thresh
391  IF( .NOT.strong )
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 srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
399  $ ir( 2, 1 ) )
400  CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
401  $ ir( 2, 1 ) )
402  CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403  $ li( 1, 1 ), li( 2, 1 ) )
404  CALL srot( 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 srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
416  $ ir( 2, 1 ) )
417  IF( wantq )
418  $ CALL srot( 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 slacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436  CALL slacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
437  $ ir( n2+1, n1+1 ), ldst )
438  CALL stgsy2( '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 sscal( n1, -one, li( 1, i ), 1 )
453  li( n1+i, i ) = scale
454  10 CONTINUE
455  CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
456  IF( linfo.NE.0 )
457  $ GO TO 70
458  CALL sorg2r( 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 sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
472  IF( linfo.NE.0 )
473  $ GO TO 70
474  CALL sorgr2( 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 sgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
481  $ work, m )
482  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
483  $ ldst )
484  CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
485  $ work, m )
486  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
487  $ ldst )
488  CALL slacpy( 'F', m, m, s, ldst, scpy, ldst )
489  CALL slacpy( 'F', m, m, t, ldst, tcpy, ldst )
490  CALL slacpy( 'F', m, m, ir, ldst, ircop, ldst )
491  CALL slacpy( '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 sgerq2( m, m, t, ldst, taur, work, linfo )
497  IF( linfo.NE.0 )
498  $ GO TO 70
499  CALL sormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
500  $ linfo )
501  IF( linfo.NE.0 )
502  $ GO TO 70
503  CALL sormr2( '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 slassq( 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 sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
521  IF( linfo.NE.0 )
522  $ GO TO 70
523  CALL sorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
524  $ work, info )
525  CALL sorm2r( '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 slassq( 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 slacpy( 'F', m, m, scpy, ldst, s, ldst )
545  CALL slacpy( 'F', m, m, tcpy, ldst, t, ldst )
546  CALL slacpy( 'F', m, m, ircop, ldst, ir, ldst )
547  CALL slacpy( '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 slaset( '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 slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
562  $ m )
563  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
564  $ work, m )
565  CALL sgemm( '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 slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
570 *
571  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
572  $ m )
573  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
574  $ work, m )
575  CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
576  $ work( m*m+1 ), m )
577  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578  ss = dscale*sqrt( dsum )
579  strong = ( ss.LE.thresh )
580  IF( .NOT.strong )
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 slaset( '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 slacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
593  CALL slacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
594  CALL slaset( 'Full', ldst, ldst, zero, zero, t, ldst )
595 *
596 * Standardize existing 2-by-2 blocks.
597 *
598  CALL slaset( '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 slagv2( 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 slagv2( 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 sgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624  $ lda, zero, work( m*m+1 ), n2 )
625  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
626  $ lda )
627  CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628  $ ldb, zero, work( m*m+1 ), n2 )
629  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
630  $ ldb )
631  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
632  $ work( m*m+1 ), m )
633  CALL slacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
634  CALL sgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636  CALL slacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637  CALL sgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639  CALL slacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640  CALL sgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
641  $ work, m )
642  CALL slacpy( 'Full', m, m, work, m, ir, ldst )
643 *
644 * Accumulate transformations into Q and Z if requested.
645 *
646  IF( wantq ) THEN
647  CALL sgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
648  $ ldst, zero, work, n )
649  CALL slacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
650 *
651  END IF
652 *
653  IF( wantz ) THEN
654  CALL sgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655  $ ldst, zero, work, n )
656  CALL slacpy( '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 sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
666  $ a( j1, i ), lda, zero, work, m )
667  CALL slacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
668  CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
669  $ b( j1, i ), ldb, zero, work, m )
670  CALL slacpy( '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 sgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
675  $ ldst, zero, work, i )
676  CALL slacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
677  CALL sgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678  $ ldst, zero, work, i )
679  CALL slacpy( '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 STGEX2
696 *
697  END
subroutine sorgr2(M, N, K, A, LDA, TAU, WORK, INFO)
SORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
Definition: sorgr2.f:116
subroutine sormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition: sormr2.f:161
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine stgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: stgsy2.f:276
subroutine slagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
Definition: slagv2.f:159
subroutine sorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
SORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition: sorg2r.f:116
subroutine sgerq2(M, N, A, LDA, TAU, WORK, INFO)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgerq2.f:125
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:161
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f:105
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgeqr2.f:123
subroutine stgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
STGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition: stgex2.f:223
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112