1: *> \brief \b DTGEXC
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTGEXC + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgexc.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgexc.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgexc.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
22: * LDZ, IFST, ILST, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL WANTQ, WANTZ
26: * INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30: * $ WORK( * ), Z( LDZ, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DTGEXC reorders the generalized real Schur decomposition of a real
40: *> matrix pair (A,B) using an orthogonal equivalence transformation
41: *>
42: *> (A, B) = Q * (A, B) * Z**T,
43: *>
44: *> so that the diagonal block of (A, B) with row index IFST is moved
45: *> to row ILST.
46: *>
47: *> (A, B) must be in generalized real Schur canonical form (as returned
48: *> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
49: *> diagonal blocks. B is upper triangular.
50: *>
51: *> Optionally, the matrices Q and Z of generalized Schur vectors are
52: *> updated.
53: *>
54: *> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
55: *> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
56: *>
57: *> \endverbatim
58: *
59: * Arguments:
60: * ==========
61: *
62: *> \param[in] WANTQ
63: *> \verbatim
64: *> WANTQ is LOGICAL
65: *> .TRUE. : update the left transformation matrix Q;
66: *> .FALSE.: do not update Q.
67: *> \endverbatim
68: *>
69: *> \param[in] WANTZ
70: *> \verbatim
71: *> WANTZ is LOGICAL
72: *> .TRUE. : update the right transformation matrix Z;
73: *> .FALSE.: do not update Z.
74: *> \endverbatim
75: *>
76: *> \param[in] N
77: *> \verbatim
78: *> N is INTEGER
79: *> The order of the matrices A and B. N >= 0.
80: *> \endverbatim
81: *>
82: *> \param[in,out] A
83: *> \verbatim
84: *> A is DOUBLE PRECISION array, dimension (LDA,N)
85: *> On entry, the matrix A in generalized real Schur canonical
86: *> form.
87: *> On exit, the updated matrix A, again in generalized
88: *> real Schur canonical form.
89: *> \endverbatim
90: *>
91: *> \param[in] LDA
92: *> \verbatim
93: *> LDA is INTEGER
94: *> The leading dimension of the array A. LDA >= max(1,N).
95: *> \endverbatim
96: *>
97: *> \param[in,out] B
98: *> \verbatim
99: *> B is DOUBLE PRECISION array, dimension (LDB,N)
100: *> On entry, the matrix B in generalized real Schur canonical
101: *> form (A,B).
102: *> On exit, the updated matrix B, again in generalized
103: *> real Schur canonical form (A,B).
104: *> \endverbatim
105: *>
106: *> \param[in] LDB
107: *> \verbatim
108: *> LDB is INTEGER
109: *> The leading dimension of the array B. LDB >= max(1,N).
110: *> \endverbatim
111: *>
112: *> \param[in,out] Q
113: *> \verbatim
114: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
115: *> On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
116: *> On exit, the updated matrix Q.
117: *> If WANTQ = .FALSE., Q is not referenced.
118: *> \endverbatim
119: *>
120: *> \param[in] LDQ
121: *> \verbatim
122: *> LDQ is INTEGER
123: *> The leading dimension of the array Q. LDQ >= 1.
124: *> If WANTQ = .TRUE., LDQ >= N.
125: *> \endverbatim
126: *>
127: *> \param[in,out] Z
128: *> \verbatim
129: *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
130: *> On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
131: *> On exit, the updated matrix Z.
132: *> If WANTZ = .FALSE., Z is not referenced.
133: *> \endverbatim
134: *>
135: *> \param[in] LDZ
136: *> \verbatim
137: *> LDZ is INTEGER
138: *> The leading dimension of the array Z. LDZ >= 1.
139: *> If WANTZ = .TRUE., LDZ >= N.
140: *> \endverbatim
141: *>
142: *> \param[in,out] IFST
143: *> \verbatim
144: *> IFST is INTEGER
145: *> \endverbatim
146: *>
147: *> \param[in,out] ILST
148: *> \verbatim
149: *> ILST is INTEGER
150: *> Specify the reordering of the diagonal blocks of (A, B).
151: *> The block with row index IFST is moved to row ILST, by a
152: *> sequence of swapping between adjacent blocks.
153: *> On exit, if IFST pointed on entry to the second row of
154: *> a 2-by-2 block, it is changed to point to the first row;
155: *> ILST always points to the first row of the block in its
156: *> final position (which may differ from its input value by
157: *> +1 or -1). 1 <= IFST, ILST <= N.
158: *> \endverbatim
159: *>
160: *> \param[out] WORK
161: *> \verbatim
162: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
163: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164: *> \endverbatim
165: *>
166: *> \param[in] LWORK
167: *> \verbatim
168: *> LWORK is INTEGER
169: *> The dimension of the array WORK.
170: *> LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
171: *>
172: *> If LWORK = -1, then a workspace query is assumed; the routine
173: *> only calculates the optimal size of the WORK array, returns
174: *> this value as the first entry of the WORK array, and no error
175: *> message related to LWORK is issued by XERBLA.
176: *> \endverbatim
177: *>
178: *> \param[out] INFO
179: *> \verbatim
180: *> INFO is INTEGER
181: *> =0: successful exit.
182: *> <0: if INFO = -i, the i-th argument had an illegal value.
183: *> =1: The transformed matrix pair (A, B) would be too far
184: *> from generalized Schur form; the problem is ill-
185: *> conditioned. (A, B) may have been partially reordered,
186: *> and ILST points to the first row of the current
187: *> position of the block being moved.
188: *> \endverbatim
189: *
190: * Authors:
191: * ========
192: *
193: *> \author Univ. of Tennessee
194: *> \author Univ. of California Berkeley
195: *> \author Univ. of Colorado Denver
196: *> \author NAG Ltd.
197: *
198: *> \ingroup doubleGEcomputational
199: *
200: *> \par Contributors:
201: * ==================
202: *>
203: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
204: *> Umea University, S-901 87 Umea, Sweden.
205: *
206: *> \par References:
207: * ================
208: *>
209: *> \verbatim
210: *>
211: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
212: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
213: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
214: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
215: *> \endverbatim
216: *>
217: * =====================================================================
218: SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
219: $ LDZ, IFST, ILST, WORK, LWORK, INFO )
220: *
221: * -- LAPACK computational routine --
222: * -- LAPACK is a software package provided by Univ. of Tennessee, --
223: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224: *
225: * .. Scalar Arguments ..
226: LOGICAL WANTQ, WANTZ
227: INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
228: * ..
229: * .. Array Arguments ..
230: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
231: $ WORK( * ), Z( LDZ, * )
232: * ..
233: *
234: * =====================================================================
235: *
236: * .. Parameters ..
237: DOUBLE PRECISION ZERO
238: PARAMETER ( ZERO = 0.0D+0 )
239: * ..
240: * .. Local Scalars ..
241: LOGICAL LQUERY
242: INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
243: * ..
244: * .. External Subroutines ..
245: EXTERNAL DTGEX2, XERBLA
246: * ..
247: * .. Intrinsic Functions ..
248: INTRINSIC MAX
249: * ..
250: * .. Executable Statements ..
251: *
252: * Decode and test input arguments.
253: *
254: INFO = 0
255: LQUERY = ( LWORK.EQ.-1 )
256: IF( N.LT.0 ) THEN
257: INFO = -3
258: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
259: INFO = -5
260: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
261: INFO = -7
262: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
263: INFO = -9
264: ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
265: INFO = -11
266: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
267: INFO = -12
268: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
269: INFO = -13
270: END IF
271: *
272: IF( INFO.EQ.0 ) THEN
273: IF( N.LE.1 ) THEN
274: LWMIN = 1
275: ELSE
276: LWMIN = 4*N + 16
277: END IF
278: WORK(1) = LWMIN
279: *
280: IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
281: INFO = -15
282: END IF
283: END IF
284: *
285: IF( INFO.NE.0 ) THEN
286: CALL XERBLA( 'DTGEXC', -INFO )
287: RETURN
288: ELSE IF( LQUERY ) THEN
289: RETURN
290: END IF
291: *
292: * Quick return if possible
293: *
294: IF( N.LE.1 )
295: $ RETURN
296: *
297: * Determine the first row of the specified block and find out
298: * if it is 1-by-1 or 2-by-2.
299: *
300: IF( IFST.GT.1 ) THEN
301: IF( A( IFST, IFST-1 ).NE.ZERO )
302: $ IFST = IFST - 1
303: END IF
304: NBF = 1
305: IF( IFST.LT.N ) THEN
306: IF( A( IFST+1, IFST ).NE.ZERO )
307: $ NBF = 2
308: END IF
309: *
310: * Determine the first row of the final block
311: * and find out if it is 1-by-1 or 2-by-2.
312: *
313: IF( ILST.GT.1 ) THEN
314: IF( A( ILST, ILST-1 ).NE.ZERO )
315: $ ILST = ILST - 1
316: END IF
317: NBL = 1
318: IF( ILST.LT.N ) THEN
319: IF( A( ILST+1, ILST ).NE.ZERO )
320: $ NBL = 2
321: END IF
322: IF( IFST.EQ.ILST )
323: $ RETURN
324: *
325: IF( IFST.LT.ILST ) THEN
326: *
327: * Update ILST.
328: *
329: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
330: $ ILST = ILST - 1
331: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
332: $ ILST = ILST + 1
333: *
334: HERE = IFST
335: *
336: 10 CONTINUE
337: *
338: * Swap with next one below.
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+NBF+1.LE.N ) THEN
346: IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
347: $ NBNEXT = 2
348: END IF
349: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
350: $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, 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( A( 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+3.LE.N ) THEN
371: IF( A( HERE+3, HERE+2 ).NE.ZERO )
372: $ NBNEXT = 2
373: END IF
374: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
375: $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, 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.
383: *
384: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
385: $ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
386: IF( INFO.NE.0 ) THEN
387: ILST = HERE
388: RETURN
389: END IF
390: HERE = HERE + 1
391: *
392: ELSE
393: *
394: * Recompute NBNEXT in case of 2-by-2 split.
395: *
396: IF( A( HERE+2, HERE+1 ).EQ.ZERO )
397: $ NBNEXT = 1
398: IF( NBNEXT.EQ.2 ) THEN
399: *
400: * 2-by-2 block did not split.
401: *
402: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
403: $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
404: $ INFO )
405: IF( INFO.NE.0 ) THEN
406: ILST = HERE
407: RETURN
408: END IF
409: HERE = HERE + 2
410: ELSE
411: *
412: * 2-by-2 block did split.
413: *
414: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
415: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
416: IF( INFO.NE.0 ) THEN
417: ILST = HERE
418: RETURN
419: END IF
420: HERE = HERE + 1
421: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
422: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
423: IF( INFO.NE.0 ) THEN
424: ILST = HERE
425: RETURN
426: END IF
427: HERE = HERE + 1
428: END IF
429: *
430: END IF
431: END IF
432: IF( HERE.LT.ILST )
433: $ GO TO 10
434: ELSE
435: HERE = IFST
436: *
437: 20 CONTINUE
438: *
439: * Swap with next one below.
440: *
441: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
442: *
443: * Current block either 1-by-1 or 2-by-2.
444: *
445: NBNEXT = 1
446: IF( HERE.GE.3 ) THEN
447: IF( A( HERE-1, HERE-2 ).NE.ZERO )
448: $ NBNEXT = 2
449: END IF
450: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
451: $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
452: $ INFO )
453: IF( INFO.NE.0 ) THEN
454: ILST = HERE
455: RETURN
456: END IF
457: HERE = HERE - NBNEXT
458: *
459: * Test if 2-by-2 block breaks into two 1-by-1 blocks.
460: *
461: IF( NBF.EQ.2 ) THEN
462: IF( A( HERE+1, HERE ).EQ.ZERO )
463: $ NBF = 3
464: END IF
465: *
466: ELSE
467: *
468: * Current block consists of two 1-by-1 blocks, each of which
469: * must be swapped individually.
470: *
471: NBNEXT = 1
472: IF( HERE.GE.3 ) THEN
473: IF( A( HERE-1, HERE-2 ).NE.ZERO )
474: $ NBNEXT = 2
475: END IF
476: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
477: $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
478: $ INFO )
479: IF( INFO.NE.0 ) THEN
480: ILST = HERE
481: RETURN
482: END IF
483: IF( NBNEXT.EQ.1 ) THEN
484: *
485: * Swap two 1-by-1 blocks.
486: *
487: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
488: $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
489: IF( INFO.NE.0 ) THEN
490: ILST = HERE
491: RETURN
492: END IF
493: HERE = HERE - 1
494: ELSE
495: *
496: * Recompute NBNEXT in case of 2-by-2 split.
497: *
498: IF( A( HERE, HERE-1 ).EQ.ZERO )
499: $ NBNEXT = 1
500: IF( NBNEXT.EQ.2 ) THEN
501: *
502: * 2-by-2 block did not split.
503: *
504: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
505: $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
506: IF( INFO.NE.0 ) THEN
507: ILST = HERE
508: RETURN
509: END IF
510: HERE = HERE - 2
511: ELSE
512: *
513: * 2-by-2 block did split.
514: *
515: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
516: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
517: IF( INFO.NE.0 ) THEN
518: ILST = HERE
519: RETURN
520: END IF
521: HERE = HERE - 1
522: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
523: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
524: IF( INFO.NE.0 ) THEN
525: ILST = HERE
526: RETURN
527: END IF
528: HERE = HERE - 1
529: END IF
530: END IF
531: END IF
532: IF( HERE.GT.ILST )
533: $ GO TO 20
534: END IF
535: ILST = HERE
536: WORK( 1 ) = LWMIN
537: RETURN
538: *
539: * End of DTGEXC
540: *
541: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack