1: SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
2: $ LDZ, IFST, ILST, WORK, LWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: LOGICAL WANTQ, WANTZ
11: INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
15: $ WORK( * ), Z( LDZ, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DTGEXC reorders the generalized real Schur decomposition of a real
22: * matrix pair (A,B) using an orthogonal equivalence transformation
23: *
24: * (A, B) = Q * (A, B) * Z',
25: *
26: * so that the diagonal block of (A, B) with row index IFST is moved
27: * to row ILST.
28: *
29: * (A, B) must be in generalized real Schur canonical form (as returned
30: * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
31: * diagonal blocks. B is upper triangular.
32: *
33: * Optionally, the matrices Q and Z of generalized Schur vectors are
34: * updated.
35: *
36: * Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
37: * Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
38: *
39: *
40: * Arguments
41: * =========
42: *
43: * WANTQ (input) LOGICAL
44: * .TRUE. : update the left transformation matrix Q;
45: * .FALSE.: do not update Q.
46: *
47: * WANTZ (input) LOGICAL
48: * .TRUE. : update the right transformation matrix Z;
49: * .FALSE.: do not update Z.
50: *
51: * N (input) INTEGER
52: * The order of the matrices A and B. N >= 0.
53: *
54: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
55: * On entry, the matrix A in generalized real Schur canonical
56: * form.
57: * On exit, the updated matrix A, again in generalized
58: * real Schur canonical form.
59: *
60: * LDA (input) INTEGER
61: * The leading dimension of the array A. LDA >= max(1,N).
62: *
63: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
64: * On entry, the matrix B in generalized real Schur canonical
65: * form (A,B).
66: * On exit, the updated matrix B, again in generalized
67: * real Schur canonical form (A,B).
68: *
69: * LDB (input) INTEGER
70: * The leading dimension of the array B. LDB >= max(1,N).
71: *
72: * Q (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
73: * On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
74: * On exit, the updated matrix Q.
75: * If WANTQ = .FALSE., Q is not referenced.
76: *
77: * LDQ (input) INTEGER
78: * The leading dimension of the array Q. LDQ >= 1.
79: * If WANTQ = .TRUE., LDQ >= N.
80: *
81: * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
82: * On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
83: * On exit, the updated matrix Z.
84: * If WANTZ = .FALSE., Z is not referenced.
85: *
86: * LDZ (input) INTEGER
87: * The leading dimension of the array Z. LDZ >= 1.
88: * If WANTZ = .TRUE., LDZ >= N.
89: *
90: * IFST (input/output) INTEGER
91: * ILST (input/output) INTEGER
92: * Specify the reordering of the diagonal blocks of (A, B).
93: * The block with row index IFST is moved to row ILST, by a
94: * sequence of swapping between adjacent blocks.
95: * On exit, if IFST pointed on entry to the second row of
96: * a 2-by-2 block, it is changed to point to the first row;
97: * ILST always points to the first row of the block in its
98: * final position (which may differ from its input value by
99: * +1 or -1). 1 <= IFST, ILST <= N.
100: *
101: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
102: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103: *
104: * LWORK (input) INTEGER
105: * The dimension of the array WORK.
106: * LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
107: *
108: * If LWORK = -1, then a workspace query is assumed; the routine
109: * only calculates the optimal size of the WORK array, returns
110: * this value as the first entry of the WORK array, and no error
111: * message related to LWORK is issued by XERBLA.
112: *
113: * INFO (output) INTEGER
114: * =0: successful exit.
115: * <0: if INFO = -i, the i-th argument had an illegal value.
116: * =1: The transformed matrix pair (A, B) would be too far
117: * from generalized Schur form; the problem is ill-
118: * conditioned. (A, B) may have been partially reordered,
119: * and ILST points to the first row of the current
120: * position of the block being moved.
121: *
122: * Further Details
123: * ===============
124: *
125: * Based on contributions by
126: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
127: * Umea University, S-901 87 Umea, Sweden.
128: *
129: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
130: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
131: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
132: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
133: *
134: * =====================================================================
135: *
136: * .. Parameters ..
137: DOUBLE PRECISION ZERO
138: PARAMETER ( ZERO = 0.0D+0 )
139: * ..
140: * .. Local Scalars ..
141: LOGICAL LQUERY
142: INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
143: * ..
144: * .. External Subroutines ..
145: EXTERNAL DTGEX2, XERBLA
146: * ..
147: * .. Intrinsic Functions ..
148: INTRINSIC MAX
149: * ..
150: * .. Executable Statements ..
151: *
152: * Decode and test input arguments.
153: *
154: INFO = 0
155: LQUERY = ( LWORK.EQ.-1 )
156: IF( N.LT.0 ) THEN
157: INFO = -3
158: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
159: INFO = -5
160: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
161: INFO = -7
162: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
163: INFO = -9
164: ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
165: INFO = -11
166: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
167: INFO = -12
168: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
169: INFO = -13
170: END IF
171: *
172: IF( INFO.EQ.0 ) THEN
173: IF( N.LE.1 ) THEN
174: LWMIN = 1
175: ELSE
176: LWMIN = 4*N + 16
177: END IF
178: WORK(1) = LWMIN
179: *
180: IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
181: INFO = -15
182: END IF
183: END IF
184: *
185: IF( INFO.NE.0 ) THEN
186: CALL XERBLA( 'DTGEXC', -INFO )
187: RETURN
188: ELSE IF( LQUERY ) THEN
189: RETURN
190: END IF
191: *
192: * Quick return if possible
193: *
194: IF( N.LE.1 )
195: $ RETURN
196: *
197: * Determine the first row of the specified block and find out
198: * if it is 1-by-1 or 2-by-2.
199: *
200: IF( IFST.GT.1 ) THEN
201: IF( A( IFST, IFST-1 ).NE.ZERO )
202: $ IFST = IFST - 1
203: END IF
204: NBF = 1
205: IF( IFST.LT.N ) THEN
206: IF( A( IFST+1, IFST ).NE.ZERO )
207: $ NBF = 2
208: END IF
209: *
210: * Determine the first row of the final block
211: * and find out if it is 1-by-1 or 2-by-2.
212: *
213: IF( ILST.GT.1 ) THEN
214: IF( A( ILST, ILST-1 ).NE.ZERO )
215: $ ILST = ILST - 1
216: END IF
217: NBL = 1
218: IF( ILST.LT.N ) THEN
219: IF( A( ILST+1, ILST ).NE.ZERO )
220: $ NBL = 2
221: END IF
222: IF( IFST.EQ.ILST )
223: $ RETURN
224: *
225: IF( IFST.LT.ILST ) THEN
226: *
227: * Update ILST.
228: *
229: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
230: $ ILST = ILST - 1
231: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
232: $ ILST = ILST + 1
233: *
234: HERE = IFST
235: *
236: 10 CONTINUE
237: *
238: * Swap with next one below.
239: *
240: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
241: *
242: * Current block either 1-by-1 or 2-by-2.
243: *
244: NBNEXT = 1
245: IF( HERE+NBF+1.LE.N ) THEN
246: IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
247: $ NBNEXT = 2
248: END IF
249: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
250: $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
251: IF( INFO.NE.0 ) THEN
252: ILST = HERE
253: RETURN
254: END IF
255: HERE = HERE + NBNEXT
256: *
257: * Test if 2-by-2 block breaks into two 1-by-1 blocks.
258: *
259: IF( NBF.EQ.2 ) THEN
260: IF( A( HERE+1, HERE ).EQ.ZERO )
261: $ NBF = 3
262: END IF
263: *
264: ELSE
265: *
266: * Current block consists of two 1-by-1 blocks, each of which
267: * must be swapped individually.
268: *
269: NBNEXT = 1
270: IF( HERE+3.LE.N ) THEN
271: IF( A( HERE+3, HERE+2 ).NE.ZERO )
272: $ NBNEXT = 2
273: END IF
274: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
275: $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
276: IF( INFO.NE.0 ) THEN
277: ILST = HERE
278: RETURN
279: END IF
280: IF( NBNEXT.EQ.1 ) THEN
281: *
282: * Swap two 1-by-1 blocks.
283: *
284: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
285: $ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
286: IF( INFO.NE.0 ) THEN
287: ILST = HERE
288: RETURN
289: END IF
290: HERE = HERE + 1
291: *
292: ELSE
293: *
294: * Recompute NBNEXT in case of 2-by-2 split.
295: *
296: IF( A( HERE+2, HERE+1 ).EQ.ZERO )
297: $ NBNEXT = 1
298: IF( NBNEXT.EQ.2 ) THEN
299: *
300: * 2-by-2 block did not split.
301: *
302: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
303: $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
304: $ INFO )
305: IF( INFO.NE.0 ) THEN
306: ILST = HERE
307: RETURN
308: END IF
309: HERE = HERE + 2
310: ELSE
311: *
312: * 2-by-2 block did split.
313: *
314: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
315: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
316: IF( INFO.NE.0 ) THEN
317: ILST = HERE
318: RETURN
319: END IF
320: HERE = HERE + 1
321: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
322: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
323: IF( INFO.NE.0 ) THEN
324: ILST = HERE
325: RETURN
326: END IF
327: HERE = HERE + 1
328: END IF
329: *
330: END IF
331: END IF
332: IF( HERE.LT.ILST )
333: $ GO TO 10
334: ELSE
335: HERE = IFST
336: *
337: 20 CONTINUE
338: *
339: * Swap with next one below.
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( A( HERE-1, HERE-2 ).NE.ZERO )
348: $ NBNEXT = 2
349: END IF
350: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
351: $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
352: $ INFO )
353: IF( INFO.NE.0 ) THEN
354: ILST = HERE
355: RETURN
356: END IF
357: HERE = HERE - NBNEXT
358: *
359: * Test if 2-by-2 block breaks into two 1-by-1 blocks.
360: *
361: IF( NBF.EQ.2 ) THEN
362: IF( A( HERE+1, HERE ).EQ.ZERO )
363: $ NBF = 3
364: END IF
365: *
366: ELSE
367: *
368: * Current block consists of two 1-by-1 blocks, each of which
369: * must be swapped individually.
370: *
371: NBNEXT = 1
372: IF( HERE.GE.3 ) THEN
373: IF( A( HERE-1, HERE-2 ).NE.ZERO )
374: $ NBNEXT = 2
375: END IF
376: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
377: $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
378: $ 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.
386: *
387: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
388: $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
389: IF( INFO.NE.0 ) THEN
390: ILST = HERE
391: RETURN
392: END IF
393: HERE = HERE - 1
394: ELSE
395: *
396: * Recompute NBNEXT in case of 2-by-2 split.
397: *
398: IF( A( HERE, HERE-1 ).EQ.ZERO )
399: $ NBNEXT = 1
400: IF( NBNEXT.EQ.2 ) THEN
401: *
402: * 2-by-2 block did not split.
403: *
404: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
405: $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
406: IF( INFO.NE.0 ) THEN
407: ILST = HERE
408: RETURN
409: END IF
410: HERE = HERE - 2
411: ELSE
412: *
413: * 2-by-2 block did split.
414: *
415: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
416: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
417: IF( INFO.NE.0 ) THEN
418: ILST = HERE
419: RETURN
420: END IF
421: HERE = HERE - 1
422: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
423: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
424: IF( INFO.NE.0 ) THEN
425: ILST = HERE
426: RETURN
427: END IF
428: HERE = HERE - 1
429: END IF
430: END IF
431: END IF
432: IF( HERE.GT.ILST )
433: $ GO TO 20
434: END IF
435: ILST = HERE
436: WORK( 1 ) = LWMIN
437: RETURN
438: *
439: * End of DTGEXC
440: *
441: 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