LAPACK  3.10.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 array, dimension (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 array, dimension (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 (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 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 *> \ingroup realGEauxiliary
185 *
186 *> \par Further Details:
187 * =====================
188 *>
189 *> In the current code both weak and strong stability tests are
190 *> performed. The user can omit the strong stability test by changing
191 *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
192 *> details.
193 *
194 *> \par Contributors:
195 * ==================
196 *>
197 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
198 *> Umea University, S-901 87 Umea, Sweden.
199 *
200 *> \par References:
201 * ================
202 *>
203 *> \verbatim
204 *>
205 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
206 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
207 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
208 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
209 *>
210 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
211 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
212 *> Estimation: Theory, Algorithms and Software,
213 *> Report UMINF - 94.04, Department of Computing Science, Umea
214 *> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
215 *> Note 87. To appear in Numerical Algorithms, 1996.
216 *> \endverbatim
217 *>
218 * =====================================================================
219  SUBROUTINE stgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
220  $ LDZ, J1, N1, N2, WORK, LWORK, INFO )
221 *
222 * -- LAPACK auxiliary routine --
223 * -- LAPACK is a software package provided by Univ. of Tennessee, --
224 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225 *
226 * .. Scalar Arguments ..
227  LOGICAL WANTQ, WANTZ
228  INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
229 * ..
230 * .. Array Arguments ..
231  REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
232  $ work( * ), z( ldz, * )
233 * ..
234 *
235 * =====================================================================
236 * Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
237 * loops. Sven Hammarling, 1/5/02.
238 *
239 * .. Parameters ..
240  REAL ZERO, ONE
241  parameter( zero = 0.0e+0, one = 1.0e+0 )
242  REAL TWENTY
243  parameter( twenty = 2.0e+01 )
244  INTEGER LDST
245  parameter( ldst = 4 )
246  LOGICAL WANDS
247  parameter( wands = .true. )
248 * ..
249 * .. Local Scalars ..
250  LOGICAL STRONG, WEAK
251  INTEGER I, IDUM, LINFO, M
252  REAL BQRA21, BRQA21, DDUM, DNORMA, DNORMB,
253  $ dscale,
254  $ dsum, eps, f, g, sa, sb, scale, smlnum,
255  $ thresha, threshb
256 * ..
257 * .. Local Arrays ..
258  INTEGER IWORK( LDST )
259  REAL AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
260  $ ircop( ldst, ldst ), li( ldst, ldst ),
261  $ licop( ldst, ldst ), s( ldst, ldst ),
262  $ scpy( ldst, ldst ), t( ldst, ldst ),
263  $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
264 * ..
265 * .. External Functions ..
266  REAL SLAMCH
267  EXTERNAL slamch
268 * ..
269 * .. External Subroutines ..
270  EXTERNAL sgemm, sgeqr2, sgerq2, slacpy, slagv2, slartg,
272  $ srot, sscal, stgsy2
273 * ..
274 * .. Intrinsic Functions ..
275  INTRINSIC abs, max, sqrt
276 * ..
277 * .. Executable Statements ..
278 *
279  info = 0
280 *
281 * Quick return if possible
282 *
283  IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
284  $ RETURN
285  IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
286  $ RETURN
287  m = n1 + n2
288  IF( lwork.LT.max( n*m, m*m*2 ) ) THEN
289  info = -16
290  work( 1 ) = max( n*m, m*m*2 )
291  RETURN
292  END IF
293 *
294  weak = .false.
295  strong = .false.
296 *
297 * Make a local copy of selected block
298 *
299  CALL slaset( 'Full', ldst, ldst, zero, zero, li, ldst )
300  CALL slaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
301  CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
302  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
303 *
304 * Compute threshold for testing acceptance of swapping.
305 *
306  eps = slamch( 'P' )
307  smlnum = slamch( 'S' ) / eps
308  dscale = zero
309  dsum = one
310  CALL slacpy( 'Full', m, m, s, ldst, work, m )
311  CALL slassq( m*m, work, 1, dscale, dsum )
312  dnorma = dscale*sqrt( dsum )
313  dscale = zero
314  dsum = one
315  CALL slacpy( 'Full', m, m, t, ldst, work, m )
316  CALL slassq( m*m, work, 1, dscale, dsum )
317  dnormb = dscale*sqrt( dsum )
318 *
319 * THRES has been changed from
320 * THRESH = MAX( TEN*EPS*SA, SMLNUM )
321 * to
322 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
323 * on 04/01/10.
324 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
325 * Jim Demmel and Guillaume Revy. See forum post 1783.
326 *
327  thresha = max( twenty*eps*dnorma, smlnum )
328  threshb = max( twenty*eps*dnormb, smlnum )
329 *
330  IF( m.EQ.2 ) THEN
331 *
332 * CASE 1: Swap 1-by-1 and 1-by-1 blocks.
333 *
334 * Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
335 * using Givens rotations and perform the swap tentatively.
336 *
337  f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
338  g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
339  sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
340  sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
341  CALL slartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
342  ir( 2, 1 ) = -ir( 1, 2 )
343  ir( 2, 2 ) = ir( 1, 1 )
344  CALL srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
345  $ ir( 2, 1 ) )
346  CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
347  $ ir( 2, 1 ) )
348  IF( sa.GE.sb ) THEN
349  CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  ELSE
352  CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
353  $ ddum )
354  END IF
355  CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
356  $ li( 2, 1 ) )
357  CALL srot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
358  $ li( 2, 1 ) )
359  li( 2, 2 ) = li( 1, 1 )
360  li( 1, 2 ) = -li( 2, 1 )
361 *
362 * Weak stability test: |S21| <= O(EPS F-norm((A)))
363 * and |T21| <= O(EPS F-norm((B)))
364 *
365  weak = abs( s( 2, 1 ) ) .LE. thresha .AND.
366  $ abs( t( 2, 1 ) ) .LE. threshb
367  IF( .NOT.weak )
368  $ GO TO 70
369 *
370  IF( wands ) THEN
371 *
372 * Strong stability test:
373 * F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
374 * and
375 * F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
376 *
377  CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
378  $ m )
379  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
380  $ work, m )
381  CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
382  $ work( m*m+1 ), m )
383  dscale = zero
384  dsum = one
385  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
386  sa = dscale*sqrt( dsum )
387 *
388  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
389  $ m )
390  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
391  $ work, m )
392  CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
393  $ work( m*m+1 ), m )
394  dscale = zero
395  dsum = one
396  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
397  sb = dscale*sqrt( dsum )
398  strong = sa.LE.thresha .AND. sb.LE.threshb
399  IF( .NOT.strong )
400  $ GO TO 70
401  END IF
402 *
403 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
404 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
405 *
406  CALL srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
407  $ ir( 2, 1 ) )
408  CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
409  $ ir( 2, 1 ) )
410  CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
411  $ li( 1, 1 ), li( 2, 1 ) )
412  CALL srot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
413  $ li( 1, 1 ), li( 2, 1 ) )
414 *
415 * Set N1-by-N2 (2,1) - blocks to ZERO.
416 *
417  a( j1+1, j1 ) = zero
418  b( j1+1, j1 ) = zero
419 *
420 * Accumulate transformations into Q and Z if requested.
421 *
422  IF( wantz )
423  $ CALL srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
424  $ ir( 2, 1 ) )
425  IF( wantq )
426  $ CALL srot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
427  $ li( 2, 1 ) )
428 *
429 * Exit with INFO = 0 if swap was successfully performed.
430 *
431  RETURN
432 *
433  ELSE
434 *
435 * CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
436 * and 2-by-2 blocks.
437 *
438 * Solve the generalized Sylvester equation
439 * S11 * R - L * S22 = SCALE * S12
440 * T11 * R - L * T22 = SCALE * T12
441 * for R and L. Solutions in LI and IR.
442 *
443  CALL slacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
444  CALL slacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
445  $ ir( n2+1, n1+1 ), ldst )
446  CALL stgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
447  $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
448  $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
449  $ linfo )
450  IF( linfo.NE.0 )
451  $ GO TO 70
452 *
453 * Compute orthogonal matrix QL:
454 *
455 * QL**T * LI = [ TL ]
456 * [ 0 ]
457 * where
458 * LI = [ -L ]
459 * [ SCALE * identity(N2) ]
460 *
461  DO 10 i = 1, n2
462  CALL sscal( n1, -one, li( 1, i ), 1 )
463  li( n1+i, i ) = scale
464  10 CONTINUE
465  CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
466  IF( linfo.NE.0 )
467  $ GO TO 70
468  CALL sorg2r( m, m, n2, li, ldst, taul, work, linfo )
469  IF( linfo.NE.0 )
470  $ GO TO 70
471 *
472 * Compute orthogonal matrix RQ:
473 *
474 * IR * RQ**T = [ 0 TR],
475 *
476 * where IR = [ SCALE * identity(N1), R ]
477 *
478  DO 20 i = 1, n1
479  ir( n2+i, i ) = scale
480  20 CONTINUE
481  CALL sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
482  IF( linfo.NE.0 )
483  $ GO TO 70
484  CALL sorgr2( m, m, n1, ir, ldst, taur, work, linfo )
485  IF( linfo.NE.0 )
486  $ GO TO 70
487 *
488 * Perform the swapping tentatively:
489 *
490  CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
491  $ work, m )
492  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
493  $ ldst )
494  CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
495  $ work, m )
496  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
497  $ ldst )
498  CALL slacpy( 'F', m, m, s, ldst, scpy, ldst )
499  CALL slacpy( 'F', m, m, t, ldst, tcpy, ldst )
500  CALL slacpy( 'F', m, m, ir, ldst, ircop, ldst )
501  CALL slacpy( 'F', m, m, li, ldst, licop, ldst )
502 *
503 * Triangularize the B-part by an RQ factorization.
504 * Apply transformation (from left) to A-part, giving S.
505 *
506  CALL sgerq2( m, m, t, ldst, taur, work, linfo )
507  IF( linfo.NE.0 )
508  $ GO TO 70
509  CALL sormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
510  $ linfo )
511  IF( linfo.NE.0 )
512  $ GO TO 70
513  CALL sormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
514  $ linfo )
515  IF( linfo.NE.0 )
516  $ GO TO 70
517 *
518 * Compute F-norm(S21) in BRQA21. (T21 is 0.)
519 *
520  dscale = zero
521  dsum = one
522  DO 30 i = 1, n2
523  CALL slassq( n1, s( n2+1, i ), 1, dscale, dsum )
524  30 CONTINUE
525  brqa21 = dscale*sqrt( dsum )
526 *
527 * Triangularize the B-part by a QR factorization.
528 * Apply transformation (from right) to A-part, giving S.
529 *
530  CALL sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
531  IF( linfo.NE.0 )
532  $ GO TO 70
533  CALL sorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
534  $ work, info )
535  CALL sorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
536  $ work, info )
537  IF( linfo.NE.0 )
538  $ GO TO 70
539 *
540 * Compute F-norm(S21) in BQRA21. (T21 is 0.)
541 *
542  dscale = zero
543  dsum = one
544  DO 40 i = 1, n2
545  CALL slassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
546  40 CONTINUE
547  bqra21 = dscale*sqrt( dsum )
548 *
549 * Decide which method to use.
550 * Weak stability test:
551 * F-norm(S21) <= O(EPS * F-norm((S)))
552 *
553  IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresha ) THEN
554  CALL slacpy( 'F', m, m, scpy, ldst, s, ldst )
555  CALL slacpy( 'F', m, m, tcpy, ldst, t, ldst )
556  CALL slacpy( 'F', m, m, ircop, ldst, ir, ldst )
557  CALL slacpy( 'F', m, m, licop, ldst, li, ldst )
558  ELSE IF( brqa21.GE.thresha ) THEN
559  GO TO 70
560  END IF
561 *
562 * Set lower triangle of B-part to zero
563 *
564  CALL slaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
565 *
566  IF( wands ) THEN
567 *
568 * Strong stability test:
569 * F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
570 * and
571 * F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
572 *
573  CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
574  $ m )
575  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
576  $ work, m )
577  CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
578  $ work( m*m+1 ), m )
579  dscale = zero
580  dsum = one
581  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
582  sa = dscale*sqrt( dsum )
583 *
584  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
585  $ m )
586  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
587  $ work, m )
588  CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
589  $ work( m*m+1 ), m )
590  dscale = zero
591  dsum = one
592  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
593  sb = dscale*sqrt( dsum )
594  strong = sa.LE.thresha .AND. sb.LE.threshb
595  IF( .NOT.strong )
596  $ GO TO 70
597 *
598  END IF
599 *
600 * If the swap is accepted ("weakly" and "strongly"), apply the
601 * transformations and set N1-by-N2 (2,1)-block to zero.
602 *
603  CALL slaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
604 *
605 * copy back M-by-M diagonal block starting at index J1 of (A, B)
606 *
607  CALL slacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
608  CALL slacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
609  CALL slaset( 'Full', ldst, ldst, zero, zero, t, ldst )
610 *
611 * Standardize existing 2-by-2 blocks.
612 *
613  CALL slaset( 'Full', m, m, zero, zero, work, m )
614  work( 1 ) = one
615  t( 1, 1 ) = one
616  idum = lwork - m*m - 2
617  IF( n2.GT.1 ) THEN
618  CALL slagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
619  $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
620  work( m+1 ) = -work( 2 )
621  work( m+2 ) = work( 1 )
622  t( n2, n2 ) = t( 1, 1 )
623  t( 1, 2 ) = -t( 2, 1 )
624  END IF
625  work( m*m ) = one
626  t( m, m ) = one
627 *
628  IF( n1.GT.1 ) THEN
629  CALL slagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
630  $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
631  $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
632  $ t( m, m-1 ) )
633  work( m*m ) = work( n2*m+n2+1 )
634  work( m*m-1 ) = -work( n2*m+n2+2 )
635  t( m, m ) = t( n2+1, n2+1 )
636  t( m-1, m ) = -t( m, m-1 )
637  END IF
638  CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
639  $ lda, zero, work( m*m+1 ), n2 )
640  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
641  $ lda )
642  CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
643  $ ldb, zero, work( m*m+1 ), n2 )
644  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
645  $ ldb )
646  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
647  $ work( m*m+1 ), m )
648  CALL slacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
649  CALL sgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
650  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
651  CALL slacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
652  CALL sgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
653  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
654  CALL slacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
655  CALL sgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
656  $ work, m )
657  CALL slacpy( 'Full', m, m, work, m, ir, ldst )
658 *
659 * Accumulate transformations into Q and Z if requested.
660 *
661  IF( wantq ) THEN
662  CALL sgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
663  $ ldst, zero, work, n )
664  CALL slacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
665 *
666  END IF
667 *
668  IF( wantz ) THEN
669  CALL sgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
670  $ ldst, zero, work, n )
671  CALL slacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
672 *
673  END IF
674 *
675 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
676 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
677 *
678  i = j1 + m
679  IF( i.LE.n ) THEN
680  CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
681  $ a( j1, i ), lda, zero, work, m )
682  CALL slacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
683  CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
684  $ b( j1, i ), ldb, zero, work, m )
685  CALL slacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
686  END IF
687  i = j1 - 1
688  IF( i.GT.0 ) THEN
689  CALL sgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
690  $ ldst, zero, work, i )
691  CALL slacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
692  CALL sgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
693  $ ldst, zero, work, i )
694  CALL slacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
695  END IF
696 *
697 * Exit with INFO = 0 if swap was successfully performed.
698 *
699  RETURN
700 *
701  END IF
702 *
703 * Exit with INFO = 1 if swap was rejected.
704 *
705  70 CONTINUE
706 *
707  info = 1
708  RETURN
709 *
710 * End of STGEX2
711 *
712  END
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f90:126
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:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:113
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:221
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:123
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:130
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:157
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:159
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: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:114
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:114
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:274
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187