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