Annotation of rpl/lapack/lapack/dtgexc.f, revision 1.10
1.10 ! bertrand 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 November 2011
! 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: * =====================================================================
1.1 bertrand 220: SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
221: $ LDZ, IFST, ILST, WORK, LWORK, INFO )
222: *
1.10 ! bertrand 223: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 224: * -- LAPACK is a software package provided by Univ. of Tennessee, --
225: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 226: * November 2011
1.1 bertrand 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