LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztgex2.f
Go to the documentation of this file.
00001 *> \brief \b ZTGEX2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZTGEX2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgex2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgex2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgex2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00022 *                          LDZ, J1, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       LOGICAL            WANTQ, WANTZ
00026 *       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00030 *      $                   Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
00040 *> in an upper triangular matrix pair (A, B) by an unitary equivalence
00041 *> transformation.
00042 *>
00043 *> (A, B) must be in generalized Schur canonical form, that is, A and
00044 *> B are both upper triangular.
00045 *>
00046 *> Optionally, the matrices Q and Z of generalized Schur vectors are
00047 *> updated.
00048 *>
00049 *>        Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
00050 *>        Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
00051 *>
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] WANTQ
00058 *> \verbatim
00059 *>          WANTQ is LOGICAL
00060 *>          .TRUE. : update the left transformation matrix Q;
00061 *>          .FALSE.: do not update Q.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] WANTZ
00065 *> \verbatim
00066 *>          WANTZ is LOGICAL
00067 *>          .TRUE. : update the right transformation matrix Z;
00068 *>          .FALSE.: do not update Z.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] N
00072 *> \verbatim
00073 *>          N is INTEGER
00074 *>          The order of the matrices A and B. N >= 0.
00075 *> \endverbatim
00076 *>
00077 *> \param[in,out] A
00078 *> \verbatim
00079 *>          A is COMPLEX*16 arrays, dimensions (LDA,N)
00080 *>          On entry, the matrix A in the pair (A, B).
00081 *>          On exit, the updated matrix A.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] LDA
00085 *> \verbatim
00086 *>          LDA is INTEGER
00087 *>          The leading dimension of the array A. LDA >= max(1,N).
00088 *> \endverbatim
00089 *>
00090 *> \param[in,out] B
00091 *> \verbatim
00092 *>          B is COMPLEX*16 arrays, dimensions (LDB,N)
00093 *>          On entry, the matrix B in the pair (A, B).
00094 *>          On exit, the updated matrix B.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] LDB
00098 *> \verbatim
00099 *>          LDB is INTEGER
00100 *>          The leading dimension of the array B. LDB >= max(1,N).
00101 *> \endverbatim
00102 *>
00103 *> \param[in,out] Q
00104 *> \verbatim
00105 *>          Q is COMPLEX*16 array, dimension (LDZ,N)
00106 *>          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
00107 *>          the updated matrix Q.
00108 *>          Not referenced if WANTQ = .FALSE..
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDQ
00112 *> \verbatim
00113 *>          LDQ is INTEGER
00114 *>          The leading dimension of the array Q. LDQ >= 1;
00115 *>          If WANTQ = .TRUE., LDQ >= N.
00116 *> \endverbatim
00117 *>
00118 *> \param[in,out] Z
00119 *> \verbatim
00120 *>          Z is COMPLEX*16 array, dimension (LDZ,N)
00121 *>          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
00122 *>          the updated matrix Z.
00123 *>          Not referenced if WANTZ = .FALSE..
00124 *> \endverbatim
00125 *>
00126 *> \param[in] LDZ
00127 *> \verbatim
00128 *>          LDZ is INTEGER
00129 *>          The leading dimension of the array Z. LDZ >= 1;
00130 *>          If WANTZ = .TRUE., LDZ >= N.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] J1
00134 *> \verbatim
00135 *>          J1 is INTEGER
00136 *>          The index to the first block (A11, B11).
00137 *> \endverbatim
00138 *>
00139 *> \param[out] INFO
00140 *> \verbatim
00141 *>          INFO is INTEGER
00142 *>           =0:  Successful exit.
00143 *>           =1:  The transformed matrix pair (A, B) would be too far
00144 *>                from generalized Schur form; the problem is ill-
00145 *>                conditioned. 
00146 *> \endverbatim
00147 *
00148 *  Authors:
00149 *  ========
00150 *
00151 *> \author Univ. of Tennessee 
00152 *> \author Univ. of California Berkeley 
00153 *> \author Univ. of Colorado Denver 
00154 *> \author NAG Ltd. 
00155 *
00156 *> \date November 2011
00157 *
00158 *> \ingroup complex16GEauxiliary
00159 *
00160 *> \par Further Details:
00161 *  =====================
00162 *>
00163 *>  In the current code both weak and strong stability tests are
00164 *>  performed. The user can omit the strong stability test by changing
00165 *>  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
00166 *>  details.
00167 *
00168 *> \par Contributors:
00169 *  ==================
00170 *>
00171 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00172 *>     Umea University, S-901 87 Umea, Sweden.
00173 *
00174 *> \par References:
00175 *  ================
00176 *>
00177 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00178 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00179 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00180 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00181 *> \n
00182 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00183 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00184 *>      Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
00185 *>      Department of Computing Science, Umea University, S-901 87 Umea,
00186 *>      Sweden, 1994. Also as LAPACK Working Note 87. To appear in
00187 *>      Numerical Algorithms, 1996.
00188 *>
00189 *  =====================================================================
00190       SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00191      $                   LDZ, J1, INFO )
00192 *
00193 *  -- LAPACK auxiliary routine (version 3.4.0) --
00194 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00195 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00196 *     November 2011
00197 *
00198 *     .. Scalar Arguments ..
00199       LOGICAL            WANTQ, WANTZ
00200       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, N
00201 *     ..
00202 *     .. Array Arguments ..
00203       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00204      $                   Z( LDZ, * )
00205 *     ..
00206 *
00207 *  =====================================================================
00208 *
00209 *     .. Parameters ..
00210       COMPLEX*16         CZERO, CONE
00211       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00212      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00213       DOUBLE PRECISION   TWENTY
00214       PARAMETER          ( TWENTY = 2.0D+1 )
00215       INTEGER            LDST
00216       PARAMETER          ( LDST = 2 )
00217       LOGICAL            WANDS
00218       PARAMETER          ( WANDS = .TRUE. )
00219 *     ..
00220 *     .. Local Scalars ..
00221       LOGICAL            DTRONG, WEAK
00222       INTEGER            I, M
00223       DOUBLE PRECISION   CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
00224      $                   THRESH, WS
00225       COMPLEX*16         CDUM, F, G, SQ, SZ
00226 *     ..
00227 *     .. Local Arrays ..
00228       COMPLEX*16         S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
00229 *     ..
00230 *     .. External Functions ..
00231       DOUBLE PRECISION   DLAMCH
00232       EXTERNAL           DLAMCH
00233 *     ..
00234 *     .. External Subroutines ..
00235       EXTERNAL           ZLACPY, ZLARTG, ZLASSQ, ZROT
00236 *     ..
00237 *     .. Intrinsic Functions ..
00238       INTRINSIC          ABS, DBLE, DCONJG, MAX, SQRT
00239 *     ..
00240 *     .. Executable Statements ..
00241 *
00242       INFO = 0
00243 *
00244 *     Quick return if possible
00245 *
00246       IF( N.LE.1 )
00247      $   RETURN
00248 *
00249       M = LDST
00250       WEAK = .FALSE.
00251       DTRONG = .FALSE.
00252 *
00253 *     Make a local copy of selected block in (A, B)
00254 *
00255       CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
00256       CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
00257 *
00258 *     Compute the threshold for testing the acceptance of swapping.
00259 *
00260       EPS = DLAMCH( 'P' )
00261       SMLNUM = DLAMCH( 'S' ) / EPS
00262       SCALE = DBLE( CZERO )
00263       SUM = DBLE( CONE )
00264       CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
00265       CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
00266       CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
00267       SA = SCALE*SQRT( SUM )
00268 *
00269 *     THRES has been changed from 
00270 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
00271 *     to
00272 *        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
00273 *     on 04/01/10.
00274 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
00275 *     Jim Demmel and Guillaume Revy. See forum post 1783.
00276 *
00277       THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
00278 *
00279 *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
00280 *     using Givens rotations and perform the swap tentatively.
00281 *
00282       F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
00283       G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
00284       SA = ABS( S( 2, 2 ) )
00285       SB = ABS( T( 2, 2 ) )
00286       CALL ZLARTG( G, F, CZ, SZ, CDUM )
00287       SZ = -SZ
00288       CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
00289       CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
00290       IF( SA.GE.SB ) THEN
00291          CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
00292       ELSE
00293          CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
00294       END IF
00295       CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
00296       CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
00297 *
00298 *     Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
00299 *
00300       WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
00301       WEAK = WS.LE.THRESH
00302       IF( .NOT.WEAK )
00303      $   GO TO 20
00304 *
00305       IF( WANDS ) THEN
00306 *
00307 *        Strong stability test:
00308 *           F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
00309 *
00310          CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
00311          CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
00312          CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
00313          CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
00314          CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
00315          CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
00316          DO 10 I = 1, 2
00317             WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
00318             WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
00319             WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
00320             WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
00321    10    CONTINUE
00322          SCALE = DBLE( CZERO )
00323          SUM = DBLE( CONE )
00324          CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
00325          SS = SCALE*SQRT( SUM )
00326          DTRONG = SS.LE.THRESH
00327          IF( .NOT.DTRONG )
00328      $      GO TO 20
00329       END IF
00330 *
00331 *     If the swap is accepted ("weakly" and "strongly"), apply the
00332 *     equivalence transformations to the original matrix pair (A,B)
00333 *
00334       CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
00335      $           DCONJG( SZ ) )
00336       CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
00337      $           DCONJG( SZ ) )
00338       CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
00339       CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
00340 *
00341 *     Set  N1 by N2 (2,1) blocks to 0
00342 *
00343       A( J1+1, J1 ) = CZERO
00344       B( J1+1, J1 ) = CZERO
00345 *
00346 *     Accumulate transformations into Q and Z if requested.
00347 *
00348       IF( WANTZ )
00349      $   CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
00350      $              DCONJG( SZ ) )
00351       IF( WANTQ )
00352      $   CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
00353      $              DCONJG( SQ ) )
00354 *
00355 *     Exit with INFO = 0 if swap was successfully performed.
00356 *
00357       RETURN
00358 *
00359 *     Exit with INFO = 1 if swap was rejected.
00360 *
00361    20 CONTINUE
00362       INFO = 1
00363       RETURN
00364 *
00365 *     End of ZTGEX2
00366 *
00367       END
 All Files Functions