LAPACK: SRC/ztgex2.f Source File

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