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: *> \date December 2016
144: *
145: *> \ingroup doubleOTHERcomputational
146: *
147: * =====================================================================
148: SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
149: $ INFO )
150: *
151: * -- LAPACK computational routine (version 3.7.0) --
152: * -- LAPACK is a software package provided by Univ. of Tennessee, --
153: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154: * December 2016
155: *
156: * .. Scalar Arguments ..
157: CHARACTER COMPQ
158: INTEGER IFST, ILST, INFO, LDQ, LDT, N
159: * ..
160: * .. Array Arguments ..
161: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
162: * ..
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION ZERO
168: PARAMETER ( ZERO = 0.0D+0 )
169: * ..
170: * .. Local Scalars ..
171: LOGICAL WANTQ
172: INTEGER HERE, NBF, NBL, NBNEXT
173: * ..
174: * .. External Functions ..
175: LOGICAL LSAME
176: EXTERNAL LSAME
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL DLAEXC, XERBLA
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC MAX
183: * ..
184: * .. Executable Statements ..
185: *
186: * Decode and test the input arguments.
187: *
188: INFO = 0
189: WANTQ = LSAME( COMPQ, 'V' )
190: IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
191: INFO = -1
192: ELSE IF( N.LT.0 ) THEN
193: INFO = -2
194: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
195: INFO = -4
196: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
197: INFO = -6
198: ELSE IF(( IFST.LT.1 .OR. IFST.GT.N ).AND.( N.GT.0 )) THEN
199: INFO = -7
200: ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) THEN
201: INFO = -8
202: END IF
203: IF( INFO.NE.0 ) THEN
204: CALL XERBLA( 'DTREXC', -INFO )
205: RETURN
206: END IF
207: *
208: * Quick return if possible
209: *
210: IF( N.LE.1 )
211: $ RETURN
212: *
213: * Determine the first row of specified block
214: * and find out it is 1 by 1 or 2 by 2.
215: *
216: IF( IFST.GT.1 ) THEN
217: IF( T( IFST, IFST-1 ).NE.ZERO )
218: $ IFST = IFST - 1
219: END IF
220: NBF = 1
221: IF( IFST.LT.N ) THEN
222: IF( T( IFST+1, IFST ).NE.ZERO )
223: $ NBF = 2
224: END IF
225: *
226: * Determine the first row of the final block
227: * and find out it is 1 by 1 or 2 by 2.
228: *
229: IF( ILST.GT.1 ) THEN
230: IF( T( ILST, ILST-1 ).NE.ZERO )
231: $ ILST = ILST - 1
232: END IF
233: NBL = 1
234: IF( ILST.LT.N ) THEN
235: IF( T( ILST+1, ILST ).NE.ZERO )
236: $ NBL = 2
237: END IF
238: *
239: IF( IFST.EQ.ILST )
240: $ RETURN
241: *
242: IF( IFST.LT.ILST ) THEN
243: *
244: * Update ILST
245: *
246: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
247: $ ILST = ILST - 1
248: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
249: $ ILST = ILST + 1
250: *
251: HERE = IFST
252: *
253: 10 CONTINUE
254: *
255: * Swap block with next one below
256: *
257: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
258: *
259: * Current block either 1 by 1 or 2 by 2
260: *
261: NBNEXT = 1
262: IF( HERE+NBF+1.LE.N ) THEN
263: IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
264: $ NBNEXT = 2
265: END IF
266: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
267: $ WORK, INFO )
268: IF( INFO.NE.0 ) THEN
269: ILST = HERE
270: RETURN
271: END IF
272: HERE = HERE + NBNEXT
273: *
274: * Test if 2 by 2 block breaks into two 1 by 1 blocks
275: *
276: IF( NBF.EQ.2 ) THEN
277: IF( T( HERE+1, HERE ).EQ.ZERO )
278: $ NBF = 3
279: END IF
280: *
281: ELSE
282: *
283: * Current block consists of two 1 by 1 blocks each of which
284: * must be swapped individually
285: *
286: NBNEXT = 1
287: IF( HERE+3.LE.N ) THEN
288: IF( T( HERE+3, HERE+2 ).NE.ZERO )
289: $ NBNEXT = 2
290: END IF
291: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
292: $ WORK, INFO )
293: IF( INFO.NE.0 ) THEN
294: ILST = HERE
295: RETURN
296: END IF
297: IF( NBNEXT.EQ.1 ) THEN
298: *
299: * Swap two 1 by 1 blocks, no problems possible
300: *
301: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
302: $ WORK, INFO )
303: HERE = HERE + 1
304: ELSE
305: *
306: * Recompute NBNEXT in case 2 by 2 split
307: *
308: IF( T( HERE+2, HERE+1 ).EQ.ZERO )
309: $ NBNEXT = 1
310: IF( NBNEXT.EQ.2 ) THEN
311: *
312: * 2 by 2 Block did not split
313: *
314: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
315: $ NBNEXT, WORK, INFO )
316: IF( INFO.NE.0 ) THEN
317: ILST = HERE
318: RETURN
319: END IF
320: HERE = HERE + 2
321: ELSE
322: *
323: * 2 by 2 Block did split
324: *
325: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
326: $ WORK, INFO )
327: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
328: $ WORK, INFO )
329: HERE = HERE + 2
330: END IF
331: END IF
332: END IF
333: IF( HERE.LT.ILST )
334: $ GO TO 10
335: *
336: ELSE
337: *
338: HERE = IFST
339: 20 CONTINUE
340: *
341: * Swap block with next one above
342: *
343: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
344: *
345: * Current block either 1 by 1 or 2 by 2
346: *
347: NBNEXT = 1
348: IF( HERE.GE.3 ) THEN
349: IF( T( HERE-1, HERE-2 ).NE.ZERO )
350: $ NBNEXT = 2
351: END IF
352: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
353: $ NBF, WORK, INFO )
354: IF( INFO.NE.0 ) THEN
355: ILST = HERE
356: RETURN
357: END IF
358: HERE = HERE - NBNEXT
359: *
360: * Test if 2 by 2 block breaks into two 1 by 1 blocks
361: *
362: IF( NBF.EQ.2 ) THEN
363: IF( T( HERE+1, HERE ).EQ.ZERO )
364: $ NBF = 3
365: END IF
366: *
367: ELSE
368: *
369: * Current block consists of two 1 by 1 blocks each of which
370: * must be swapped individually
371: *
372: NBNEXT = 1
373: IF( HERE.GE.3 ) THEN
374: IF( T( HERE-1, HERE-2 ).NE.ZERO )
375: $ NBNEXT = 2
376: END IF
377: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
378: $ 1, WORK, INFO )
379: IF( INFO.NE.0 ) THEN
380: ILST = HERE
381: RETURN
382: END IF
383: IF( NBNEXT.EQ.1 ) THEN
384: *
385: * Swap two 1 by 1 blocks, no problems possible
386: *
387: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
388: $ WORK, INFO )
389: HERE = HERE - 1
390: ELSE
391: *
392: * Recompute NBNEXT in case 2 by 2 split
393: *
394: IF( T( HERE, HERE-1 ).EQ.ZERO )
395: $ NBNEXT = 1
396: IF( NBNEXT.EQ.2 ) THEN
397: *
398: * 2 by 2 Block did not split
399: *
400: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
401: $ WORK, INFO )
402: IF( INFO.NE.0 ) THEN
403: ILST = HERE
404: RETURN
405: END IF
406: HERE = HERE - 2
407: ELSE
408: *
409: * 2 by 2 Block did split
410: *
411: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
412: $ WORK, INFO )
413: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
414: $ WORK, INFO )
415: HERE = HERE - 2
416: END IF
417: END IF
418: END IF
419: IF( HERE.GT.ILST )
420: $ GO TO 20
421: END IF
422: ILST = HERE
423: *
424: RETURN
425: *
426: * End of DTREXC
427: *
428: 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