001:       SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
002:      $                   LDZ, J1, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            WANTQ, WANTZ
011:       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
012: *     ..
013: *     .. Array Arguments ..
014:       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
015:      $                   Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
022: *  in an upper triangular matrix pair (A, B) by an unitary equivalence
023: *  transformation.
024: *
025: *  (A, B) must be in generalized Schur canonical form, that is, A and
026: *  B are both upper triangular.
027: *
028: *  Optionally, the matrices Q and Z of generalized Schur vectors are
029: *  updated.
030: *
031: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
032: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
033: *
034: *
035: *  Arguments
036: *  =========
037: *
038: *  WANTQ   (input) LOGICAL
039: *          .TRUE. : update the left transformation matrix Q;
040: *          .FALSE.: do not update Q.
041: *
042: *  WANTZ   (input) LOGICAL
043: *          .TRUE. : update the right transformation matrix Z;
044: *          .FALSE.: do not update Z.
045: *
046: *  N       (input) INTEGER
047: *          The order of the matrices A and B. N >= 0.
048: *
049: *  A       (input/output) COMPLEX*16 arrays, dimensions (LDA,N)
050: *          On entry, the matrix A in the pair (A, B).
051: *          On exit, the updated matrix A.
052: *
053: *  LDA     (input)  INTEGER
054: *          The leading dimension of the array A. LDA >= max(1,N).
055: *
056: *  B       (input/output) COMPLEX*16 arrays, dimensions (LDB,N)
057: *          On entry, the matrix B in the pair (A, B).
058: *          On exit, the updated matrix B.
059: *
060: *  LDB     (input)  INTEGER
061: *          The leading dimension of the array B. LDB >= max(1,N).
062: *
063: *  Q       (input/output) COMPLEX*16 array, dimension (LDZ,N)
064: *          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
065: *          the updated matrix Q.
066: *          Not referenced if WANTQ = .FALSE..
067: *
068: *  LDQ     (input) INTEGER
069: *          The leading dimension of the array Q. LDQ >= 1;
070: *          If WANTQ = .TRUE., LDQ >= N.
071: *
072: *  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
073: *          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
074: *          the updated matrix Z.
075: *          Not referenced if WANTZ = .FALSE..
076: *
077: *  LDZ     (input) INTEGER
078: *          The leading dimension of the array Z. LDZ >= 1;
079: *          If WANTZ = .TRUE., LDZ >= N.
080: *
081: *  J1      (input) INTEGER
082: *          The index to the first block (A11, B11).
083: *
084: *  INFO    (output) INTEGER
085: *           =0:  Successful exit.
086: *           =1:  The transformed matrix pair (A, B) would be too far
087: *                from generalized Schur form; the problem is ill-
088: *                conditioned. 
089: *
090: *
091: *  Further Details
092: *  ===============
093: *
094: *  Based on contributions by
095: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
096: *     Umea University, S-901 87 Umea, Sweden.
097: *
098: *  In the current code both weak and strong stability tests are
099: *  performed. The user can omit the strong stability test by changing
100: *  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
101: *  details.
102: *
103: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
104: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
105: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
106: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
107: *
108: *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
109: *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
110: *      Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
111: *      Department of Computing Science, Umea University, S-901 87 Umea,
112: *      Sweden, 1994. Also as LAPACK Working Note 87. To appear in
113: *      Numerical Algorithms, 1996.
114: *
115: *  =====================================================================
116: *
117: *     .. Parameters ..
118:       COMPLEX*16         CZERO, CONE
119:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
120:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
121:       DOUBLE PRECISION   TEN
122:       PARAMETER          ( TEN = 10.0D+0 )
123:       INTEGER            LDST
124:       PARAMETER          ( LDST = 2 )
125:       LOGICAL            WANDS
126:       PARAMETER          ( WANDS = .TRUE. )
127: *     ..
128: *     .. Local Scalars ..
129:       LOGICAL            DTRONG, WEAK
130:       INTEGER            I, M
131:       DOUBLE PRECISION   CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
132:      $                   THRESH, WS
133:       COMPLEX*16         CDUM, F, G, SQ, SZ
134: *     ..
135: *     .. Local Arrays ..
136:       COMPLEX*16         S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
137: *     ..
138: *     .. External Functions ..
139:       DOUBLE PRECISION   DLAMCH
140:       EXTERNAL           DLAMCH
141: *     ..
142: *     .. External Subroutines ..
143:       EXTERNAL           ZLACPY, ZLARTG, ZLASSQ, ZROT
144: *     ..
145: *     .. Intrinsic Functions ..
146:       INTRINSIC          ABS, DBLE, DCONJG, MAX, SQRT
147: *     ..
148: *     .. Executable Statements ..
149: *
150:       INFO = 0
151: *
152: *     Quick return if possible
153: *
154:       IF( N.LE.1 )
155:      $   RETURN
156: *
157:       M = LDST
158:       WEAK = .FALSE.
159:       DTRONG = .FALSE.
160: *
161: *     Make a local copy of selected block in (A, B)
162: *
163:       CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
164:       CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
165: *
166: *     Compute the threshold for testing the acceptance of swapping.
167: *
168:       EPS = DLAMCH( 'P' )
169:       SMLNUM = DLAMCH( 'S' ) / EPS
170:       SCALE = DBLE( CZERO )
171:       SUM = DBLE( CONE )
172:       CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
173:       CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
174:       CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
175:       SA = SCALE*SQRT( SUM )
176:       THRESH = MAX( TEN*EPS*SA, SMLNUM )
177: *
178: *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
179: *     using Givens rotations and perform the swap tentatively.
180: *
181:       F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
182:       G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
183:       SA = ABS( S( 2, 2 ) )
184:       SB = ABS( T( 2, 2 ) )
185:       CALL ZLARTG( G, F, CZ, SZ, CDUM )
186:       SZ = -SZ
187:       CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
188:       CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
189:       IF( SA.GE.SB ) THEN
190:          CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
191:       ELSE
192:          CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
193:       END IF
194:       CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
195:       CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
196: *
197: *     Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
198: *
199:       WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
200:       WEAK = WS.LE.THRESH
201:       IF( .NOT.WEAK )
202:      $   GO TO 20
203: *
204:       IF( WANDS ) THEN
205: *
206: *        Strong stability test:
207: *           F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A, B)))
208: *
209:          CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
210:          CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
211:          CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
212:          CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
213:          CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
214:          CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
215:          DO 10 I = 1, 2
216:             WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
217:             WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
218:             WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
219:             WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
220:    10    CONTINUE
221:          SCALE = DBLE( CZERO )
222:          SUM = DBLE( CONE )
223:          CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
224:          SS = SCALE*SQRT( SUM )
225:          DTRONG = SS.LE.THRESH
226:          IF( .NOT.DTRONG )
227:      $      GO TO 20
228:       END IF
229: *
230: *     If the swap is accepted ("weakly" and "strongly"), apply the
231: *     equivalence transformations to the original matrix pair (A,B)
232: *
233:       CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
234:      $           DCONJG( SZ ) )
235:       CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
236:      $           DCONJG( SZ ) )
237:       CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
238:       CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
239: *
240: *     Set  N1 by N2 (2,1) blocks to 0
241: *
242:       A( J1+1, J1 ) = CZERO
243:       B( J1+1, J1 ) = CZERO
244: *
245: *     Accumulate transformations into Q and Z if requested.
246: *
247:       IF( WANTZ )
248:      $   CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
249:      $              DCONJG( SZ ) )
250:       IF( WANTQ )
251:      $   CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
252:      $              DCONJG( SQ ) )
253: *
254: *     Exit with INFO = 0 if swap was successfully performed.
255: *
256:       RETURN
257: *
258: *     Exit with INFO = 1 if swap was rejected.
259: *
260:    20 CONTINUE
261:       INFO = 1
262:       RETURN
263: *
264: *     End of ZTGEX2
265: *
266:       END
267: