1: *> \brief \b DTREXC
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTREXC + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
22: * INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER COMPQ
26: * INTEGER IFST, ILST, INFO, LDQ, LDT, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DTREXC reorders the real Schur factorization of a real matrix
39: *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
40: *> moved to row ILST.
41: *>
42: *> The real Schur form T is reordered by an orthogonal similarity
43: *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
44: *> is updated by postmultiplying it with Z.
45: *>
46: *> T must be in Schur canonical form (as returned by DHSEQR), that is,
47: *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
48: *> 2-by-2 diagonal block has its diagonal elements equal and its
49: *> off-diagonal elements of opposite sign.
50: *> \endverbatim
51: *
52: * Arguments:
53: * ==========
54: *
55: *> \param[in] COMPQ
56: *> \verbatim
57: *> COMPQ is CHARACTER*1
58: *> = 'V': update the matrix Q of Schur vectors;
59: *> = 'N': do not update Q.
60: *> \endverbatim
61: *>
62: *> \param[in] N
63: *> \verbatim
64: *> N is INTEGER
65: *> The order of the matrix T. N >= 0.
66: *> If N == 0 arguments ILST and IFST may be any value.
67: *> \endverbatim
68: *>
69: *> \param[in,out] T
70: *> \verbatim
71: *> T is DOUBLE PRECISION array, dimension (LDT,N)
72: *> On entry, the upper quasi-triangular matrix T, in Schur
73: *> Schur canonical form.
74: *> On exit, the reordered upper quasi-triangular matrix, again
75: *> in Schur canonical form.
76: *> \endverbatim
77: *>
78: *> \param[in] LDT
79: *> \verbatim
80: *> LDT is INTEGER
81: *> The leading dimension of the array T. LDT >= max(1,N).
82: *> \endverbatim
83: *>
84: *> \param[in,out] Q
85: *> \verbatim
86: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
87: *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
88: *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
89: *> orthogonal transformation matrix Z which reorders T.
90: *> If COMPQ = 'N', Q is not referenced.
91: *> \endverbatim
92: *>
93: *> \param[in] LDQ
94: *> \verbatim
95: *> LDQ is INTEGER
96: *> The leading dimension of the array Q. LDQ >= 1, and if
97: *> COMPQ = 'V', LDQ >= max(1,N).
98: *> \endverbatim
99: *>
100: *> \param[in,out] IFST
101: *> \verbatim
102: *> IFST is INTEGER
103: *> \endverbatim
104: *>
105: *> \param[in,out] ILST
106: *> \verbatim
107: *> ILST is INTEGER
108: *>
109: *> Specify the reordering of the diagonal blocks of T.
110: *> The block with row index IFST is moved to row ILST, by a
111: *> sequence of transpositions between adjacent blocks.
112: *> On exit, if IFST pointed on entry to the second row of a
113: *> 2-by-2 block, it is changed to point to the first row; ILST
114: *> always points to the first row of the block in its final
115: *> position (which may differ from its input value by +1 or -1).
116: *> 1 <= IFST <= N; 1 <= ILST <= N.
117: *> \endverbatim
118: *>
119: *> \param[out] WORK
120: *> \verbatim
121: *> WORK is DOUBLE PRECISION array, dimension (N)
122: *> \endverbatim
123: *>
124: *> \param[out] INFO
125: *> \verbatim
126: *> INFO is INTEGER
127: *> = 0: successful exit
128: *> < 0: if INFO = -i, the i-th argument had an illegal value
129: *> = 1: two adjacent blocks were too close to swap (the problem
130: *> is very ill-conditioned); T may have been partially
131: *> reordered, and ILST points to the first row of the
132: *> current position of the block being moved.
133: *> \endverbatim
134: *
135: * Authors:
136: * ========
137: *
138: *> \author Univ. of Tennessee
139: *> \author Univ. of California Berkeley
140: *> \author Univ. of Colorado Denver
141: *> \author NAG Ltd.
142: *
143: *> \ingroup doubleOTHERcomputational
144: *
145: * =====================================================================
146: SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
147: $ INFO )
148: *
149: * -- LAPACK computational routine --
150: * -- LAPACK is a software package provided by Univ. of Tennessee, --
151: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152: *
153: * .. Scalar Arguments ..
154: CHARACTER COMPQ
155: INTEGER IFST, ILST, INFO, LDQ, LDT, N
156: * ..
157: * .. Array Arguments ..
158: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
159: * ..
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION ZERO
165: PARAMETER ( ZERO = 0.0D+0 )
166: * ..
167: * .. Local Scalars ..
168: LOGICAL WANTQ
169: INTEGER HERE, NBF, NBL, NBNEXT
170: * ..
171: * .. External Functions ..
172: LOGICAL LSAME
173: EXTERNAL LSAME
174: * ..
175: * .. External Subroutines ..
176: EXTERNAL DLAEXC, XERBLA
177: * ..
178: * .. Intrinsic Functions ..
179: INTRINSIC MAX
180: * ..
181: * .. Executable Statements ..
182: *
183: * Decode and test the input arguments.
184: *
185: INFO = 0
186: WANTQ = LSAME( COMPQ, 'V' )
187: IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
188: INFO = -1
189: ELSE IF( N.LT.0 ) THEN
190: INFO = -2
191: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
192: INFO = -4
193: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
194: INFO = -6
195: ELSE IF(( IFST.LT.1 .OR. IFST.GT.N ).AND.( N.GT.0 )) THEN
196: INFO = -7
197: ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) THEN
198: INFO = -8
199: END IF
200: IF( INFO.NE.0 ) THEN
201: CALL XERBLA( 'DTREXC', -INFO )
202: RETURN
203: END IF
204: *
205: * Quick return if possible
206: *
207: IF( N.LE.1 )
208: $ RETURN
209: *
210: * Determine the first row of specified block
211: * and find out it is 1 by 1 or 2 by 2.
212: *
213: IF( IFST.GT.1 ) THEN
214: IF( T( IFST, IFST-1 ).NE.ZERO )
215: $ IFST = IFST - 1
216: END IF
217: NBF = 1
218: IF( IFST.LT.N ) THEN
219: IF( T( IFST+1, IFST ).NE.ZERO )
220: $ NBF = 2
221: END IF
222: *
223: * Determine the first row of the final block
224: * and find out it is 1 by 1 or 2 by 2.
225: *
226: IF( ILST.GT.1 ) THEN
227: IF( T( ILST, ILST-1 ).NE.ZERO )
228: $ ILST = ILST - 1
229: END IF
230: NBL = 1
231: IF( ILST.LT.N ) THEN
232: IF( T( ILST+1, ILST ).NE.ZERO )
233: $ NBL = 2
234: END IF
235: *
236: IF( IFST.EQ.ILST )
237: $ RETURN
238: *
239: IF( IFST.LT.ILST ) THEN
240: *
241: * Update ILST
242: *
243: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
244: $ ILST = ILST - 1
245: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
246: $ ILST = ILST + 1
247: *
248: HERE = IFST
249: *
250: 10 CONTINUE
251: *
252: * Swap block with next one below
253: *
254: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
255: *
256: * Current block either 1 by 1 or 2 by 2
257: *
258: NBNEXT = 1
259: IF( HERE+NBF+1.LE.N ) THEN
260: IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
261: $ NBNEXT = 2
262: END IF
263: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
264: $ WORK, INFO )
265: IF( INFO.NE.0 ) THEN
266: ILST = HERE
267: RETURN
268: END IF
269: HERE = HERE + NBNEXT
270: *
271: * Test if 2 by 2 block breaks into two 1 by 1 blocks
272: *
273: IF( NBF.EQ.2 ) THEN
274: IF( T( HERE+1, HERE ).EQ.ZERO )
275: $ NBF = 3
276: END IF
277: *
278: ELSE
279: *
280: * Current block consists of two 1 by 1 blocks each of which
281: * must be swapped individually
282: *
283: NBNEXT = 1
284: IF( HERE+3.LE.N ) THEN
285: IF( T( HERE+3, HERE+2 ).NE.ZERO )
286: $ NBNEXT = 2
287: END IF
288: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
289: $ WORK, INFO )
290: IF( INFO.NE.0 ) THEN
291: ILST = HERE
292: RETURN
293: END IF
294: IF( NBNEXT.EQ.1 ) THEN
295: *
296: * Swap two 1 by 1 blocks, no problems possible
297: *
298: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
299: $ WORK, INFO )
300: HERE = HERE + 1
301: ELSE
302: *
303: * Recompute NBNEXT in case 2 by 2 split
304: *
305: IF( T( HERE+2, HERE+1 ).EQ.ZERO )
306: $ NBNEXT = 1
307: IF( NBNEXT.EQ.2 ) THEN
308: *
309: * 2 by 2 Block did not split
310: *
311: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
312: $ NBNEXT, WORK, INFO )
313: IF( INFO.NE.0 ) THEN
314: ILST = HERE
315: RETURN
316: END IF
317: HERE = HERE + 2
318: ELSE
319: *
320: * 2 by 2 Block did split
321: *
322: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
323: $ WORK, INFO )
324: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
325: $ WORK, INFO )
326: HERE = HERE + 2
327: END IF
328: END IF
329: END IF
330: IF( HERE.LT.ILST )
331: $ GO TO 10
332: *
333: ELSE
334: *
335: HERE = IFST
336: 20 CONTINUE
337: *
338: * Swap block with next one above
339: *
340: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
341: *
342: * Current block either 1 by 1 or 2 by 2
343: *
344: NBNEXT = 1
345: IF( HERE.GE.3 ) THEN
346: IF( T( HERE-1, HERE-2 ).NE.ZERO )
347: $ NBNEXT = 2
348: END IF
349: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
350: $ NBF, WORK, INFO )
351: IF( INFO.NE.0 ) THEN
352: ILST = HERE
353: RETURN
354: END IF
355: HERE = HERE - NBNEXT
356: *
357: * Test if 2 by 2 block breaks into two 1 by 1 blocks
358: *
359: IF( NBF.EQ.2 ) THEN
360: IF( T( HERE+1, HERE ).EQ.ZERO )
361: $ NBF = 3
362: END IF
363: *
364: ELSE
365: *
366: * Current block consists of two 1 by 1 blocks each of which
367: * must be swapped individually
368: *
369: NBNEXT = 1
370: IF( HERE.GE.3 ) THEN
371: IF( T( HERE-1, HERE-2 ).NE.ZERO )
372: $ NBNEXT = 2
373: END IF
374: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
375: $ 1, WORK, INFO )
376: IF( INFO.NE.0 ) THEN
377: ILST = HERE
378: RETURN
379: END IF
380: IF( NBNEXT.EQ.1 ) THEN
381: *
382: * Swap two 1 by 1 blocks, no problems possible
383: *
384: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
385: $ WORK, INFO )
386: HERE = HERE - 1
387: ELSE
388: *
389: * Recompute NBNEXT in case 2 by 2 split
390: *
391: IF( T( HERE, HERE-1 ).EQ.ZERO )
392: $ NBNEXT = 1
393: IF( NBNEXT.EQ.2 ) THEN
394: *
395: * 2 by 2 Block did not split
396: *
397: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
398: $ WORK, INFO )
399: IF( INFO.NE.0 ) THEN
400: ILST = HERE
401: RETURN
402: END IF
403: HERE = HERE - 2
404: ELSE
405: *
406: * 2 by 2 Block did split
407: *
408: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
409: $ WORK, INFO )
410: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
411: $ WORK, INFO )
412: HERE = HERE - 2
413: END IF
414: END IF
415: END IF
416: IF( HERE.GT.ILST )
417: $ GO TO 20
418: END IF
419: ILST = HERE
420: *
421: RETURN
422: *
423: * End of DTREXC
424: *
425: 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