Annotation of rpl/lapack/lapack/dlaexc.f, revision 1.19
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: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
22: * INFO )
1.16 bertrand 23: *
1.9 bertrand 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: * ..
1.16 bertrand 31: *
1.9 bertrand 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
1.19 ! bertrand 44: *> has its diagonal elements equal and its off-diagonal elements of
1.9 bertrand 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: *
1.16 bertrand 128: *> \author Univ. of Tennessee
129: *> \author Univ. of California Berkeley
130: *> \author Univ. of Colorado Denver
131: *> \author NAG Ltd.
1.9 bertrand 132: *
133: *> \ingroup doubleOTHERauxiliary
134: *
135: * =====================================================================
1.1 bertrand 136: SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
137: $ INFO )
138: *
1.19 ! bertrand 139: * -- LAPACK auxiliary routine --
1.1 bertrand 140: * -- LAPACK is a software package provided by Univ. of Tennessee, --
141: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142: *
143: * .. Scalar Arguments ..
144: LOGICAL WANTQ
145: INTEGER INFO, J1, LDQ, LDT, N, N1, N2
146: * ..
147: * .. Array Arguments ..
148: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
149: * ..
150: *
151: * =====================================================================
152: *
153: * .. Parameters ..
154: DOUBLE PRECISION ZERO, ONE
155: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
156: DOUBLE PRECISION TEN
157: PARAMETER ( TEN = 1.0D+1 )
158: INTEGER LDD, LDX
159: PARAMETER ( LDD = 4, LDX = 2 )
160: * ..
161: * .. Local Scalars ..
162: INTEGER IERR, J2, J3, J4, K, ND
163: DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
164: $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
165: $ WR1, WR2, XNORM
166: * ..
167: * .. Local Arrays ..
168: DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
169: $ X( LDX, 2 )
170: * ..
171: * .. External Functions ..
172: DOUBLE PRECISION DLAMCH, DLANGE
173: EXTERNAL DLAMCH, DLANGE
174: * ..
175: * .. External Subroutines ..
176: EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
177: $ DROT
178: * ..
179: * .. Intrinsic Functions ..
180: INTRINSIC ABS, MAX
181: * ..
182: * .. Executable Statements ..
183: *
184: INFO = 0
185: *
186: * Quick return if possible
187: *
188: IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
189: $ RETURN
190: IF( J1+N1.GT.N )
191: $ RETURN
192: *
193: J2 = J1 + 1
194: J3 = J1 + 2
195: J4 = J1 + 3
196: *
197: IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
198: *
199: * Swap two 1-by-1 blocks.
200: *
201: T11 = T( J1, J1 )
202: T22 = T( J2, J2 )
203: *
204: * Determine the transformation to perform the interchange.
205: *
206: CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
207: *
208: * Apply transformation to the matrix T.
209: *
210: IF( J3.LE.N )
211: $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
212: $ SN )
213: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
214: *
215: T( J1, J1 ) = T22
216: T( J2, J2 ) = T11
217: *
218: IF( WANTQ ) THEN
219: *
220: * Accumulate transformation in the matrix Q.
221: *
222: CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
223: END IF
224: *
225: ELSE
226: *
227: * Swapping involves at least one 2-by-2 block.
228: *
229: * Copy the diagonal block of order N1+N2 to the local array D
230: * and compute its norm.
231: *
232: ND = N1 + N2
233: CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
234: DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
235: *
236: * Compute machine-dependent threshold for test for accepting
237: * swap.
238: *
239: EPS = DLAMCH( 'P' )
240: SMLNUM = DLAMCH( 'S' ) / EPS
241: THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
242: *
243: * Solve T11*X - X*T22 = scale*T12 for X.
244: *
245: CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
246: $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
247: $ LDX, XNORM, IERR )
248: *
249: * Swap the adjacent diagonal blocks.
250: *
251: K = N1 + N1 + N2 - 3
252: GO TO ( 10, 20, 30 )K
253: *
254: 10 CONTINUE
255: *
256: * N1 = 1, N2 = 2: generate elementary reflector H so that:
257: *
258: * ( scale, X11, X12 ) H = ( 0, 0, * )
259: *
260: U( 1 ) = SCALE
261: U( 2 ) = X( 1, 1 )
262: U( 3 ) = X( 1, 2 )
263: CALL DLARFG( 3, U( 3 ), U, 1, TAU )
264: U( 3 ) = ONE
265: T11 = T( J1, J1 )
266: *
267: * Perform swap provisionally on diagonal block in D.
268: *
269: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
270: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
271: *
272: * Test whether to reject swap.
273: *
274: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
275: $ 3 )-T11 ) ).GT.THRESH )GO TO 50
276: *
277: * Accept swap: apply transformation to the entire matrix T.
278: *
279: CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
280: CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
281: *
282: T( J3, J1 ) = ZERO
283: T( J3, J2 ) = ZERO
284: T( J3, J3 ) = T11
285: *
286: IF( WANTQ ) THEN
287: *
288: * Accumulate transformation in the matrix Q.
289: *
290: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
291: END IF
292: GO TO 40
293: *
294: 20 CONTINUE
295: *
296: * N1 = 2, N2 = 1: generate elementary reflector H so that:
297: *
298: * H ( -X11 ) = ( * )
299: * ( -X21 ) = ( 0 )
300: * ( scale ) = ( 0 )
301: *
302: U( 1 ) = -X( 1, 1 )
303: U( 2 ) = -X( 2, 1 )
304: U( 3 ) = SCALE
305: CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
306: U( 1 ) = ONE
307: T33 = T( J3, J3 )
308: *
309: * Perform swap provisionally on diagonal block in D.
310: *
311: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
312: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
313: *
314: * Test whether to reject swap.
315: *
316: IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
317: $ 1 )-T33 ) ).GT.THRESH )GO TO 50
318: *
319: * Accept swap: apply transformation to the entire matrix T.
320: *
321: CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
322: CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
323: *
324: T( J1, J1 ) = T33
325: T( J2, J1 ) = ZERO
326: T( J3, J1 ) = ZERO
327: *
328: IF( WANTQ ) THEN
329: *
330: * Accumulate transformation in the matrix Q.
331: *
332: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
333: END IF
334: GO TO 40
335: *
336: 30 CONTINUE
337: *
338: * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
339: * that:
340: *
341: * H(2) H(1) ( -X11 -X12 ) = ( * * )
342: * ( -X21 -X22 ) ( 0 * )
343: * ( scale 0 ) ( 0 0 )
344: * ( 0 scale ) ( 0 0 )
345: *
346: U1( 1 ) = -X( 1, 1 )
347: U1( 2 ) = -X( 2, 1 )
348: U1( 3 ) = SCALE
349: CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
350: U1( 1 ) = ONE
351: *
352: TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
353: U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
354: U2( 2 ) = -TEMP*U1( 3 )
355: U2( 3 ) = SCALE
356: CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
357: U2( 1 ) = ONE
358: *
359: * Perform swap provisionally on diagonal block in D.
360: *
361: CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
362: CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
363: CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
364: CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
365: *
366: * Test whether to reject swap.
367: *
368: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
369: $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
370: *
371: * Accept swap: apply transformation to the entire matrix T.
372: *
373: CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
374: CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
375: CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
376: CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
377: *
378: T( J3, J1 ) = ZERO
379: T( J3, J2 ) = ZERO
380: T( J4, J1 ) = ZERO
381: T( J4, J2 ) = ZERO
382: *
383: IF( WANTQ ) THEN
384: *
385: * Accumulate transformation in the matrix Q.
386: *
387: CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
388: CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
389: END IF
390: *
391: 40 CONTINUE
392: *
393: IF( N2.EQ.2 ) THEN
394: *
395: * Standardize new 2-by-2 block T11
396: *
397: CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
398: $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
399: CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
400: $ CS, SN )
401: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
402: IF( WANTQ )
403: $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
404: END IF
405: *
406: IF( N1.EQ.2 ) THEN
407: *
408: * Standardize new 2-by-2 block T22
409: *
410: J3 = J1 + N2
411: J4 = J3 + 1
412: CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
413: $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
414: IF( J3+2.LE.N )
415: $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
416: $ LDT, CS, SN )
417: CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
418: IF( WANTQ )
419: $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
420: END IF
421: *
422: END IF
423: RETURN
424: *
425: * Exit with INFO = 1 if swap was rejected.
426: *
427: 50 CONTINUE
428: INFO = 1
429: RETURN
430: *
431: * End of DLAEXC
432: *
433: 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