001:       SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
002:      $                   LDZ, J1, N1, N2, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            WANTQ, WANTZ
011:       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
015:      $                   WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  STGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
022: *  of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
023: *  (A, B) by an orthogonal equivalence transformation.
024: *
025: *  (A, B) must be in generalized real Schur canonical form (as returned
026: *  by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
027: *  diagonal blocks. B is upper triangular.
028: *
029: *  Optionally, the matrices Q and Z of generalized Schur vectors are
030: *  updated.
031: *
032: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
033: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
034: *
035: *
036: *  Arguments
037: *  =========
038: *
039: *  WANTQ   (input) LOGICAL
040: *          .TRUE. : update the left transformation matrix Q;
041: *          .FALSE.: do not update Q.
042: *
043: *  WANTZ   (input) LOGICAL
044: *          .TRUE. : update the right transformation matrix Z;
045: *          .FALSE.: do not update Z.
046: *
047: *  N       (input) INTEGER
048: *          The order of the matrices A and B. N >= 0.
049: *
050: *  A      (input/output) REAL arrays, dimensions (LDA,N)
051: *          On entry, the matrix A in the pair (A, B).
052: *          On exit, the updated matrix A.
053: *
054: *  LDA     (input)  INTEGER
055: *          The leading dimension of the array A. LDA >= max(1,N).
056: *
057: *  B      (input/output) REAL arrays, dimensions (LDB,N)
058: *          On entry, the matrix B in the pair (A, B).
059: *          On exit, the updated matrix B.
060: *
061: *  LDB     (input)  INTEGER
062: *          The leading dimension of the array B. LDB >= max(1,N).
063: *
064: *  Q       (input/output) REAL array, dimension (LDZ,N)
065: *          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
066: *          On exit, the updated matrix Q.
067: *          Not referenced if WANTQ = .FALSE..
068: *
069: *  LDQ     (input) INTEGER
070: *          The leading dimension of the array Q. LDQ >= 1.
071: *          If WANTQ = .TRUE., LDQ >= N.
072: *
073: *  Z       (input/output) REAL array, dimension (LDZ,N)
074: *          On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
075: *          On exit, the updated matrix Z.
076: *          Not referenced if WANTZ = .FALSE..
077: *
078: *  LDZ     (input) INTEGER
079: *          The leading dimension of the array Z. LDZ >= 1.
080: *          If WANTZ = .TRUE., LDZ >= N.
081: *
082: *  J1      (input) INTEGER
083: *          The index to the first block (A11, B11). 1 <= J1 <= N.
084: *
085: *  N1      (input) INTEGER
086: *          The order of the first block (A11, B11). N1 = 0, 1 or 2.
087: *
088: *  N2      (input) INTEGER
089: *          The order of the second block (A22, B22). N2 = 0, 1 or 2.
090: *
091: *  WORK    (workspace) REAL array, dimension (MAX(1,LWORK)).
092: *
093: *  LWORK   (input) INTEGER
094: *          The dimension of the array WORK.
095: *          LWORK >=  MAX( N*(N2+N1), (N2+N1)*(N2+N1)*2 )
096: *
097: *  INFO    (output) INTEGER
098: *            =0: Successful exit
099: *            >0: If INFO = 1, the transformed matrix (A, B) would be
100: *                too far from generalized Schur form; the blocks are
101: *                not swapped and (A, B) and (Q, Z) are unchanged.
102: *                The problem of swapping is too ill-conditioned.
103: *            <0: If INFO = -16: LWORK is too small. Appropriate value
104: *                for LWORK is returned in WORK(1).
105: *
106: *  Further Details
107: *  ===============
108: *
109: *  Based on contributions by
110: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
111: *     Umea University, S-901 87 Umea, Sweden.
112: *
113: *  In the current code both weak and strong stability tests are
114: *  performed. The user can omit the strong stability test by changing
115: *  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
116: *  details.
117: *
118: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
119: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
120: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
121: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
122: *
123: *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
124: *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
125: *      Estimation: Theory, Algorithms and Software,
126: *      Report UMINF - 94.04, Department of Computing Science, Umea
127: *      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
128: *      Note 87. To appear in Numerical Algorithms, 1996.
129: *
130: *  =====================================================================
131: *  Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
132: *  loops. Sven Hammarling, 1/5/02.
133: *
134: *     .. Parameters ..
135:       REAL               ZERO, ONE
136:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
137:       REAL               TEN
138:       PARAMETER          ( TEN = 1.0E+01 )
139:       INTEGER            LDST
140:       PARAMETER          ( LDST = 4 )
141:       LOGICAL            WANDS
142:       PARAMETER          ( WANDS = .TRUE. )
143: *     ..
144: *     .. Local Scalars ..
145:       LOGICAL            STRONG, WEAK
146:       INTEGER            I, IDUM, LINFO, M
147:       REAL               BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
148:      $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
149: *     ..
150: *     .. Local Arrays ..
151:       INTEGER            IWORK( LDST )
152:       REAL               AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
153:      $                   IRCOP( LDST, LDST ), LI( LDST, LDST ),
154:      $                   LICOP( LDST, LDST ), S( LDST, LDST ),
155:      $                   SCPY( LDST, LDST ), T( LDST, LDST ),
156:      $                   TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
157: *     ..
158: *     .. External Functions ..
159:       REAL               SLAMCH
160:       EXTERNAL           SLAMCH
161: *     ..
162: *     .. External Subroutines ..
163:       EXTERNAL           SGEMM, SGEQR2, SGERQ2, SLACPY, SLAGV2, SLARTG,
164:      $                   SLASET, SLASSQ, SORG2R, SORGR2, SORM2R, SORMR2,
165:      $                   SROT, SSCAL, STGSY2
166: *     ..
167: *     .. Intrinsic Functions ..
168:       INTRINSIC          ABS, MAX, SQRT
169: *     ..
170: *     .. Executable Statements ..
171: *
172:       INFO = 0
173: *
174: *     Quick return if possible
175: *
176:       IF( N.LE.1 .OR. N1.LE.0 .OR. N2.LE.0 )
177:      $   RETURN
178:       IF( N1.GT.N .OR. ( J1+N1 ).GT.N )
179:      $   RETURN
180:       M = N1 + N2
181:       IF( LWORK.LT.MAX( N*M, M*M*2 ) ) THEN
182:          INFO = -16
183:          WORK( 1 ) = MAX( N*M, M*M*2 )
184:          RETURN
185:       END IF
186: *
187:       WEAK = .FALSE.
188:       STRONG = .FALSE.
189: *
190: *     Make a local copy of selected block
191: *
192:       CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, LI, LDST )
193:       CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, IR, LDST )
194:       CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
195:       CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
196: *
197: *     Compute threshold for testing acceptance of swapping.
198: *
199:       EPS = SLAMCH( 'P' )
200:       SMLNUM = SLAMCH( 'S' ) / EPS
201:       DSCALE = ZERO
202:       DSUM = ONE
203:       CALL SLACPY( 'Full', M, M, S, LDST, WORK, M )
204:       CALL SLASSQ( M*M, WORK, 1, DSCALE, DSUM )
205:       CALL SLACPY( 'Full', M, M, T, LDST, WORK, M )
206:       CALL SLASSQ( M*M, WORK, 1, DSCALE, DSUM )
207:       DNORM = DSCALE*SQRT( DSUM )
208:       THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
209: *
210:       IF( M.EQ.2 ) THEN
211: *
212: *        CASE 1: Swap 1-by-1 and 1-by-1 blocks.
213: *
214: *        Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
215: *        using Givens rotations and perform the swap tentatively.
216: *
217:          F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
218:          G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
219:          SB = ABS( T( 2, 2 ) )
220:          SA = ABS( S( 2, 2 ) )
221:          CALL SLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
222:          IR( 2, 1 ) = -IR( 1, 2 )
223:          IR( 2, 2 ) = IR( 1, 1 )
224:          CALL SROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ),
225:      $              IR( 2, 1 ) )
226:          CALL SROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ),
227:      $              IR( 2, 1 ) )
228:          IF( SA.GE.SB ) THEN
229:             CALL SLARTG( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
230:      $                   DDUM )
231:          ELSE
232:             CALL SLARTG( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
233:      $                   DDUM )
234:          END IF
235:          CALL SROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ),
236:      $              LI( 2, 1 ) )
237:          CALL SROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, LI( 1, 1 ),
238:      $              LI( 2, 1 ) )
239:          LI( 2, 2 ) = LI( 1, 1 )
240:          LI( 1, 2 ) = -LI( 2, 1 )
241: *
242: *        Weak stability test:
243: *           |S21| + |T21| <= O(EPS * F-norm((S, T)))
244: *
245:          WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
246:          WEAK = WS.LE.THRESH
247:          IF( .NOT.WEAK )
248:      $      GO TO 70
249: *
250:          IF( WANDS ) THEN
251: *
252: *           Strong stability test:
253: *             F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A,B)))
254: *
255:             CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
256:      $                   M )
257:             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
258:      $                  WORK, M )
259:             CALL SGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
260:      $                  WORK( M*M+1 ), M )
261:             DSCALE = ZERO
262:             DSUM = ONE
263:             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
264: *
265:             CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
266:      $                   M )
267:             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
268:      $                  WORK, M )
269:             CALL SGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
270:      $                  WORK( M*M+1 ), M )
271:             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
272:             SS = DSCALE*SQRT( DSUM )
273:             STRONG = SS.LE.THRESH
274:             IF( .NOT.STRONG )
275:      $         GO TO 70
276:          END IF
277: *
278: *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
279: *               (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
280: *
281:          CALL SROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ),
282:      $              IR( 2, 1 ) )
283:          CALL SROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ),
284:      $              IR( 2, 1 ) )
285:          CALL SROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA,
286:      $              LI( 1, 1 ), LI( 2, 1 ) )
287:          CALL SROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB,
288:      $              LI( 1, 1 ), LI( 2, 1 ) )
289: *
290: *        Set  N1-by-N2 (2,1) - blocks to ZERO.
291: *
292:          A( J1+1, J1 ) = ZERO
293:          B( J1+1, J1 ) = ZERO
294: *
295: *        Accumulate transformations into Q and Z if requested.
296: *
297:          IF( WANTZ )
298:      $      CALL SROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
299:      $                 IR( 2, 1 ) )
300:          IF( WANTQ )
301:      $      CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, LI( 1, 1 ),
302:      $                 LI( 2, 1 ) )
303: *
304: *        Exit with INFO = 0 if swap was successfully performed.
305: *
306:          RETURN
307: *
308:       ELSE
309: *
310: *        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
311: *                and 2-by-2 blocks.
312: *
313: *        Solve the generalized Sylvester equation
314: *                 S11 * R - L * S22 = SCALE * S12
315: *                 T11 * R - L * T22 = SCALE * T12
316: *        for R and L. Solutions in LI and IR.
317: *
318:          CALL SLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
319:          CALL SLACPY( 'Full', N1, N2, S( 1, N1+1 ), LDST,
320:      $                IR( N2+1, N1+1 ), LDST )
321:          CALL STGSY2( 'N', 0, N1, N2, S, LDST, S( N1+1, N1+1 ), LDST,
322:      $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
323:      $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
324:      $                LINFO )
325: *
326: *        Compute orthogonal matrix QL:
327: *
328: *                    QL' * LI = [ TL ]
329: *                               [ 0  ]
330: *        where
331: *                    LI =  [      -L              ]
332: *                          [ SCALE * identity(N2) ]
333: *
334:          DO 10 I = 1, N2
335:             CALL SSCAL( N1, -ONE, LI( 1, I ), 1 )
336:             LI( N1+I, I ) = SCALE
337:    10    CONTINUE
338:          CALL SGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO )
339:          IF( LINFO.NE.0 )
340:      $      GO TO 70
341:          CALL SORG2R( M, M, N2, LI, LDST, TAUL, WORK, LINFO )
342:          IF( LINFO.NE.0 )
343:      $      GO TO 70
344: *
345: *        Compute orthogonal matrix RQ:
346: *
347: *                    IR * RQ' =   [ 0  TR],
348: *
349: *         where IR = [ SCALE * identity(N1), R ]
350: *
351:          DO 20 I = 1, N1
352:             IR( N2+I, I ) = SCALE
353:    20    CONTINUE
354:          CALL SGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
355:          IF( LINFO.NE.0 )
356:      $      GO TO 70
357:          CALL SORGR2( M, M, N1, IR, LDST, TAUR, WORK, LINFO )
358:          IF( LINFO.NE.0 )
359:      $      GO TO 70
360: *
361: *        Perform the swapping tentatively:
362: *
363:          CALL SGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
364:      $               WORK, M )
365:          CALL SGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, S,
366:      $               LDST )
367:          CALL SGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
368:      $               WORK, M )
369:          CALL SGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, T,
370:      $               LDST )
371:          CALL SLACPY( 'F', M, M, S, LDST, SCPY, LDST )
372:          CALL SLACPY( 'F', M, M, T, LDST, TCPY, LDST )
373:          CALL SLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
374:          CALL SLACPY( 'F', M, M, LI, LDST, LICOP, LDST )
375: *
376: *        Triangularize the B-part by an RQ factorization.
377: *        Apply transformation (from left) to A-part, giving S.
378: *
379:          CALL SGERQ2( M, M, T, LDST, TAUR, WORK, LINFO )
380:          IF( LINFO.NE.0 )
381:      $      GO TO 70
382:          CALL SORMR2( 'R', 'T', M, M, M, T, LDST, TAUR, S, LDST, WORK,
383:      $                LINFO )
384:          IF( LINFO.NE.0 )
385:      $      GO TO 70
386:          CALL SORMR2( 'L', 'N', M, M, M, T, LDST, TAUR, IR, LDST, WORK,
387:      $                LINFO )
388:          IF( LINFO.NE.0 )
389:      $      GO TO 70
390: *
391: *        Compute F-norm(S21) in BRQA21. (T21 is 0.)
392: *
393:          DSCALE = ZERO
394:          DSUM = ONE
395:          DO 30 I = 1, N2
396:             CALL SLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
397:    30    CONTINUE
398:          BRQA21 = DSCALE*SQRT( DSUM )
399: *
400: *        Triangularize the B-part by a QR factorization.
401: *        Apply transformation (from right) to A-part, giving S.
402: *
403:          CALL SGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO )
404:          IF( LINFO.NE.0 )
405:      $      GO TO 70
406:          CALL SORM2R( 'L', 'T', M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
407:      $                WORK, INFO )
408:          CALL SORM2R( 'R', 'N', M, M, M, TCPY, LDST, TAUL, LICOP, LDST,
409:      $                WORK, INFO )
410:          IF( LINFO.NE.0 )
411:      $      GO TO 70
412: *
413: *        Compute F-norm(S21) in BQRA21. (T21 is 0.)
414: *
415:          DSCALE = ZERO
416:          DSUM = ONE
417:          DO 40 I = 1, N2
418:             CALL SLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
419:    40    CONTINUE
420:          BQRA21 = DSCALE*SQRT( DSUM )
421: *
422: *        Decide which method to use.
423: *          Weak stability test:
424: *             F-norm(S21) <= O(EPS * F-norm((S, T)))
425: *
426:          IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
427:             CALL SLACPY( 'F', M, M, SCPY, LDST, S, LDST )
428:             CALL SLACPY( 'F', M, M, TCPY, LDST, T, LDST )
429:             CALL SLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
430:             CALL SLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
431:          ELSE IF( BRQA21.GE.THRESH ) THEN
432:             GO TO 70
433:          END IF
434: *
435: *        Set lower triangle of B-part to zero
436: *
437:          CALL SLASET( 'Lower', M-1, M-1, ZERO, ZERO, T(2,1), LDST )
438: *
439:          IF( WANDS ) THEN
440: *
441: *           Strong stability test:
442: *              F-norm((A-QL*S*QR', B-QL*T*QR')) <= O(EPS*F-norm((A,B)))
443: *
444:             CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
445:      $                   M )
446:             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
447:      $                  WORK, M )
448:             CALL SGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
449:      $                  WORK( M*M+1 ), M )
450:             DSCALE = ZERO
451:             DSUM = ONE
452:             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
453: *
454:             CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
455:      $                   M )
456:             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
457:      $                  WORK, M )
458:             CALL SGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
459:      $                  WORK( M*M+1 ), M )
460:             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
461:             SS = DSCALE*SQRT( DSUM )
462:             STRONG = ( SS.LE.THRESH )
463:             IF( .NOT.STRONG )
464:      $         GO TO 70
465: *
466:          END IF
467: *
468: *        If the swap is accepted ("weakly" and "strongly"), apply the
469: *        transformations and set N1-by-N2 (2,1)-block to zero.
470: *
471:          CALL SLASET( 'Full', N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
472: *
473: *        copy back M-by-M diagonal block starting at index J1 of (A, B)
474: *
475:          CALL SLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA )
476:          CALL SLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB )
477:          CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, T, LDST )
478: *
479: *        Standardize existing 2-by-2 blocks.
480: *
481:          DO 50 I = 1, M*M
482:             WORK(I) = ZERO
483:    50    CONTINUE
484:          WORK( 1 ) = ONE
485:          T( 1, 1 ) = ONE
486:          IDUM = LWORK - M*M - 2
487:          IF( N2.GT.1 ) THEN
488:             CALL SLAGV2( A( J1, J1 ), LDA, B( J1, J1 ), LDB, AR, AI, BE,
489:      $                   WORK( 1 ), WORK( 2 ), T( 1, 1 ), T( 2, 1 ) )
490:             WORK( M+1 ) = -WORK( 2 )
491:             WORK( M+2 ) = WORK( 1 )
492:             T( N2, N2 ) = T( 1, 1 )
493:             T( 1, 2 ) = -T( 2, 1 )
494:          END IF
495:          WORK( M*M ) = ONE
496:          T( M, M ) = ONE
497: *
498:          IF( N1.GT.1 ) THEN
499:             CALL SLAGV2( A( J1+N2, J1+N2 ), LDA, B( J1+N2, J1+N2 ), LDB,
500:      $                   TAUR, TAUL, WORK( M*M+1 ), WORK( N2*M+N2+1 ),
501:      $                   WORK( N2*M+N2+2 ), T( N2+1, N2+1 ),
502:      $                   T( M, M-1 ) )
503:             WORK( M*M ) = WORK( N2*M+N2+1 )
504:             WORK( M*M-1 ) = -WORK( N2*M+N2+2 )
505:             T( M, M ) = T( N2+1, N2+1 )
506:             T( M-1, M ) = -T( M, M-1 )
507:          END IF
508:          CALL SGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, A( J1, J1+N2 ),
509:      $               LDA, ZERO, WORK( M*M+1 ), N2 )
510:          CALL SLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ),
511:      $                LDA )
512:          CALL SGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, B( J1, J1+N2 ),
513:      $               LDB, ZERO, WORK( M*M+1 ), N2 )
514:          CALL SLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ),
515:      $                LDB )
516:          CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, WORK, M, ZERO,
517:      $               WORK( M*M+1 ), M )
518:          CALL SLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST )
519:          CALL SGEMM( 'N', 'N', N2, N1, N1, ONE, A( J1, J1+N2 ), LDA,
520:      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
521:          CALL SLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA )
522:          CALL SGEMM( 'N', 'N', N2, N1, N1, ONE, B( J1, J1+N2 ), LDB,
523:      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
524:          CALL SLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB )
525:          CALL SGEMM( 'T', 'N', M, M, M, ONE, IR, LDST, T, LDST, ZERO,
526:      $               WORK, M )
527:          CALL SLACPY( 'Full', M, M, WORK, M, IR, LDST )
528: *
529: *        Accumulate transformations into Q and Z if requested.
530: *
531:          IF( WANTQ ) THEN
532:             CALL SGEMM( 'N', 'N', N, M, M, ONE, Q( 1, J1 ), LDQ, LI,
533:      $                  LDST, ZERO, WORK, N )
534:             CALL SLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ )
535: *
536:          END IF
537: *
538:          IF( WANTZ ) THEN
539:             CALL SGEMM( 'N', 'N', N, M, M, ONE, Z( 1, J1 ), LDZ, IR,
540:      $                  LDST, ZERO, WORK, N )
541:             CALL SLACPY( 'Full', N, M, WORK, N, Z( 1, J1 ), LDZ )
542: *
543:          END IF
544: *
545: *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
546: *                (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
547: *
548:          I = J1 + M
549:          IF( I.LE.N ) THEN
550:             CALL SGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
551:      $                  A( J1, I ), LDA, ZERO, WORK, M )
552:             CALL SLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA )
553:             CALL SGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
554:      $                  B( J1, I ), LDB, ZERO, WORK, M )
555:             CALL SLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB )
556:          END IF
557:          I = J1 - 1
558:          IF( I.GT.0 ) THEN
559:             CALL SGEMM( 'N', 'N', I, M, M, ONE, A( 1, J1 ), LDA, IR,
560:      $                  LDST, ZERO, WORK, I )
561:             CALL SLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA )
562:             CALL SGEMM( 'N', 'N', I, M, M, ONE, B( 1, J1 ), LDB, IR,
563:      $                  LDST, ZERO, WORK, I )
564:             CALL SLACPY( 'Full', I, M, WORK, I, B( 1, J1 ), LDB )
565:          END IF
566: *
567: *        Exit with INFO = 0 if swap was successfully performed.
568: *
569:          RETURN
570: *
571:       END IF
572: *
573: *     Exit with INFO = 1 if swap was rejected.
574: *
575:    70 CONTINUE
576: *
577:       INFO = 1
578:       RETURN
579: *
580: *     End of STGEX2
581: *
582:       END
583: