Annotation of rpl/lapack/lapack/dlaexc.f, revision 1.15
1.12 bertrand 1: *> \brief \b DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAEXC + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
22: * INFO )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL WANTQ
26: * INTEGER INFO, J1, LDQ, LDT, N, N1, N2
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
39: *> an upper quasi-triangular matrix T by an orthogonal similarity
40: *> transformation.
41: *>
42: *> T must be in Schur canonical form, that is, block upper triangular
43: *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
44: *> has its diagonal elemnts equal and its off-diagonal elements of
45: *> opposite sign.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] WANTQ
52: *> \verbatim
53: *> WANTQ is LOGICAL
54: *> = .TRUE. : accumulate the transformation in the matrix Q;
55: *> = .FALSE.: do not accumulate the transformation.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix T. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in,out] T
65: *> \verbatim
66: *> T is DOUBLE PRECISION array, dimension (LDT,N)
67: *> On entry, the upper quasi-triangular matrix T, in Schur
68: *> canonical form.
69: *> On exit, the updated matrix T, again in Schur canonical form.
70: *> \endverbatim
71: *>
72: *> \param[in] LDT
73: *> \verbatim
74: *> LDT is INTEGER
75: *> The leading dimension of the array T. LDT >= max(1,N).
76: *> \endverbatim
77: *>
78: *> \param[in,out] Q
79: *> \verbatim
80: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
81: *> On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
82: *> On exit, if WANTQ is .TRUE., the updated matrix Q.
83: *> If WANTQ is .FALSE., Q is not referenced.
84: *> \endverbatim
85: *>
86: *> \param[in] LDQ
87: *> \verbatim
88: *> LDQ is INTEGER
89: *> The leading dimension of the array Q.
90: *> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
91: *> \endverbatim
92: *>
93: *> \param[in] J1
94: *> \verbatim
95: *> J1 is INTEGER
96: *> The index of the first row of the first block T11.
97: *> \endverbatim
98: *>
99: *> \param[in] N1
100: *> \verbatim
101: *> N1 is INTEGER
102: *> The order of the first block T11. N1 = 0, 1 or 2.
103: *> \endverbatim
104: *>
105: *> \param[in] N2
106: *> \verbatim
107: *> N2 is INTEGER
108: *> The order of the second block T22. N2 = 0, 1 or 2.
109: *> \endverbatim
110: *>
111: *> \param[out] WORK
112: *> \verbatim
113: *> WORK is DOUBLE PRECISION array, dimension (N)
114: *> \endverbatim
115: *>
116: *> \param[out] INFO
117: *> \verbatim
118: *> INFO is INTEGER
119: *> = 0: successful exit
120: *> = 1: the transformed matrix T would be too far from Schur
121: *> form; the blocks are not swapped and T and Q are
122: *> unchanged.
123: *> \endverbatim
124: *
125: * Authors:
126: * ========
127: *
128: *> \author Univ. of Tennessee
129: *> \author Univ. of California Berkeley
130: *> \author Univ. of Colorado Denver
131: *> \author NAG Ltd.
132: *
1.12 bertrand 133: *> \date September 2012
1.9 bertrand 134: *
135: *> \ingroup doubleOTHERauxiliary
136: *
137: * =====================================================================
1.1 bertrand 138: SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
139: $ INFO )
140: *
1.12 bertrand 141: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 142: * -- LAPACK is a software package provided by Univ. of Tennessee, --
143: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 bertrand 144: * September 2012
1.1 bertrand 145: *
146: * .. Scalar Arguments ..
147: LOGICAL WANTQ
148: INTEGER INFO, J1, LDQ, LDT, N, N1, N2
149: * ..
150: * .. Array Arguments ..
151: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
152: * ..
153: *
154: * =====================================================================
155: *
156: * .. Parameters ..
157: DOUBLE PRECISION ZERO, ONE
158: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
159: DOUBLE PRECISION TEN
160: PARAMETER ( TEN = 1.0D+1 )
161: INTEGER LDD, LDX
162: PARAMETER ( LDD = 4, LDX = 2 )
163: * ..
164: * .. Local Scalars ..
165: INTEGER IERR, J2, J3, J4, K, ND
166: DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
167: $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
168: $ WR1, WR2, XNORM
169: * ..
170: * .. Local Arrays ..
171: DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
172: $ X( LDX, 2 )
173: * ..
174: * .. External Functions ..
175: DOUBLE PRECISION DLAMCH, DLANGE
176: EXTERNAL DLAMCH, DLANGE
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
180: $ DROT
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC ABS, MAX
184: * ..
185: * .. Executable Statements ..
186: *
187: INFO = 0
188: *
189: * Quick return if possible
190: *
191: IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
192: $ RETURN
193: IF( J1+N1.GT.N )
194: $ RETURN
195: *
196: J2 = J1 + 1
197: J3 = J1 + 2
198: J4 = J1 + 3
199: *
200: IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
201: *
202: * Swap two 1-by-1 blocks.
203: *
204: T11 = T( J1, J1 )
205: T22 = T( J2, J2 )
206: *
207: * Determine the transformation to perform the interchange.
208: *
209: CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
210: *
211: * Apply transformation to the matrix T.
212: *
213: IF( J3.LE.N )
214: $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
215: $ SN )
216: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
217: *
218: T( J1, J1 ) = T22
219: T( J2, J2 ) = T11
220: *
221: IF( WANTQ ) THEN
222: *
223: * Accumulate transformation in the matrix Q.
224: *
225: CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
226: END IF
227: *
228: ELSE
229: *
230: * Swapping involves at least one 2-by-2 block.
231: *
232: * Copy the diagonal block of order N1+N2 to the local array D
233: * and compute its norm.
234: *
235: ND = N1 + N2
236: CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
237: DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
238: *
239: * Compute machine-dependent threshold for test for accepting
240: * swap.
241: *
242: EPS = DLAMCH( 'P' )
243: SMLNUM = DLAMCH( 'S' ) / EPS
244: THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
245: *
246: * Solve T11*X - X*T22 = scale*T12 for X.
247: *
248: CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
249: $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
250: $ LDX, XNORM, IERR )
251: *
252: * Swap the adjacent diagonal blocks.
253: *
254: K = N1 + N1 + N2 - 3
255: GO TO ( 10, 20, 30 )K
256: *
257: 10 CONTINUE
258: *
259: * N1 = 1, N2 = 2: generate elementary reflector H so that:
260: *
261: * ( scale, X11, X12 ) H = ( 0, 0, * )
262: *
263: U( 1 ) = SCALE
264: U( 2 ) = X( 1, 1 )
265: U( 3 ) = X( 1, 2 )
266: CALL DLARFG( 3, U( 3 ), U, 1, TAU )
267: U( 3 ) = ONE
268: T11 = T( J1, J1 )
269: *
270: * Perform swap provisionally on diagonal block in D.
271: *
272: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
273: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
274: *
275: * Test whether to reject swap.
276: *
277: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
278: $ 3 )-T11 ) ).GT.THRESH )GO TO 50
279: *
280: * Accept swap: apply transformation to the entire matrix T.
281: *
282: CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
283: CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
284: *
285: T( J3, J1 ) = ZERO
286: T( J3, J2 ) = ZERO
287: T( J3, J3 ) = T11
288: *
289: IF( WANTQ ) THEN
290: *
291: * Accumulate transformation in the matrix Q.
292: *
293: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
294: END IF
295: GO TO 40
296: *
297: 20 CONTINUE
298: *
299: * N1 = 2, N2 = 1: generate elementary reflector H so that:
300: *
301: * H ( -X11 ) = ( * )
302: * ( -X21 ) = ( 0 )
303: * ( scale ) = ( 0 )
304: *
305: U( 1 ) = -X( 1, 1 )
306: U( 2 ) = -X( 2, 1 )
307: U( 3 ) = SCALE
308: CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
309: U( 1 ) = ONE
310: T33 = T( J3, J3 )
311: *
312: * Perform swap provisionally on diagonal block in D.
313: *
314: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
315: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
316: *
317: * Test whether to reject swap.
318: *
319: IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
320: $ 1 )-T33 ) ).GT.THRESH )GO TO 50
321: *
322: * Accept swap: apply transformation to the entire matrix T.
323: *
324: CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
325: CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
326: *
327: T( J1, J1 ) = T33
328: T( J2, J1 ) = ZERO
329: T( J3, J1 ) = ZERO
330: *
331: IF( WANTQ ) THEN
332: *
333: * Accumulate transformation in the matrix Q.
334: *
335: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
336: END IF
337: GO TO 40
338: *
339: 30 CONTINUE
340: *
341: * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
342: * that:
343: *
344: * H(2) H(1) ( -X11 -X12 ) = ( * * )
345: * ( -X21 -X22 ) ( 0 * )
346: * ( scale 0 ) ( 0 0 )
347: * ( 0 scale ) ( 0 0 )
348: *
349: U1( 1 ) = -X( 1, 1 )
350: U1( 2 ) = -X( 2, 1 )
351: U1( 3 ) = SCALE
352: CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
353: U1( 1 ) = ONE
354: *
355: TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
356: U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
357: U2( 2 ) = -TEMP*U1( 3 )
358: U2( 3 ) = SCALE
359: CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
360: U2( 1 ) = ONE
361: *
362: * Perform swap provisionally on diagonal block in D.
363: *
364: CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
365: CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
366: CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
367: CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
368: *
369: * Test whether to reject swap.
370: *
371: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
372: $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
373: *
374: * Accept swap: apply transformation to the entire matrix T.
375: *
376: CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
377: CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
378: CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
379: CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
380: *
381: T( J3, J1 ) = ZERO
382: T( J3, J2 ) = ZERO
383: T( J4, J1 ) = ZERO
384: T( J4, J2 ) = ZERO
385: *
386: IF( WANTQ ) THEN
387: *
388: * Accumulate transformation in the matrix Q.
389: *
390: CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
391: CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
392: END IF
393: *
394: 40 CONTINUE
395: *
396: IF( N2.EQ.2 ) THEN
397: *
398: * Standardize new 2-by-2 block T11
399: *
400: CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
401: $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
402: CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
403: $ CS, SN )
404: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
405: IF( WANTQ )
406: $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
407: END IF
408: *
409: IF( N1.EQ.2 ) THEN
410: *
411: * Standardize new 2-by-2 block T22
412: *
413: J3 = J1 + N2
414: J4 = J3 + 1
415: CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
416: $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
417: IF( J3+2.LE.N )
418: $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
419: $ LDT, CS, SN )
420: CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
421: IF( WANTQ )
422: $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
423: END IF
424: *
425: END IF
426: RETURN
427: *
428: * Exit with INFO = 1 if swap was rejected.
429: *
430: 50 CONTINUE
431: INFO = 1
432: RETURN
433: *
434: * End of DLAEXC
435: *
436: 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