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