 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ ctgex2()

 subroutine ctgex2 ( logical WANTQ, logical WANTZ, integer N, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldq, * ) Q, integer LDQ, complex, dimension( ldz, * ) Z, integer LDZ, integer J1, integer INFO )

CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.

Purpose:
``` CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
in an upper triangular matrix pair (A, B) by an unitary equivalence
transformation.

(A, B) must be in generalized Schur canonical form, that is, A and
B are both upper triangular.

Optionally, the matrices Q and Z of generalized Schur vectors are
updated.

Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H```
Parameters
 [in] WANTQ ``` WANTQ is LOGICAL .TRUE. : update the left transformation matrix Q; .FALSE.: do not update Q.``` [in] WANTZ ``` WANTZ is LOGICAL .TRUE. : update the right transformation matrix Z; .FALSE.: do not update Z.``` [in] N ``` N is INTEGER The order of the matrices A and B. N >= 0.``` [in,out] A ``` A is COMPLEX array, dimension (LDA,N) On entry, the matrix A in the pair (A, B). On exit, the updated matrix A.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in,out] B ``` B is COMPLEX array, dimension (LDB,N) On entry, the matrix B in the pair (A, B). On exit, the updated matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [in,out] Q ``` Q is COMPLEX array, dimension (LDQ,N) If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit, the updated matrix Q. Not referenced if WANTQ = .FALSE..``` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= 1; If WANTQ = .TRUE., LDQ >= N.``` [in,out] Z ``` Z is COMPLEX array, dimension (LDZ,N) If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit, the updated matrix Z. Not referenced if WANTZ = .FALSE..``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1; If WANTZ = .TRUE., LDZ >= N.``` [in] J1 ``` J1 is INTEGER The index to the first block (A11, B11).``` [out] INFO ``` INFO is INTEGER =0: Successful exit. =1: The transformed matrix pair (A, B) would be too far from generalized Schur form; the problem is ill- conditioned.```
Further Details:
In the current code both weak and strong stability tests are performed. The user can omit the strong stability test by changing the internal logical parameter WANDS to .FALSE.. See ref.  for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
 B. Kagstrom; A Direct Method for Reordering Eigenvalues in the Generalized Real Schur Form of a Regular Matrix Pair (A, B), in M.S. Moonen et al (eds), Linear Algebra for Large Scale and Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
 B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified Eigenvalues of a Regular Matrix Pair (A, B) and Condition Estimation: Theory, Algorithms and Software, Report UMINF-94.04, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. To appear in Numerical Algorithms, 1996.

Definition at line 188 of file ctgex2.f.

190*
191* -- LAPACK auxiliary routine --
192* -- LAPACK is a software package provided by Univ. of Tennessee, --
193* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194*
195* .. Scalar Arguments ..
196 LOGICAL WANTQ, WANTZ
197 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
198* ..
199* .. Array Arguments ..
200 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201 \$ Z( LDZ, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 COMPLEX CZERO, CONE
208 parameter( czero = ( 0.0e+0, 0.0e+0 ),
209 \$ cone = ( 1.0e+0, 0.0e+0 ) )
210 REAL TWENTY
211 parameter( twenty = 2.0e+1 )
212 INTEGER LDST
213 parameter( ldst = 2 )
214 LOGICAL WANDS
215 parameter( wands = .true. )
216* ..
217* .. Local Scalars ..
218 LOGICAL STRONG, WEAK
219 INTEGER I, M
220 REAL CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221 \$ THRESHA, THRESHB
222 COMPLEX CDUM, F, G, SQ, SZ
223* ..
224* .. Local Arrays ..
225 COMPLEX S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226* ..
227* .. External Functions ..
228 REAL SLAMCH
229 EXTERNAL slamch
230* ..
231* .. External Subroutines ..
232 EXTERNAL clacpy, clartg, classq, crot
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC abs, conjg, max, real, sqrt
236* ..
237* .. Executable Statements ..
238*
239 info = 0
240*
241* Quick return if possible
242*
243 IF( n.LE.1 )
244 \$ RETURN
245*
246 m = ldst
247 weak = .false.
248 strong = .false.
249*
250* Make a local copy of selected block in (A, B)
251*
252 CALL clacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
253 CALL clacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
254*
255* Compute the threshold for testing the acceptance of swapping.
256*
257 eps = slamch( 'P' )
258 smlnum = slamch( 'S' ) / eps
259 scale = real( czero )
260 sum = real( cone )
261 CALL clacpy( 'Full', m, m, s, ldst, work, m )
262 CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
263 CALL classq( m*m, work, 1, scale, sum )
264 sa = scale*sqrt( sum )
265 scale = dble( czero )
266 sum = dble( cone )
267 CALL classq( m*m, work(m*m+1), 1, scale, sum )
268 sb = scale*sqrt( sum )
269*
270* THRES has been changed from
271* THRESH = MAX( TEN*EPS*SA, SMLNUM )
272* to
273* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
274* on 04/01/10.
275* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
276* Jim Demmel and Guillaume Revy. See forum post 1783.
277*
278 thresha = max( twenty*eps*sa, smlnum )
279 threshb = max( twenty*eps*sb, smlnum )
280*
281* Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
282* using Givens rotations and perform the swap tentatively.
283*
284 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
285 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
286 sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
287 sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
288 CALL clartg( g, f, cz, sz, cdum )
289 sz = -sz
290 CALL crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
291 CALL crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
292 IF( sa.GE.sb ) THEN
293 CALL clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
294 ELSE
295 CALL clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
296 END IF
297 CALL crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
298 CALL crot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, cq, sq )
299*
300* Weak stability test: |S21| <= O(EPS F-norm((A)))
301* and |T21| <= O(EPS F-norm((B)))
302*
303 weak = abs( s( 2, 1 ) ).LE.thresha .AND.
304 \$ abs( t( 2, 1 ) ).LE.threshb
305 IF( .NOT.weak )
306 \$ GO TO 20
307*
308 IF( wands ) THEN
309*
310* Strong stability test:
311* F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
312*
313 CALL clacpy( 'Full', m, m, s, ldst, work, m )
314 CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
315 CALL crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
316 CALL crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
317 CALL crot( 2, work, 2, work( 2 ), 2, cq, -sq )
318 CALL crot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
319 DO 10 i = 1, 2
320 work( i ) = work( i ) - a( j1+i-1, j1 )
321 work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
322 work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
323 work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
324 10 CONTINUE
325 scale = dble( czero )
326 sum = dble( cone )
327 CALL classq( m*m, work, 1, scale, sum )
328 sa = scale*sqrt( sum )
329 scale = dble( czero )
330 sum = dble( cone )
331 CALL classq( m*m, work(m*m+1), 1, scale, sum )
332 sb = scale*sqrt( sum )
333 strong = sa.LE.thresha .AND. sb.LE.threshb
334 IF( .NOT.strong )
335 \$ GO TO 20
336 END IF
337*
338* If the swap is accepted ("weakly" and "strongly"), apply the
339* equivalence transformations to the original matrix pair (A,B)
340*
341 CALL crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
342 CALL crot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz, conjg( sz ) )
343 CALL crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
344 CALL crot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
345*
346* Set N1 by N2 (2,1) blocks to 0
347*
348 a( j1+1, j1 ) = czero
349 b( j1+1, j1 ) = czero
350*
351* Accumulate transformations into Q and Z if requested.
352*
353 IF( wantz )
354 \$ CALL crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
355 IF( wantq )
356 \$ CALL crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
357*
358* Exit with INFO = 0 if swap was successfully performed.
359*
360 RETURN
361*
362* Exit with INFO = 1 if swap was rejected.
363*
364 20 CONTINUE
365 info = 1
366 RETURN
367*
368* End of CTGEX2
369*
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:116
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:137
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: crot.f:103
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: