LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
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.

Download CTGEX2 + dependencies [TGZ] [ZIP] [TXT]

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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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. [2] for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
[1] 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.
[2] 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 186 of file ctgex2.f.

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