001:       SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
002:      $                   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
011:       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
021: *  an upper quasi-triangular matrix T by an orthogonal similarity
022: *  transformation.
023: *
024: *  T must be in Schur canonical form, that is, block upper triangular
025: *  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
026: *  has its diagonal elemnts equal and its off-diagonal elements of
027: *  opposite sign.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  WANTQ   (input) LOGICAL
033: *          = .TRUE. : accumulate the transformation in the matrix Q;
034: *          = .FALSE.: do not accumulate the transformation.
035: *
036: *  N       (input) INTEGER
037: *          The order of the matrix T. N >= 0.
038: *
039: *  T       (input/output) REAL array, dimension (LDT,N)
040: *          On entry, the upper quasi-triangular matrix T, in Schur
041: *          canonical form.
042: *          On exit, the updated matrix T, again in Schur canonical form.
043: *
044: *  LDT     (input)  INTEGER
045: *          The leading dimension of the array T. LDT >= max(1,N).
046: *
047: *  Q       (input/output) REAL array, dimension (LDQ,N)
048: *          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
049: *          On exit, if WANTQ is .TRUE., the updated matrix Q.
050: *          If WANTQ is .FALSE., Q is not referenced.
051: *
052: *  LDQ     (input) INTEGER
053: *          The leading dimension of the array Q.
054: *          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
055: *
056: *  J1      (input) INTEGER
057: *          The index of the first row of the first block T11.
058: *
059: *  N1      (input) INTEGER
060: *          The order of the first block T11. N1 = 0, 1 or 2.
061: *
062: *  N2      (input) INTEGER
063: *          The order of the second block T22. N2 = 0, 1 or 2.
064: *
065: *  WORK    (workspace) REAL array, dimension (N)
066: *
067: *  INFO    (output) INTEGER
068: *          = 0: successful exit
069: *          = 1: the transformed matrix T would be too far from Schur
070: *               form; the blocks are not swapped and T and Q are
071: *               unchanged.
072: *
073: *  =====================================================================
074: *
075: *     .. Parameters ..
076:       REAL               ZERO, ONE
077:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
078:       REAL               TEN
079:       PARAMETER          ( TEN = 1.0E+1 )
080:       INTEGER            LDD, LDX
081:       PARAMETER          ( LDD = 4, LDX = 2 )
082: *     ..
083: *     .. Local Scalars ..
084:       INTEGER            IERR, J2, J3, J4, K, ND
085:       REAL               CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
086:      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
087:      $                   WR1, WR2, XNORM
088: *     ..
089: *     .. Local Arrays ..
090:       REAL               D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
091:      $                   X( LDX, 2 )
092: *     ..
093: *     .. External Functions ..
094:       REAL               SLAMCH, SLANGE
095:       EXTERNAL           SLAMCH, SLANGE
096: *     ..
097: *     .. External Subroutines ..
098:       EXTERNAL           SLACPY, SLANV2, SLARFG, SLARFX, SLARTG, SLASY2,
099:      $                   SROT
100: *     ..
101: *     .. Intrinsic Functions ..
102:       INTRINSIC          ABS, MAX
103: *     ..
104: *     .. Executable Statements ..
105: *
106:       INFO = 0
107: *
108: *     Quick return if possible
109: *
110:       IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
111:      $   RETURN
112:       IF( J1+N1.GT.N )
113:      $   RETURN
114: *
115:       J2 = J1 + 1
116:       J3 = J1 + 2
117:       J4 = J1 + 3
118: *
119:       IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
120: *
121: *        Swap two 1-by-1 blocks.
122: *
123:          T11 = T( J1, J1 )
124:          T22 = T( J2, J2 )
125: *
126: *        Determine the transformation to perform the interchange.
127: *
128:          CALL SLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
129: *
130: *        Apply transformation to the matrix T.
131: *
132:          IF( J3.LE.N )
133:      $      CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
134:      $                 SN )
135:          CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
136: *
137:          T( J1, J1 ) = T22
138:          T( J2, J2 ) = T11
139: *
140:          IF( WANTQ ) THEN
141: *
142: *           Accumulate transformation in the matrix Q.
143: *
144:             CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
145:          END IF
146: *
147:       ELSE
148: *
149: *        Swapping involves at least one 2-by-2 block.
150: *
151: *        Copy the diagonal block of order N1+N2 to the local array D
152: *        and compute its norm.
153: *
154:          ND = N1 + N2
155:          CALL SLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
156:          DNORM = SLANGE( 'Max', ND, ND, D, LDD, WORK )
157: *
158: *        Compute machine-dependent threshold for test for accepting
159: *        swap.
160: *
161:          EPS = SLAMCH( 'P' )
162:          SMLNUM = SLAMCH( 'S' ) / EPS
163:          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
164: *
165: *        Solve T11*X - X*T22 = scale*T12 for X.
166: *
167:          CALL SLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
168:      $                D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
169:      $                LDX, XNORM, IERR )
170: *
171: *        Swap the adjacent diagonal blocks.
172: *
173:          K = N1 + N1 + N2 - 3
174:          GO TO ( 10, 20, 30 )K
175: *
176:    10    CONTINUE
177: *
178: *        N1 = 1, N2 = 2: generate elementary reflector H so that:
179: *
180: *        ( scale, X11, X12 ) H = ( 0, 0, * )
181: *
182:          U( 1 ) = SCALE
183:          U( 2 ) = X( 1, 1 )
184:          U( 3 ) = X( 1, 2 )
185:          CALL SLARFG( 3, U( 3 ), U, 1, TAU )
186:          U( 3 ) = ONE
187:          T11 = T( J1, J1 )
188: *
189: *        Perform swap provisionally on diagonal block in D.
190: *
191:          CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
192:          CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
193: *
194: *        Test whether to reject swap.
195: *
196:          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
197:      $       3 )-T11 ) ).GT.THRESH )GO TO 50
198: *
199: *        Accept swap: apply transformation to the entire matrix T.
200: *
201:          CALL SLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
202:          CALL SLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
203: *
204:          T( J3, J1 ) = ZERO
205:          T( J3, J2 ) = ZERO
206:          T( J3, J3 ) = T11
207: *
208:          IF( WANTQ ) THEN
209: *
210: *           Accumulate transformation in the matrix Q.
211: *
212:             CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
213:          END IF
214:          GO TO 40
215: *
216:    20    CONTINUE
217: *
218: *        N1 = 2, N2 = 1: generate elementary reflector H so that:
219: *
220: *        H (  -X11 ) = ( * )
221: *          (  -X21 ) = ( 0 )
222: *          ( scale ) = ( 0 )
223: *
224:          U( 1 ) = -X( 1, 1 )
225:          U( 2 ) = -X( 2, 1 )
226:          U( 3 ) = SCALE
227:          CALL SLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
228:          U( 1 ) = ONE
229:          T33 = T( J3, J3 )
230: *
231: *        Perform swap provisionally on diagonal block in D.
232: *
233:          CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
234:          CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
235: *
236: *        Test whether to reject swap.
237: *
238:          IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
239:      $       1 )-T33 ) ).GT.THRESH )GO TO 50
240: *
241: *        Accept swap: apply transformation to the entire matrix T.
242: *
243:          CALL SLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
244:          CALL SLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
245: *
246:          T( J1, J1 ) = T33
247:          T( J2, J1 ) = ZERO
248:          T( J3, J1 ) = ZERO
249: *
250:          IF( WANTQ ) THEN
251: *
252: *           Accumulate transformation in the matrix Q.
253: *
254:             CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
255:          END IF
256:          GO TO 40
257: *
258:    30    CONTINUE
259: *
260: *        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
261: *        that:
262: *
263: *        H(2) H(1) (  -X11  -X12 ) = (  *  * )
264: *                  (  -X21  -X22 )   (  0  * )
265: *                  ( scale    0  )   (  0  0 )
266: *                  (    0  scale )   (  0  0 )
267: *
268:          U1( 1 ) = -X( 1, 1 )
269:          U1( 2 ) = -X( 2, 1 )
270:          U1( 3 ) = SCALE
271:          CALL SLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
272:          U1( 1 ) = ONE
273: *
274:          TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
275:          U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
276:          U2( 2 ) = -TEMP*U1( 3 )
277:          U2( 3 ) = SCALE
278:          CALL SLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
279:          U2( 1 ) = ONE
280: *
281: *        Perform swap provisionally on diagonal block in D.
282: *
283:          CALL SLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
284:          CALL SLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
285:          CALL SLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
286:          CALL SLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
287: *
288: *        Test whether to reject swap.
289: *
290:          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
291:      $       ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
292: *
293: *        Accept swap: apply transformation to the entire matrix T.
294: *
295:          CALL SLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
296:          CALL SLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
297:          CALL SLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
298:          CALL SLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
299: *
300:          T( J3, J1 ) = ZERO
301:          T( J3, J2 ) = ZERO
302:          T( J4, J1 ) = ZERO
303:          T( J4, J2 ) = ZERO
304: *
305:          IF( WANTQ ) THEN
306: *
307: *           Accumulate transformation in the matrix Q.
308: *
309:             CALL SLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
310:             CALL SLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
311:          END IF
312: *
313:    40    CONTINUE
314: *
315:          IF( N2.EQ.2 ) THEN
316: *
317: *           Standardize new 2-by-2 block T11
318: *
319:             CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
320:      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
321:             CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
322:      $                 CS, SN )
323:             CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
324:             IF( WANTQ )
325:      $         CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
326:          END IF
327: *
328:          IF( N1.EQ.2 ) THEN
329: *
330: *           Standardize new 2-by-2 block T22
331: *
332:             J3 = J1 + N2
333:             J4 = J3 + 1
334:             CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
335:      $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
336:             IF( J3+2.LE.N )
337:      $         CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
338:      $                    LDT, CS, SN )
339:             CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
340:             IF( WANTQ )
341:      $         CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
342:          END IF
343: *
344:       END IF
345:       RETURN
346: *
347: *     Exit with INFO = 1 if swap was rejected.
348: *
349:    50 INFO = 1
350:       RETURN
351: *
352: *     End of SLAEXC
353: *
354:       END
355: