LAPACK 3.12.0
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 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*
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:103
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:124
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:103
Here is the call graph for this function:
Here is the caller graph for this function: