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: *> \date December 2016
199: *
200: *> \ingroup doubleGEcomputational
201: *
202: *> \par Contributors:
203: * ==================
204: *>
205: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
206: *> Umea University, S-901 87 Umea, Sweden.
207: *
208: *> \par References:
209: * ================
210: *>
211: *> \verbatim
212: *>
213: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
214: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
215: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
216: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
217: *> \endverbatim
218: *>
219: * =====================================================================
220: SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
221: $ LDZ, IFST, ILST, WORK, LWORK, INFO )
222: *
223: * -- LAPACK computational routine (version 3.7.0) --
224: * -- LAPACK is a software package provided by Univ. of Tennessee, --
225: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226: * December 2016
227: *
228: * .. Scalar Arguments ..
229: LOGICAL WANTQ, WANTZ
230: INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
231: * ..
232: * .. Array Arguments ..
233: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
234: $ WORK( * ), Z( LDZ, * )
235: * ..
236: *
237: * =====================================================================
238: *
239: * .. Parameters ..
240: DOUBLE PRECISION ZERO
241: PARAMETER ( ZERO = 0.0D+0 )
242: * ..
243: * .. Local Scalars ..
244: LOGICAL LQUERY
245: INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
246: * ..
247: * .. External Subroutines ..
248: EXTERNAL DTGEX2, XERBLA
249: * ..
250: * .. Intrinsic Functions ..
251: INTRINSIC MAX
252: * ..
253: * .. Executable Statements ..
254: *
255: * Decode and test input arguments.
256: *
257: INFO = 0
258: LQUERY = ( LWORK.EQ.-1 )
259: IF( N.LT.0 ) THEN
260: INFO = -3
261: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
262: INFO = -5
263: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
264: INFO = -7
265: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
266: INFO = -9
267: ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
268: INFO = -11
269: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
270: INFO = -12
271: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
272: INFO = -13
273: END IF
274: *
275: IF( INFO.EQ.0 ) THEN
276: IF( N.LE.1 ) THEN
277: LWMIN = 1
278: ELSE
279: LWMIN = 4*N + 16
280: END IF
281: WORK(1) = LWMIN
282: *
283: IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
284: INFO = -15
285: END IF
286: END IF
287: *
288: IF( INFO.NE.0 ) THEN
289: CALL XERBLA( 'DTGEXC', -INFO )
290: RETURN
291: ELSE IF( LQUERY ) THEN
292: RETURN
293: END IF
294: *
295: * Quick return if possible
296: *
297: IF( N.LE.1 )
298: $ RETURN
299: *
300: * Determine the first row of the specified block and find out
301: * if it is 1-by-1 or 2-by-2.
302: *
303: IF( IFST.GT.1 ) THEN
304: IF( A( IFST, IFST-1 ).NE.ZERO )
305: $ IFST = IFST - 1
306: END IF
307: NBF = 1
308: IF( IFST.LT.N ) THEN
309: IF( A( IFST+1, IFST ).NE.ZERO )
310: $ NBF = 2
311: END IF
312: *
313: * Determine the first row of the final block
314: * and find out if it is 1-by-1 or 2-by-2.
315: *
316: IF( ILST.GT.1 ) THEN
317: IF( A( ILST, ILST-1 ).NE.ZERO )
318: $ ILST = ILST - 1
319: END IF
320: NBL = 1
321: IF( ILST.LT.N ) THEN
322: IF( A( ILST+1, ILST ).NE.ZERO )
323: $ NBL = 2
324: END IF
325: IF( IFST.EQ.ILST )
326: $ RETURN
327: *
328: IF( IFST.LT.ILST ) THEN
329: *
330: * Update ILST.
331: *
332: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
333: $ ILST = ILST - 1
334: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
335: $ ILST = ILST + 1
336: *
337: HERE = IFST
338: *
339: 10 CONTINUE
340: *
341: * Swap with next one below.
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+NBF+1.LE.N ) THEN
349: IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
350: $ NBNEXT = 2
351: END IF
352: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
353: $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, 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( A( 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+3.LE.N ) THEN
374: IF( A( HERE+3, HERE+2 ).NE.ZERO )
375: $ NBNEXT = 2
376: END IF
377: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
378: $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, 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, 1, 1, WORK, LWORK, INFO )
389: IF( INFO.NE.0 ) THEN
390: ILST = HERE
391: RETURN
392: END IF
393: HERE = HERE + 1
394: *
395: ELSE
396: *
397: * Recompute NBNEXT in case of 2-by-2 split.
398: *
399: IF( A( HERE+2, HERE+1 ).EQ.ZERO )
400: $ NBNEXT = 1
401: IF( NBNEXT.EQ.2 ) THEN
402: *
403: * 2-by-2 block did not split.
404: *
405: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
406: $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
407: $ INFO )
408: IF( INFO.NE.0 ) THEN
409: ILST = HERE
410: RETURN
411: END IF
412: HERE = HERE + 2
413: ELSE
414: *
415: * 2-by-2 block did split.
416: *
417: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
418: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
419: IF( INFO.NE.0 ) THEN
420: ILST = HERE
421: RETURN
422: END IF
423: HERE = HERE + 1
424: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
425: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
426: IF( INFO.NE.0 ) THEN
427: ILST = HERE
428: RETURN
429: END IF
430: HERE = HERE + 1
431: END IF
432: *
433: END IF
434: END IF
435: IF( HERE.LT.ILST )
436: $ GO TO 10
437: ELSE
438: HERE = IFST
439: *
440: 20 CONTINUE
441: *
442: * Swap with next one below.
443: *
444: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
445: *
446: * Current block either 1-by-1 or 2-by-2.
447: *
448: NBNEXT = 1
449: IF( HERE.GE.3 ) THEN
450: IF( A( HERE-1, HERE-2 ).NE.ZERO )
451: $ NBNEXT = 2
452: END IF
453: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
454: $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
455: $ INFO )
456: IF( INFO.NE.0 ) THEN
457: ILST = HERE
458: RETURN
459: END IF
460: HERE = HERE - NBNEXT
461: *
462: * Test if 2-by-2 block breaks into two 1-by-1 blocks.
463: *
464: IF( NBF.EQ.2 ) THEN
465: IF( A( HERE+1, HERE ).EQ.ZERO )
466: $ NBF = 3
467: END IF
468: *
469: ELSE
470: *
471: * Current block consists of two 1-by-1 blocks, each of which
472: * must be swapped individually.
473: *
474: NBNEXT = 1
475: IF( HERE.GE.3 ) THEN
476: IF( A( HERE-1, HERE-2 ).NE.ZERO )
477: $ NBNEXT = 2
478: END IF
479: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
480: $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
481: $ INFO )
482: IF( INFO.NE.0 ) THEN
483: ILST = HERE
484: RETURN
485: END IF
486: IF( NBNEXT.EQ.1 ) THEN
487: *
488: * Swap two 1-by-1 blocks.
489: *
490: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
491: $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
492: IF( INFO.NE.0 ) THEN
493: ILST = HERE
494: RETURN
495: END IF
496: HERE = HERE - 1
497: ELSE
498: *
499: * Recompute NBNEXT in case of 2-by-2 split.
500: *
501: IF( A( HERE, HERE-1 ).EQ.ZERO )
502: $ NBNEXT = 1
503: IF( NBNEXT.EQ.2 ) THEN
504: *
505: * 2-by-2 block did not split.
506: *
507: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
508: $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
509: IF( INFO.NE.0 ) THEN
510: ILST = HERE
511: RETURN
512: END IF
513: HERE = HERE - 2
514: ELSE
515: *
516: * 2-by-2 block did split.
517: *
518: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
519: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
520: IF( INFO.NE.0 ) THEN
521: ILST = HERE
522: RETURN
523: END IF
524: HERE = HERE - 1
525: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
526: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
527: IF( INFO.NE.0 ) THEN
528: ILST = HERE
529: RETURN
530: END IF
531: HERE = HERE - 1
532: END IF
533: END IF
534: END IF
535: IF( HERE.GT.ILST )
536: $ GO TO 20
537: END IF
538: ILST = HERE
539: WORK( 1 ) = LWMIN
540: RETURN
541: *
542: * End of DTGEXC
543: *
544: 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