1       SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
  2      $                   LDZ, J1, INFO )
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.3.1) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *  -- April 2011                                                      --
  8 *
  9 *     .. Scalar Arguments ..
 10       LOGICAL            WANTQ, WANTZ
 11       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
 12 *     ..
 13 *     .. Array Arguments ..
 14       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 15      $                   Z( LDZ, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
 22 *  in an upper triangular matrix pair (A, B) by an unitary equivalence
 23 *  transformation.
 24 *
 25 *  (A, B) must be in generalized Schur canonical form, that is, A and
 26 *  B are both upper triangular.
 27 *
 28 *  Optionally, the matrices Q and Z of generalized Schur vectors are
 29 *  updated.
 30 *
 31 *         Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
 32 *         Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
 33 *
 34 *
 35 *  Arguments
 36 *  =========
 37 *
 38 *  WANTQ   (input) LOGICAL
 39 *          .TRUE. : update the left transformation matrix Q;
 40 *          .FALSE.: do not update Q.
 41 *
 42 *  WANTZ   (input) LOGICAL
 43 *          .TRUE. : update the right transformation matrix Z;
 44 *          .FALSE.: do not update Z.
 45 *
 46 *  N       (input) INTEGER
 47 *          The order of the matrices A and B. N >= 0.
 48 *
 49 *  A       (input/output) COMPLEX*16 arrays, dimensions (LDA,N)
 50 *          On entry, the matrix A in the pair (A, B).
 51 *          On exit, the updated matrix A.
 52 *
 53 *  LDA     (input)  INTEGER
 54 *          The leading dimension of the array A. LDA >= max(1,N).
 55 *
 56 *  B       (input/output) COMPLEX*16 arrays, dimensions (LDB,N)
 57 *          On entry, the matrix B in the pair (A, B).
 58 *          On exit, the updated matrix B.
 59 *
 60 *  LDB     (input)  INTEGER
 61 *          The leading dimension of the array B. LDB >= max(1,N).
 62 *
 63 *  Q       (input/output) COMPLEX*16 array, dimension (LDZ,N)
 64 *          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
 65 *          the updated matrix Q.
 66 *          Not referenced if WANTQ = .FALSE..
 67 *
 68 *  LDQ     (input) INTEGER
 69 *          The leading dimension of the array Q. LDQ >= 1;
 70 *          If WANTQ = .TRUE., LDQ >= N.
 71 *
 72 *  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
 73 *          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
 74 *          the updated matrix Z.
 75 *          Not referenced if WANTZ = .FALSE..
 76 *
 77 *  LDZ     (input) INTEGER
 78 *          The leading dimension of the array Z. LDZ >= 1;
 79 *          If WANTZ = .TRUE., LDZ >= N.
 80 *
 81 *  J1      (input) INTEGER
 82 *          The index to the first block (A11, B11).
 83 *
 84 *  INFO    (output) INTEGER
 85 *           =0:  Successful exit.
 86 *           =1:  The transformed matrix pair (A, B) would be too far
 87 *                from generalized Schur form; the problem is ill-
 88 *                conditioned. 
 89 *
 90 *
 91 *  Further Details
 92 *  ===============
 93 *
 94 *  Based on contributions by
 95 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 96 *     Umea University, S-901 87 Umea, Sweden.
 97 *
 98 *  In the current code both weak and strong stability tests are
 99 *  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+00.0D+0 ),
120      $                   CONE = ( 1.0D+00.0D+0 ) )
121       DOUBLE PRECISION   TWENTY
122       PARAMETER          ( TWENTY = 2.0D+1 )
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          ABSDBLEDCONJGMAXSQRT
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, 1SCALESUM )
175       SA = SCALE*SQRTSUM )
176 *
177 *     THRES has been changed from 
178 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
179 *     to
180 *        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
181 *     on 04/01/10.
182 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
183 *     Jim Demmel and Guillaume Revy. See forum post 1783.
184 *
185       THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
186 *
187 *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
188 *     using Givens rotations and perform the swap tentatively.
189 *
190       F = S( 22 )*T( 11 ) - T( 22 )*S( 11 )
191       G = S( 22 )*T( 12 ) - T( 22 )*S( 12 )
192       SA = ABS( S( 22 ) )
193       SB = ABS( T( 22 ) )
194       CALL ZLARTG( G, F, CZ, SZ, CDUM )
195       SZ = -SZ
196       CALL ZROT( 2, S( 11 ), 1, S( 12 ), 1, CZ, DCONJG( SZ ) )
197       CALL ZROT( 2, T( 11 ), 1, T( 12 ), 1, CZ, DCONJG( SZ ) )
198       IF( SA.GE.SB ) THEN
199          CALL ZLARTG( S( 11 ), S( 21 ), CQ, SQ, CDUM )
200       ELSE
201          CALL ZLARTG( T( 11 ), T( 21 ), CQ, SQ, CDUM )
202       END IF
203       CALL ZROT( 2, S( 11 ), LDST, S( 21 ), LDST, CQ, SQ )
204       CALL ZROT( 2, T( 11 ), LDST, T( 21 ), LDST, CQ, SQ )
205 *
206 *     Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
207 *
208       WS = ABS( S( 21 ) ) + ABS( T( 21 ) )
209       WEAK = WS.LE.THRESH
210       IF.NOT.WEAK )
211      $   GO TO 20
212 *
213       IF( WANDS ) THEN
214 *
215 *        Strong stability test:
216 *           F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
217 *
218          CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
219          CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
220          CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
221          CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
222          CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
223          CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
224          DO 10 I = 12
225             WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
226             WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
227             WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
228             WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
229    10    CONTINUE
230          SCALE = DBLE( CZERO )
231          SUM = DBLE( CONE )
232          CALL ZLASSQ( 2*M*M, WORK, 1SCALESUM )
233          SS = SCALE*SQRTSUM )
234          DTRONG = SS.LE.THRESH
235          IF.NOT.DTRONG )
236      $      GO TO 20
237       END IF
238 *
239 *     If the swap is accepted ("weakly" and "strongly"), apply the
240 *     equivalence transformations to the original matrix pair (A,B)
241 *
242       CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
243      $           DCONJG( SZ ) )
244       CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
245      $           DCONJG( SZ ) )
246       CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
247       CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
248 *
249 *     Set  N1 by N2 (2,1) blocks to 0
250 *
251       A( J1+1, J1 ) = CZERO
252       B( J1+1, J1 ) = CZERO
253 *
254 *     Accumulate transformations into Q and Z if requested.
255 *
256       IF( WANTZ )
257      $   CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
258      $              DCONJG( SZ ) )
259       IF( WANTQ )
260      $   CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
261      $              DCONJG( SQ ) )
262 *
263 *     Exit with INFO = 0 if swap was successfully performed.
264 *
265       RETURN
266 *
267 *     Exit with INFO = 1 if swap was rejected.
268 *
269    20 CONTINUE
270       INFO = 1
271       RETURN
272 *
273 *     End of ZTGEX2
274 *
275       END