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