ctgex2.f

Go to the documentation of this file.

*> \brief \b CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CTGEX2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgex2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgex2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgex2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
*                          LDZ, J1, INFO )
*
*       .. Scalar Arguments ..
*       LOGICAL            WANTQ, WANTZ
*       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
*       ..
*       .. Array Arguments ..
*       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
*      $                   Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] WANTQ
*> \verbatim
*>          WANTQ is LOGICAL
*>          .TRUE. : update the left transformation matrix Q;
*>          .FALSE.: do not update Q.
*> \endverbatim
*>
*> \param[in] WANTZ
*> \verbatim
*>          WANTZ is LOGICAL
*>          .TRUE. : update the right transformation matrix Z;
*>          .FALSE.: do not update Z.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A and B. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX array, dimension (LDA,N)
*>          On entry, the matrix A in the pair (A, B).
*>          On exit, the updated matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is COMPLEX array, dimension (LDB,N)
*>          On entry, the matrix B in the pair (A, B).
*>          On exit, the updated matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          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..
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q. LDQ >= 1;
*>          If WANTQ = .TRUE., LDQ >= N.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          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..
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z. LDZ >= 1;
*>          If WANTZ = .TRUE., LDZ >= N.
*> \endverbatim
*>
*> \param[in] J1
*> \verbatim
*>          J1 is INTEGER
*>          The index to the first block (A11, B11).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tgex2
*
*> \par 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.
*
*> \par Contributors:
*  ==================
*>
*>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*>     Umea University, S-901 87 Umea, Sweden.
*
*> \par 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.
*> \n
*>  [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.
*>
*  =====================================================================
      SUBROUTINE ctgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, J1, INFO )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            WANTQ, WANTZ
      INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $                   z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            CZERO, CONE
      parameter( czero = ( 0.0e+0, 0.0e+0 ),
     $                   cone = ( 1.0e+0, 0.0e+0 ) )
      REAL               TWENTY
      parameter( twenty = 2.0e+1 )
      INTEGER            LDST
      parameter( ldst = 2 )
      LOGICAL            WANDS
      parameter( wands = .true. )
*     ..
*     .. Local Scalars ..
      LOGICAL            STRONG, WEAK
      INTEGER            I, M
      REAL               CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
     $                   thresha, threshb
      COMPLEX            CDUM, F, G, SQ, SZ
*     ..
*     .. Local Arrays ..
      COMPLEX            S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
*     ..
*     .. External Functions ..
      REAL               SLAMCH
      EXTERNAL           slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           clacpy, clartg, classq, crot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, conjg, max, real, sqrt
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.LE.1 )
     $   RETURN
*
      m = ldst
      weak = .false.
      strong = .false.
*
*     Make a local copy of selected block in (A, B)
*
      CALL clacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
      CALL clacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
*
*     Compute the threshold for testing the acceptance of swapping.
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' ) / eps
      scale = real( czero )
      sum = real( cone )
      CALL clacpy( 'Full', m, m, s, ldst, work, m )
      CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
      CALL classq( m*m, work, 1, scale, sum )
      sa = scale*sqrt( sum )
      scale = dble( czero )
      sum = dble( cone )
      CALL classq( m*m, work(m*m+1), 1, scale, sum )
      sb = scale*sqrt( sum )
*
*     THRES has been changed from
*        THRESH = MAX( TEN*EPS*SA, SMLNUM )
*     to
*        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
*     on 04/01/10.
*     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
*     Jim Demmel and Guillaume Revy. See forum post 1783.
*
      thresha = max( twenty*eps*sa, smlnum )
      threshb = max( twenty*eps*sb, smlnum )
*
*     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
*     using Givens rotations and perform the swap tentatively.
*
      f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
      g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
      sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
      sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
      CALL clartg( g, f, cz, sz, cdum )
      sz = -sz
      CALL crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
      CALL crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
      IF( sa.GE.sb ) THEN
         CALL clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
      ELSE
         CALL clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
      END IF
      CALL crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
      CALL crot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, cq, sq )
*
*     Weak stability test: |S21| <= O(EPS F-norm((A)))
*                          and  |T21| <= O(EPS F-norm((B)))
*
      weak = abs( s( 2, 1 ) ).LE.thresha .AND. 
     $ abs( t( 2, 1 ) ).LE.threshb
      IF( .NOT.weak )
     $   GO TO 20
*
      IF( wands ) THEN
*
*        Strong stability test:
*           F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
*
         CALL clacpy( 'Full', m, m, s, ldst, work, m )
         CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
         CALL crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
         CALL crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
         CALL crot( 2, work, 2, work( 2 ), 2, cq, -sq )
         CALL crot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
         DO 10 i = 1, 2
            work( i ) = work( i ) - a( j1+i-1, j1 )
            work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
            work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
            work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
   10    CONTINUE
         scale = dble( czero )
         sum = dble( cone )
         CALL classq( m*m, work, 1, scale, sum )
         sa = scale*sqrt( sum )
         scale = dble( czero )
         sum = dble( cone )
         CALL classq( m*m, work(m*m+1), 1, scale, sum )
         sb = scale*sqrt( sum )
         strong = sa.LE.thresha .AND. sb.LE.threshb
         IF( .NOT.strong )
     $      GO TO 20
      END IF
*
*     If the swap is accepted ("weakly" and "strongly"), apply the
*     equivalence transformations to the original matrix pair (A,B)
*
      CALL crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
      CALL crot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz, conjg( sz ) )
      CALL crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
      CALL crot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
*
*     Set  N1 by N2 (2,1) blocks to 0
*
      a( j1+1, j1 ) = czero
      b( j1+1, j1 ) = czero
*
*     Accumulate transformations into Q and Z if requested.
*
      IF( wantz )
     $   CALL crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
      IF( wantq )
     $   CALL crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
*
*     Exit with INFO = 0 if swap was successfully performed.
*
      RETURN
*
*     Exit with INFO = 1 if swap was rejected.
*
   20 CONTINUE
      info = 1
      RETURN
*
*     End of CTGEX2
*
      END