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