1: SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.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: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER COMPQ
11: INTEGER IFST, ILST, INFO, LDQ, LDT, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DTREXC reorders the real Schur factorization of a real matrix
21: * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
22: * moved to row ILST.
23: *
24: * The real Schur form T is reordered by an orthogonal similarity
25: * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
26: * is updated by postmultiplying it with Z.
27: *
28: * T must be in Schur canonical form (as returned by DHSEQR), that is,
29: * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
30: * 2-by-2 diagonal block has its diagonal elements equal and its
31: * off-diagonal elements of opposite sign.
32: *
33: * Arguments
34: * =========
35: *
36: * COMPQ (input) CHARACTER*1
37: * = 'V': update the matrix Q of Schur vectors;
38: * = 'N': do not update Q.
39: *
40: * N (input) INTEGER
41: * The order of the matrix T. N >= 0.
42: *
43: * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
44: * On entry, the upper quasi-triangular matrix T, in Schur
45: * Schur canonical form.
46: * On exit, the reordered upper quasi-triangular matrix, again
47: * in Schur canonical form.
48: *
49: * LDT (input) INTEGER
50: * The leading dimension of the array T. LDT >= max(1,N).
51: *
52: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
53: * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
54: * On exit, if COMPQ = 'V', Q has been postmultiplied by the
55: * orthogonal transformation matrix Z which reorders T.
56: * If COMPQ = 'N', Q is not referenced.
57: *
58: * LDQ (input) INTEGER
59: * The leading dimension of the array Q. LDQ >= max(1,N).
60: *
61: * IFST (input/output) INTEGER
62: * ILST (input/output) INTEGER
63: * Specify the reordering of the diagonal blocks of T.
64: * The block with row index IFST is moved to row ILST, by a
65: * sequence of transpositions between adjacent blocks.
66: * On exit, if IFST pointed on entry to the second row of a
67: * 2-by-2 block, it is changed to point to the first row; ILST
68: * always points to the first row of the block in its final
69: * position (which may differ from its input value by +1 or -1).
70: * 1 <= IFST <= N; 1 <= ILST <= N.
71: *
72: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
73: *
74: * INFO (output) INTEGER
75: * = 0: successful exit
76: * < 0: if INFO = -i, the i-th argument had an illegal value
77: * = 1: two adjacent blocks were too close to swap (the problem
78: * is very ill-conditioned); T may have been partially
79: * reordered, and ILST points to the first row of the
80: * current position of the block being moved.
81: *
82: * =====================================================================
83: *
84: * .. Parameters ..
85: DOUBLE PRECISION ZERO
86: PARAMETER ( ZERO = 0.0D+0 )
87: * ..
88: * .. Local Scalars ..
89: LOGICAL WANTQ
90: INTEGER HERE, NBF, NBL, NBNEXT
91: * ..
92: * .. External Functions ..
93: LOGICAL LSAME
94: EXTERNAL LSAME
95: * ..
96: * .. External Subroutines ..
97: EXTERNAL DLAEXC, XERBLA
98: * ..
99: * .. Intrinsic Functions ..
100: INTRINSIC MAX
101: * ..
102: * .. Executable Statements ..
103: *
104: * Decode and test the input arguments.
105: *
106: INFO = 0
107: WANTQ = LSAME( COMPQ, 'V' )
108: IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
109: INFO = -1
110: ELSE IF( N.LT.0 ) THEN
111: INFO = -2
112: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
113: INFO = -4
114: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
115: INFO = -6
116: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
117: INFO = -7
118: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
119: INFO = -8
120: END IF
121: IF( INFO.NE.0 ) THEN
122: CALL XERBLA( 'DTREXC', -INFO )
123: RETURN
124: END IF
125: *
126: * Quick return if possible
127: *
128: IF( N.LE.1 )
129: $ RETURN
130: *
131: * Determine the first row of specified block
132: * and find out it is 1 by 1 or 2 by 2.
133: *
134: IF( IFST.GT.1 ) THEN
135: IF( T( IFST, IFST-1 ).NE.ZERO )
136: $ IFST = IFST - 1
137: END IF
138: NBF = 1
139: IF( IFST.LT.N ) THEN
140: IF( T( IFST+1, IFST ).NE.ZERO )
141: $ NBF = 2
142: END IF
143: *
144: * Determine the first row of the final block
145: * and find out it is 1 by 1 or 2 by 2.
146: *
147: IF( ILST.GT.1 ) THEN
148: IF( T( ILST, ILST-1 ).NE.ZERO )
149: $ ILST = ILST - 1
150: END IF
151: NBL = 1
152: IF( ILST.LT.N ) THEN
153: IF( T( ILST+1, ILST ).NE.ZERO )
154: $ NBL = 2
155: END IF
156: *
157: IF( IFST.EQ.ILST )
158: $ RETURN
159: *
160: IF( IFST.LT.ILST ) THEN
161: *
162: * Update ILST
163: *
164: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
165: $ ILST = ILST - 1
166: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
167: $ ILST = ILST + 1
168: *
169: HERE = IFST
170: *
171: 10 CONTINUE
172: *
173: * Swap block with next one below
174: *
175: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
176: *
177: * Current block either 1 by 1 or 2 by 2
178: *
179: NBNEXT = 1
180: IF( HERE+NBF+1.LE.N ) THEN
181: IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
182: $ NBNEXT = 2
183: END IF
184: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
185: $ WORK, INFO )
186: IF( INFO.NE.0 ) THEN
187: ILST = HERE
188: RETURN
189: END IF
190: HERE = HERE + NBNEXT
191: *
192: * Test if 2 by 2 block breaks into two 1 by 1 blocks
193: *
194: IF( NBF.EQ.2 ) THEN
195: IF( T( HERE+1, HERE ).EQ.ZERO )
196: $ NBF = 3
197: END IF
198: *
199: ELSE
200: *
201: * Current block consists of two 1 by 1 blocks each of which
202: * must be swapped individually
203: *
204: NBNEXT = 1
205: IF( HERE+3.LE.N ) THEN
206: IF( T( HERE+3, HERE+2 ).NE.ZERO )
207: $ NBNEXT = 2
208: END IF
209: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
210: $ WORK, INFO )
211: IF( INFO.NE.0 ) THEN
212: ILST = HERE
213: RETURN
214: END IF
215: IF( NBNEXT.EQ.1 ) THEN
216: *
217: * Swap two 1 by 1 blocks, no problems possible
218: *
219: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
220: $ WORK, INFO )
221: HERE = HERE + 1
222: ELSE
223: *
224: * Recompute NBNEXT in case 2 by 2 split
225: *
226: IF( T( HERE+2, HERE+1 ).EQ.ZERO )
227: $ NBNEXT = 1
228: IF( NBNEXT.EQ.2 ) THEN
229: *
230: * 2 by 2 Block did not split
231: *
232: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
233: $ NBNEXT, WORK, INFO )
234: IF( INFO.NE.0 ) THEN
235: ILST = HERE
236: RETURN
237: END IF
238: HERE = HERE + 2
239: ELSE
240: *
241: * 2 by 2 Block did split
242: *
243: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
244: $ WORK, INFO )
245: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
246: $ WORK, INFO )
247: HERE = HERE + 2
248: END IF
249: END IF
250: END IF
251: IF( HERE.LT.ILST )
252: $ GO TO 10
253: *
254: ELSE
255: *
256: HERE = IFST
257: 20 CONTINUE
258: *
259: * Swap block with next one above
260: *
261: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
262: *
263: * Current block either 1 by 1 or 2 by 2
264: *
265: NBNEXT = 1
266: IF( HERE.GE.3 ) THEN
267: IF( T( HERE-1, HERE-2 ).NE.ZERO )
268: $ NBNEXT = 2
269: END IF
270: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
271: $ NBF, WORK, INFO )
272: IF( INFO.NE.0 ) THEN
273: ILST = HERE
274: RETURN
275: END IF
276: HERE = HERE - NBNEXT
277: *
278: * Test if 2 by 2 block breaks into two 1 by 1 blocks
279: *
280: IF( NBF.EQ.2 ) THEN
281: IF( T( HERE+1, HERE ).EQ.ZERO )
282: $ NBF = 3
283: END IF
284: *
285: ELSE
286: *
287: * Current block consists of two 1 by 1 blocks each of which
288: * must be swapped individually
289: *
290: NBNEXT = 1
291: IF( HERE.GE.3 ) THEN
292: IF( T( HERE-1, HERE-2 ).NE.ZERO )
293: $ NBNEXT = 2
294: END IF
295: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
296: $ 1, WORK, INFO )
297: IF( INFO.NE.0 ) THEN
298: ILST = HERE
299: RETURN
300: END IF
301: IF( NBNEXT.EQ.1 ) THEN
302: *
303: * Swap two 1 by 1 blocks, no problems possible
304: *
305: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
306: $ WORK, INFO )
307: HERE = HERE - 1
308: ELSE
309: *
310: * Recompute NBNEXT in case 2 by 2 split
311: *
312: IF( T( HERE, HERE-1 ).EQ.ZERO )
313: $ NBNEXT = 1
314: IF( NBNEXT.EQ.2 ) THEN
315: *
316: * 2 by 2 Block did not split
317: *
318: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
319: $ WORK, INFO )
320: IF( INFO.NE.0 ) THEN
321: ILST = HERE
322: RETURN
323: END IF
324: HERE = HERE - 2
325: ELSE
326: *
327: * 2 by 2 Block did split
328: *
329: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
330: $ WORK, INFO )
331: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
332: $ WORK, INFO )
333: HERE = HERE - 2
334: END IF
335: END IF
336: END IF
337: IF( HERE.GT.ILST )
338: $ GO TO 20
339: END IF
340: ILST = HERE
341: *
342: RETURN
343: *
344: * End of DTREXC
345: *
346: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack