Annotation of rpl/lapack/lapack/dtrexc.f, revision 1.8
1.8 ! bertrand 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: * =====================================================================
1.1 bertrand 146: SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
147: $ INFO )
148: *
1.8 ! bertrand 149: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 150: * -- LAPACK is a software package provided by Univ. of Tennessee, --
151: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 152: * November 2011
1.1 bertrand 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