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 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 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 ( 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 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', 3, 3, U, TAU, D, LDD, WORK )
192: CALL DLARFX( '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 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( 1, 1 )
225: U( 2 ) = -X( 2, 1 )
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', 3, 3, U, TAU, D, LDD, WORK )
234: CALL DLARFX( '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 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( 1, 1 )
269: U1( 2 ) = -X( 2, 1 )
270: U1( 3 ) = SCALE
271: CALL DLARFG( 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 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', 3, 4, U1, TAU1, D, LDD, WORK )
284: CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
285: CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
286: CALL DLARFX( '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 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack