Annotation of rpl/lapack/lapack/ztgex2.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b ZTGEX2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTGEX2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgex2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgex2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgex2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
! 22: * LDZ, J1, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * LOGICAL WANTQ, WANTZ
! 26: * INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
! 30: * $ Z( LDZ, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
! 40: *> in an upper triangular matrix pair (A, B) by an unitary equivalence
! 41: *> transformation.
! 42: *>
! 43: *> (A, B) must be in generalized Schur canonical form, that is, A and
! 44: *> B are both upper triangular.
! 45: *>
! 46: *> Optionally, the matrices Q and Z of generalized Schur vectors are
! 47: *> updated.
! 48: *>
! 49: *> Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
! 50: *> Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
! 51: *>
! 52: *> \endverbatim
! 53: *
! 54: * Arguments:
! 55: * ==========
! 56: *
! 57: *> \param[in] WANTQ
! 58: *> \verbatim
! 59: *> WANTQ is LOGICAL
! 60: *> .TRUE. : update the left transformation matrix Q;
! 61: *> .FALSE.: do not update Q.
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in] WANTZ
! 65: *> \verbatim
! 66: *> WANTZ is LOGICAL
! 67: *> .TRUE. : update the right transformation matrix Z;
! 68: *> .FALSE.: do not update Z.
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] N
! 72: *> \verbatim
! 73: *> N is INTEGER
! 74: *> The order of the matrices A and B. N >= 0.
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in,out] A
! 78: *> \verbatim
! 79: *> A is COMPLEX*16 arrays, dimensions (LDA,N)
! 80: *> On entry, the matrix A in the pair (A, B).
! 81: *> On exit, the updated matrix A.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] LDA
! 85: *> \verbatim
! 86: *> LDA is INTEGER
! 87: *> The leading dimension of the array A. LDA >= max(1,N).
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[in,out] B
! 91: *> \verbatim
! 92: *> B is COMPLEX*16 arrays, dimensions (LDB,N)
! 93: *> On entry, the matrix B in the pair (A, B).
! 94: *> On exit, the updated matrix B.
! 95: *> \endverbatim
! 96: *>
! 97: *> \param[in] LDB
! 98: *> \verbatim
! 99: *> LDB is INTEGER
! 100: *> The leading dimension of the array B. LDB >= max(1,N).
! 101: *> \endverbatim
! 102: *>
! 103: *> \param[in,out] Q
! 104: *> \verbatim
! 105: *> Q is COMPLEX*16 array, dimension (LDZ,N)
! 106: *> If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
! 107: *> the updated matrix Q.
! 108: *> Not referenced if WANTQ = .FALSE..
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in] LDQ
! 112: *> \verbatim
! 113: *> LDQ is INTEGER
! 114: *> The leading dimension of the array Q. LDQ >= 1;
! 115: *> If WANTQ = .TRUE., LDQ >= N.
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[in,out] Z
! 119: *> \verbatim
! 120: *> Z is COMPLEX*16 array, dimension (LDZ,N)
! 121: *> If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
! 122: *> the updated matrix Z.
! 123: *> Not referenced if WANTZ = .FALSE..
! 124: *> \endverbatim
! 125: *>
! 126: *> \param[in] LDZ
! 127: *> \verbatim
! 128: *> LDZ is INTEGER
! 129: *> The leading dimension of the array Z. LDZ >= 1;
! 130: *> If WANTZ = .TRUE., LDZ >= N.
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in] J1
! 134: *> \verbatim
! 135: *> J1 is INTEGER
! 136: *> The index to the first block (A11, B11).
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[out] INFO
! 140: *> \verbatim
! 141: *> INFO is INTEGER
! 142: *> =0: Successful exit.
! 143: *> =1: The transformed matrix pair (A, B) would be too far
! 144: *> from generalized Schur form; the problem is ill-
! 145: *> conditioned.
! 146: *> \endverbatim
! 147: *
! 148: * Authors:
! 149: * ========
! 150: *
! 151: *> \author Univ. of Tennessee
! 152: *> \author Univ. of California Berkeley
! 153: *> \author Univ. of Colorado Denver
! 154: *> \author NAG Ltd.
! 155: *
! 156: *> \date November 2011
! 157: *
! 158: *> \ingroup complex16GEauxiliary
! 159: *
! 160: *> \par Further Details:
! 161: * =====================
! 162: *>
! 163: *> In the current code both weak and strong stability tests are
! 164: *> performed. The user can omit the strong stability test by changing
! 165: *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
! 166: *> details.
! 167: *
! 168: *> \par Contributors:
! 169: * ==================
! 170: *>
! 171: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
! 172: *> Umea University, S-901 87 Umea, Sweden.
! 173: *
! 174: *> \par References:
! 175: * ================
! 176: *>
! 177: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
! 178: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
! 179: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
! 180: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
! 181: *> \n
! 182: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
! 183: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
! 184: *> Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
! 185: *> Department of Computing Science, Umea University, S-901 87 Umea,
! 186: *> Sweden, 1994. Also as LAPACK Working Note 87. To appear in
! 187: *> Numerical Algorithms, 1996.
! 188: *>
! 189: * =====================================================================
1.1 bertrand 190: SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
191: $ LDZ, J1, INFO )
192: *
1.10 ! bertrand 193: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 194: * -- LAPACK is a software package provided by Univ. of Tennessee, --
195: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 196: * November 2011
1.1 bertrand 197: *
198: * .. Scalar Arguments ..
199: LOGICAL WANTQ, WANTZ
200: INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
201: * ..
202: * .. Array Arguments ..
203: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
204: $ Z( LDZ, * )
205: * ..
206: *
207: * =====================================================================
208: *
209: * .. Parameters ..
210: COMPLEX*16 CZERO, CONE
211: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
212: $ CONE = ( 1.0D+0, 0.0D+0 ) )
1.5 bertrand 213: DOUBLE PRECISION TWENTY
214: PARAMETER ( TWENTY = 2.0D+1 )
1.1 bertrand 215: INTEGER LDST
216: PARAMETER ( LDST = 2 )
217: LOGICAL WANDS
218: PARAMETER ( WANDS = .TRUE. )
219: * ..
220: * .. Local Scalars ..
221: LOGICAL DTRONG, WEAK
222: INTEGER I, M
223: DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
224: $ THRESH, WS
225: COMPLEX*16 CDUM, F, G, SQ, SZ
226: * ..
227: * .. Local Arrays ..
228: COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
229: * ..
230: * .. External Functions ..
231: DOUBLE PRECISION DLAMCH
232: EXTERNAL DLAMCH
233: * ..
234: * .. External Subroutines ..
235: EXTERNAL ZLACPY, ZLARTG, ZLASSQ, ZROT
236: * ..
237: * .. Intrinsic Functions ..
238: INTRINSIC ABS, DBLE, DCONJG, MAX, SQRT
239: * ..
240: * .. Executable Statements ..
241: *
242: INFO = 0
243: *
244: * Quick return if possible
245: *
246: IF( N.LE.1 )
247: $ RETURN
248: *
249: M = LDST
250: WEAK = .FALSE.
251: DTRONG = .FALSE.
252: *
253: * Make a local copy of selected block in (A, B)
254: *
255: CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
256: CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
257: *
258: * Compute the threshold for testing the acceptance of swapping.
259: *
260: EPS = DLAMCH( 'P' )
261: SMLNUM = DLAMCH( 'S' ) / EPS
262: SCALE = DBLE( CZERO )
263: SUM = DBLE( CONE )
264: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
265: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
266: CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
267: SA = SCALE*SQRT( SUM )
1.5 bertrand 268: *
269: * THRES has been changed from
270: * THRESH = MAX( TEN*EPS*SA, SMLNUM )
271: * to
272: * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
273: * on 04/01/10.
274: * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
275: * Jim Demmel and Guillaume Revy. See forum post 1783.
276: *
277: THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
1.1 bertrand 278: *
279: * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
280: * using Givens rotations and perform the swap tentatively.
281: *
282: F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
283: G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
284: SA = ABS( S( 2, 2 ) )
285: SB = ABS( T( 2, 2 ) )
286: CALL ZLARTG( G, F, CZ, SZ, CDUM )
287: SZ = -SZ
288: CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
289: CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
290: IF( SA.GE.SB ) THEN
291: CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
292: ELSE
293: CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
294: END IF
295: CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
296: CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
297: *
298: * Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
299: *
300: WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
301: WEAK = WS.LE.THRESH
302: IF( .NOT.WEAK )
303: $ GO TO 20
304: *
305: IF( WANDS ) THEN
306: *
307: * Strong stability test:
1.9 bertrand 308: * F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
1.1 bertrand 309: *
310: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
311: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
312: CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
313: CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
314: CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
315: CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
316: DO 10 I = 1, 2
317: WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
318: WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
319: WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
320: WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
321: 10 CONTINUE
322: SCALE = DBLE( CZERO )
323: SUM = DBLE( CONE )
324: CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
325: SS = SCALE*SQRT( SUM )
326: DTRONG = SS.LE.THRESH
327: IF( .NOT.DTRONG )
328: $ GO TO 20
329: END IF
330: *
331: * If the swap is accepted ("weakly" and "strongly"), apply the
332: * equivalence transformations to the original matrix pair (A,B)
333: *
334: CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
335: $ DCONJG( SZ ) )
336: CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
337: $ DCONJG( SZ ) )
338: CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
339: CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
340: *
341: * Set N1 by N2 (2,1) blocks to 0
342: *
343: A( J1+1, J1 ) = CZERO
344: B( J1+1, J1 ) = CZERO
345: *
346: * Accumulate transformations into Q and Z if requested.
347: *
348: IF( WANTZ )
349: $ CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
350: $ DCONJG( SZ ) )
351: IF( WANTQ )
352: $ CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
353: $ DCONJG( SQ ) )
354: *
355: * Exit with INFO = 0 if swap was successfully performed.
356: *
357: RETURN
358: *
359: * Exit with INFO = 1 if swap was rejected.
360: *
361: 20 CONTINUE
362: INFO = 1
363: RETURN
364: *
365: * End of ZTGEX2
366: *
367: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgex2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack