1       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
  2      $                   INFO )
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.2.2) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *     June 2010
  8 *
  9 *     .. Scalar Arguments ..
 10       LOGICAL            WANTQ
 11       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
 21 *  an upper quasi-triangular matrix T by an orthogonal similarity
 22 *  transformation.
 23 *
 24 *  T must be in Schur canonical form, that is, block upper triangular
 25 *  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
 26 *  has its diagonal elemnts equal and its off-diagonal elements of
 27 *  opposite sign.
 28 *
 29 *  Arguments
 30 *  =========
 31 *
 32 *  WANTQ   (input) LOGICAL
 33 *          = .TRUE. : accumulate the transformation in the matrix Q;
 34 *          = .FALSE.: do not accumulate the transformation.
 35 *
 36 *  N       (input) INTEGER
 37 *          The order of the matrix T. N >= 0.
 38 *
 39 *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
 40 *          On entry, the upper quasi-triangular matrix T, in Schur
 41 *          canonical form.
 42 *          On exit, the updated matrix T, again in Schur canonical form.
 43 *
 44 *  LDT     (input) INTEGER
 45 *          The leading dimension of the array T. LDT >= max(1,N).
 46 *
 47 *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 48 *          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
 49 *          On exit, if WANTQ is .TRUE., the updated matrix Q.
 50 *          If WANTQ is .FALSE., Q is not referenced.
 51 *
 52 *  LDQ     (input) INTEGER
 53 *          The leading dimension of the array Q.
 54 *          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
 55 *
 56 *  J1      (input) INTEGER
 57 *          The index of the first row of the first block T11.
 58 *
 59 *  N1      (input) INTEGER
 60 *          The order of the first block T11. N1 = 0, 1 or 2.
 61 *
 62 *  N2      (input) INTEGER
 63 *          The order of the second block T22. N2 = 0, 1 or 2.
 64 *
 65 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
 66 *
 67 *  INFO    (output) INTEGER
 68 *          = 0: successful exit
 69 *          = 1: the transformed matrix T would be too far from Schur
 70 *               form; the blocks are not swapped and T and Q are
 71 *               unchanged.
 72 *
 73 *  =====================================================================
 74 *
 75 *     .. Parameters ..
 76       DOUBLE PRECISION   ZERO, ONE
 77       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 78       DOUBLE PRECISION   TEN
 79       PARAMETER          ( TEN = 1.0D+1 )
 80       INTEGER            LDD, LDX
 81       PARAMETER          ( LDD = 4, LDX = 2 )
 82 *     ..
 83 *     .. Local Scalars ..
 84       INTEGER            IERR, J2, J3, J4, K, ND
 85       DOUBLE PRECISION   CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
 86      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
 87      $                   WR1, WR2, XNORM
 88 *     ..
 89 *     .. Local Arrays ..
 90       DOUBLE PRECISION   D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
 91      $                   X( LDX, 2 )
 92 *     ..
 93 *     .. External Functions ..
 94       DOUBLE PRECISION   DLAMCH, DLANGE
 95       EXTERNAL           DLAMCH, DLANGE
 96 *     ..
 97 *     .. External Subroutines ..
 98       EXTERNAL           DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
 99      $                   DROT
100 *     ..
101 *     .. Intrinsic Functions ..
102       INTRINSIC          ABSMAX
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 DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
129 *
130 *        Apply transformation to the matrix T.
131 *
132          IF( J3.LE.N )
133      $      CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
134      $                 SN )
135          CALL DROT( 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 DROT( 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 DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
156          DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
157 *
158 *        Compute machine-dependent threshold for test for accepting
159 *        swap.
160 *
161          EPS = DLAMCH( 'P' )
162          SMLNUM = DLAMCH( 'S' ) / EPS
163          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
164 *
165 *        Solve T11*X - X*T22 = scale*T12 for X.
166 *
167          CALL DLASY2( .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 ( 102030 )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( 11 )
184          U( 3 ) = X( 12 )
185          CALL DLARFG( 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 DLARFX( 'L'33, U, TAU, D, LDD, WORK )
192          CALL DLARFX( 'R'33, U, TAU, D, LDD, WORK )
193 *
194 *        Test whether to reject swap.
195 *
196          IFMAXABS( D( 31 ) ), ABS( D( 32 ) ), 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 DLARFX( 'L'3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
202          CALL DLARFX( '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 DLARFX( '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( 11 )
225          U( 2 ) = -X( 21 )
226          U( 3 ) = SCALE
227          CALL DLARFG( 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 DLARFX( 'L'33, U, TAU, D, LDD, WORK )
234          CALL DLARFX( 'R'33, U, TAU, D, LDD, WORK )
235 *
236 *        Test whether to reject swap.
237 *
238          IFMAXABS( D( 21 ) ), ABS( D( 31 ) ), 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 DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
244          CALL DLARFX( '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 DLARFX( '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( 11 )
269          U1( 2 ) = -X( 21 )
270          U1( 3 ) = SCALE
271          CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
272          U1( 1 ) = ONE
273 *
274          TEMP = -TAU1*( X( 12 )+U1( 2 )*X( 22 ) )
275          U2( 1 ) = -TEMP*U1( 2 ) - X( 22 )
276          U2( 2 ) = -TEMP*U1( 3 )
277          U2( 3 ) = SCALE
278          CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
279          U2( 1 ) = ONE
280 *
281 *        Perform swap provisionally on diagonal block in D.
282 *
283          CALL DLARFX( 'L'34, U1, TAU1, D, LDD, WORK )
284          CALL DLARFX( 'R'43, U1, TAU1, D, LDD, WORK )
285          CALL DLARFX( 'L'34, U2, TAU2, D( 21 ), LDD, WORK )
286          CALL DLARFX( 'R'43, U2, TAU2, D( 12 ), LDD, WORK )
287 *
288 *        Test whether to reject swap.
289 *
290          IFMAXABS( D( 31 ) ), ABS( D( 32 ) ), ABS( D( 41 ) ),
291      $       ABS( D( 42 ) ) ).GT.THRESH )GO TO 50
292 *
293 *        Accept swap: apply transformation to the entire matrix T.
294 *
295          CALL DLARFX( 'L'3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
296          CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
297          CALL DLARFX( 'L'3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
298          CALL DLARFX( '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 DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
310             CALL DLARFX( '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 DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
320      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
321             CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
322      $                 CS, SN )
323             CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
324             IF( WANTQ )
325      $         CALL DROT( 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 DLANV2( 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 DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
338      $                    LDT, CS, SN )
339             CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
340             IF( WANTQ )
341      $         CALL DROT( 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 CONTINUE
350       INFO = 1
351       RETURN
352 *
353 *     End of DLAEXC
354 *
355       END