SRC/ctgex2.f

   1 *> \brief \b CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CTGEX2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgex2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgex2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgex2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
  22 *                          LDZ, J1, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            WANTQ, WANTZ
  26 *       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  30 *      $                   Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
  40 *> in an upper triangular matrix pair (A, B) by an unitary equivalence
  41 *> transformation.
  42 *>
  43 *> (A, B) must be in generalized Schur canonical form, that is, A and
  44 *> B are both upper triangular.
  45 *>
  46 *> Optionally, the matrices Q and Z of generalized Schur vectors are
  47 *> updated.
  48 *>
  49 *>        Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
  50 *>        Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
  51 *>
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] WANTQ
  58 *> \verbatim
  59 *>          WANTQ is LOGICAL
  60 *>          .TRUE. : update the left transformation matrix Q;
  61 *>          .FALSE.: do not update Q.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] WANTZ
  65 *> \verbatim
  66 *>          WANTZ is LOGICAL
  67 *>          .TRUE. : update the right transformation matrix Z;
  68 *>          .FALSE.: do not update Z.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] N
  72 *> \verbatim
  73 *>          N is INTEGER
  74 *>          The order of the matrices A and B. N >= 0.
  75 *> \endverbatim
  76 *>
  77 *> \param[in,out] A
  78 *> \verbatim
  79 *>          A is COMPLEX arrays, dimensions (LDA,N)
  80 *>          On entry, the matrix A in the pair (A, B).
  81 *>          On exit, the updated matrix A.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] LDA
  85 *> \verbatim
  86 *>          LDA is INTEGER
  87 *>          The leading dimension of the array A. LDA >= max(1,N).
  88 *> \endverbatim
  89 *>
  90 *> \param[in,out] B
  91 *> \verbatim
  92 *>          B is COMPLEX arrays, dimensions (LDB,N)
  93 *>          On entry, the matrix B in the pair (A, B).
  94 *>          On exit, the updated matrix B.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] LDB
  98 *> \verbatim
  99 *>          LDB is INTEGER
 100 *>          The leading dimension of the array B. LDB >= max(1,N).
 101 *> \endverbatim
 102 *>
 103 *> \param[in,out] Q
 104 *> \verbatim
 105 *>          Q is COMPLEX array, dimension (LDZ,N)
 106 *>          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
 107 *>          the updated matrix Q.
 108 *>          Not referenced if WANTQ = .FALSE..
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDQ
 112 *> \verbatim
 113 *>          LDQ is INTEGER
 114 *>          The leading dimension of the array Q. LDQ >= 1;
 115 *>          If WANTQ = .TRUE., LDQ >= N.
 116 *> \endverbatim
 117 *>
 118 *> \param[in,out] Z
 119 *> \verbatim
 120 *>          Z is COMPLEX array, dimension (LDZ,N)
 121 *>          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
 122 *>          the updated matrix Z.
 123 *>          Not referenced if WANTZ = .FALSE..
 124 *> \endverbatim
 125 *>
 126 *> \param[in] LDZ
 127 *> \verbatim
 128 *>          LDZ is INTEGER
 129 *>          The leading dimension of the array Z. LDZ >= 1;
 130 *>          If WANTZ = .TRUE., LDZ >= N.
 131 *> \endverbatim
 132 *>
 133 *> \param[in] J1
 134 *> \verbatim
 135 *>          J1 is INTEGER
 136 *>          The index to the first block (A11, B11).
 137 *> \endverbatim
 138 *>
 139 *> \param[out] INFO
 140 *> \verbatim
 141 *>          INFO is INTEGER
 142 *>           =0:  Successful exit.
 143 *>           =1:  The transformed matrix pair (A, B) would be too far
 144 *>                from generalized Schur form; the problem is ill-
 145 *>                conditioned.
 146 *> \endverbatim
 147 *
 148 *  Authors:
 149 *  ========
 150 *
 151 *> \author Univ. of Tennessee
 152 *> \author Univ. of California Berkeley
 153 *> \author Univ. of Colorado Denver
 154 *> \author NAG Ltd.
 155 *
 156 *> \date September 2012
 157 *
 158 *> \ingroup complexGEauxiliary
 159 *
 160 *> \par Further Details:
 161 *  =====================
 162 *>
 163 *>  In the current code both weak and strong stability tests are
 164 *>  performed. The user can omit the strong stability test by changing
 165 *>  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
 166 *>  details.
 167 *
 168 *> \par Contributors:
 169 *  ==================
 170 *>
 171 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 172 *>     Umea University, S-901 87 Umea, Sweden.
 173 *
 174 *> \par References:
 175 *  ================
 176 *>
 177 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
 178 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
 179 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
 180 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
 181 *> \n
 182 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
 183 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
 184 *>      Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
 185 *>      Department of Computing Science, Umea University, S-901 87 Umea,
 186 *>      Sweden, 1994. Also as LAPACK Working Note 87. To appear in
 187 *>      Numerical Algorithms, 1996.
 188 *>
 189 *  =====================================================================
 190       SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 191      $                   LDZ, J1, INFO )
 192 *
 193 *  -- LAPACK auxiliary routine (version 3.4.2) --
 194 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 195 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 196 *     September 2012
 197 *
 198 *     .. Scalar Arguments ..
 199       LOGICAL            WANTQ, WANTZ
 200       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
 201 *     ..
 202 *     .. Array Arguments ..
 203       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 204      $                   Z( LDZ, * )
 205 *     ..
 206 *
 207 *  =====================================================================
 208 *
 209 *     .. Parameters ..
 210       COMPLEX            CZERO, CONE
 211       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 212      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 213       REAL               TWENTY
 214       PARAMETER          ( TWENTY = 2.0E+1 )
 215       INTEGER            LDST
 216       PARAMETER          ( LDST = 2 )
 217       LOGICAL            WANDS
 218       PARAMETER          ( WANDS = .TRUE. )
 219 *     ..
 220 *     .. Local Scalars ..
 221       LOGICAL            STRONG, WEAK
 222       INTEGER            I, M
 223       REAL               CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
 224      $                   THRESH, WS
 225       COMPLEX            CDUM, F, G, SQ, SZ
 226 *     ..
 227 *     .. Local Arrays ..
 228       COMPLEX            S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
 229 *     ..
 230 *     .. External Functions ..
 231       REAL               SLAMCH
 232       EXTERNAL           SLAMCH
 233 *     ..
 234 *     .. External Subroutines ..
 235       EXTERNAL           CLACPY, CLARTG, CLASSQ, CROT
 236 *     ..
 237 *     .. Intrinsic Functions ..
 238       INTRINSIC          ABS, CONJG, MAX, REAL, SQRT
 239 *     ..
 240 *     .. Executable Statements ..
 241 *
 242       INFO = 0
 243 *
 244 *     Quick return if possible
 245 *
 246       IF( N.LE.1 )
 247      $   RETURN
 248 *
 249       M = LDST
 250       WEAK = .FALSE.
 251       STRONG = .FALSE.
 252 *
 253 *     Make a local copy of selected block in (A, B)
 254 *
 255       CALL CLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
 256       CALL CLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
 257 *
 258 *     Compute the threshold for testing the acceptance of swapping.
 259 *
 260       EPS = SLAMCH( 'P' )
 261       SMLNUM = SLAMCH( 'S' ) / EPS
 262       SCALE = REAL( CZERO )
 263       SUM = REAL( CONE )
 264       CALL CLACPY( 'Full', M, M, S, LDST, WORK, M )
 265       CALL CLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
 266       CALL CLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
 267       SA = SCALE*SQRT( SUM )
 268 *
 269 *     THRES has been changed from
 270 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
 271 *     to
 272 *        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
 273 *     on 04/01/10.
 274 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
 275 *     Jim Demmel and Guillaume Revy. See forum post 1783.
 276 *
 277       THRESH = MAX( TWENTY*EPS*SA, 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 ) )
 285       SB = 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| + |T21| <= O(EPS F-norm((S, T)))
 299 *
 300       WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
 301       WEAK = WS.LE.THRESH
 302       IF( .NOT.WEAK )
 303      $   GO TO 20
 304 *
 305       IF( WANDS ) THEN
 306 *
 307 *        Strong stability test:
 308 *           F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
 309 *
 310          CALL CLACPY( 'Full', M, M, S, LDST, WORK, M )
 311          CALL CLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
 312          CALL CROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -CONJG( SZ ) )
 313          CALL CROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -CONJG( SZ ) )
 314          CALL CROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
 315          CALL CROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
 316          DO 10 I = 1, 2
 317             WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
 318             WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
 319             WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
 320             WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
 321    10    CONTINUE
 322          SCALE = REAL( CZERO )
 323          SUM = REAL( CONE )
 324          CALL CLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
 325          SS = SCALE*SQRT( SUM )
 326          STRONG = SS.LE.THRESH
 327          IF( .NOT.STRONG )
 328      $      GO TO 20
 329       END IF
 330 *
 331 *     If the swap is accepted ("weakly" and "strongly"), apply the
 332 *     equivalence transformations to the original matrix pair (A,B)
 333 *
 334       CALL CROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ, CONJG( SZ ) )
 335       CALL CROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ, CONJG( SZ ) )
 336       CALL CROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
 337       CALL CROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
 338 *
 339 *     Set  N1 by N2 (2,1) blocks to 0
 340 *
 341       A( J1+1, J1 ) = CZERO
 342       B( J1+1, J1 ) = CZERO
 343 *
 344 *     Accumulate transformations into Q and Z if requested.
 345 *
 346       IF( WANTZ )
 347      $   CALL CROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ, CONJG( SZ ) )
 348       IF( WANTQ )
 349      $   CALL CROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ, CONJG( SQ ) )
 350 *
 351 *     Exit with INFO = 0 if swap was successfully performed.
 352 *
 353       RETURN
 354 *
 355 *     Exit with INFO = 1 if swap was rejected.
 356 *
 357    20 CONTINUE
 358       INFO = 1
 359       RETURN
 360 *
 361 *     End of CTGEX2
 362 *
 363       END