LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgex2.f
Go to the documentation of this file.
00001 *> \brief \b STGEX2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGEX2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgex2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgex2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgex2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00022 *                          LDZ, J1, N1, N2, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       LOGICAL            WANTQ, WANTZ
00026 *       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00030 *      $                   WORK( * ), Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> STGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
00040 *> of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
00041 *> (A, B) by an orthogonal equivalence transformation.
00042 *>
00043 *> (A, B) must be in generalized real Schur canonical form (as returned
00044 *> by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
00045 *> diagonal blocks. B is upper triangular.
00046 *>
00047 *> Optionally, the matrices Q and Z of generalized Schur vectors are
00048 *> updated.
00049 *>
00050 *>        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
00051 *>        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
00052 *>
00053 *> \endverbatim
00054 *
00055 *  Arguments:
00056 *  ==========
00057 *
00058 *> \param[in] WANTQ
00059 *> \verbatim
00060 *>          WANTQ is LOGICAL
00061 *>          .TRUE. : update the left transformation matrix Q;
00062 *>          .FALSE.: do not update Q.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] WANTZ
00066 *> \verbatim
00067 *>          WANTZ is LOGICAL
00068 *>          .TRUE. : update the right transformation matrix Z;
00069 *>          .FALSE.: do not update Z.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] N
00073 *> \verbatim
00074 *>          N is INTEGER
00075 *>          The order of the matrices A and B. N >= 0.
00076 *> \endverbatim
00077 *>
00078 *> \param[in,out] A
00079 *> \verbatim
00080 *>          A is REAL arrays, dimensions (LDA,N)
00081 *>          On entry, the matrix A in the pair (A, B).
00082 *>          On exit, the updated matrix A.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] LDA
00086 *> \verbatim
00087 *>          LDA is INTEGER
00088 *>          The leading dimension of the array A. LDA >= max(1,N).
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] B
00092 *> \verbatim
00093 *>          B is REAL arrays, dimensions (LDB,N)
00094 *>          On entry, the matrix B in the pair (A, B).
00095 *>          On exit, the updated matrix B.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] LDB
00099 *> \verbatim
00100 *>          LDB is INTEGER
00101 *>          The leading dimension of the array B. LDB >= max(1,N).
00102 *> \endverbatim
00103 *>
00104 *> \param[in,out] Q
00105 *> \verbatim
00106 *>          Q is REAL array, dimension (LDZ,N)
00107 *>          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
00108 *>          On exit, the updated matrix Q.
00109 *>          Not referenced if WANTQ = .FALSE..
00110 *> \endverbatim
00111 *>
00112 *> \param[in] LDQ
00113 *> \verbatim
00114 *>          LDQ is INTEGER
00115 *>          The leading dimension of the array Q. LDQ >= 1.
00116 *>          If WANTQ = .TRUE., LDQ >= N.
00117 *> \endverbatim
00118 *>
00119 *> \param[in,out] Z
00120 *> \verbatim
00121 *>          Z is REAL array, dimension (LDZ,N)
00122 *>          On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
00123 *>          On exit, the updated matrix Z.
00124 *>          Not referenced if WANTZ = .FALSE..
00125 *> \endverbatim
00126 *>
00127 *> \param[in] LDZ
00128 *> \verbatim
00129 *>          LDZ is INTEGER
00130 *>          The leading dimension of the array Z. LDZ >= 1.
00131 *>          If WANTZ = .TRUE., LDZ >= N.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] J1
00135 *> \verbatim
00136 *>          J1 is INTEGER
00137 *>          The index to the first block (A11, B11). 1 <= J1 <= N.
00138 *> \endverbatim
00139 *>
00140 *> \param[in] N1
00141 *> \verbatim
00142 *>          N1 is INTEGER
00143 *>          The order of the first block (A11, B11). N1 = 0, 1 or 2.
00144 *> \endverbatim
00145 *>
00146 *> \param[in] N2
00147 *> \verbatim
00148 *>          N2 is INTEGER
00149 *>          The order of the second block (A22, B22). N2 = 0, 1 or 2.
00150 *> \endverbatim
00151 *>
00152 *> \param[out] WORK
00153 *> \verbatim
00154 *>          WORK is REAL array, dimension (MAX(1,LWORK)).
00155 *> \endverbatim
00156 *>
00157 *> \param[in] LWORK
00158 *> \verbatim
00159 *>          LWORK is INTEGER
00160 *>          The dimension of the array WORK.
00161 *>          LWORK >=  MAX( N*(N2+N1), (N2+N1)*(N2+N1)*2 )
00162 *> \endverbatim
00163 *>
00164 *> \param[out] INFO
00165 *> \verbatim
00166 *>          INFO is INTEGER
00167 *>            =0: Successful exit
00168 *>            >0: If INFO = 1, the transformed matrix (A, B) would be
00169 *>                too far from generalized Schur form; the blocks are
00170 *>                not swapped and (A, B) and (Q, Z) are unchanged.
00171 *>                The problem of swapping is too ill-conditioned.
00172 *>            <0: If INFO = -16: LWORK is too small. Appropriate value
00173 *>                for LWORK is returned in WORK(1).
00174 *> \endverbatim
00175 *
00176 *  Authors:
00177 *  ========
00178 *
00179 *> \author Univ. of Tennessee 
00180 *> \author Univ. of California Berkeley 
00181 *> \author Univ. of Colorado Denver 
00182 *> \author NAG Ltd. 
00183 *
00184 *> \date November 2011
00185 *
00186 *> \ingroup realGEauxiliary
00187 *
00188 *> \par Further Details:
00189 *  =====================
00190 *>
00191 *>  In the current code both weak and strong stability tests are
00192 *>  performed. The user can omit the strong stability test by changing
00193 *>  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
00194 *>  details.
00195 *
00196 *> \par Contributors:
00197 *  ==================
00198 *>
00199 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00200 *>     Umea University, S-901 87 Umea, Sweden.
00201 *
00202 *> \par References:
00203 *  ================
00204 *>
00205 *> \verbatim
00206 *>
00207 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00208 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00209 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00210 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00211 *>
00212 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00213 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00214 *>      Estimation: Theory, Algorithms and Software,
00215 *>      Report UMINF - 94.04, Department of Computing Science, Umea
00216 *>      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
00217 *>      Note 87. To appear in Numerical Algorithms, 1996.
00218 *> \endverbatim
00219 *>
00220 *  =====================================================================
00221       SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00222      $                   LDZ, J1, N1, N2, WORK, LWORK, INFO )
00223 *
00224 *  -- LAPACK auxiliary routine (version 3.4.0) --
00225 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00226 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00227 *     November 2011
00228 *
00229 *     .. Scalar Arguments ..
00230       LOGICAL            WANTQ, WANTZ
00231       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
00232 *     ..
00233 *     .. Array Arguments ..
00234       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00235      $                   WORK( * ), Z( LDZ, * )
00236 *     ..
00237 *
00238 *  =====================================================================
00239 *  Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
00240 *  loops. Sven Hammarling, 1/5/02.
00241 *
00242 *     .. Parameters ..
00243       REAL               ZERO, ONE
00244       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00245       REAL               TWENTY
00246       PARAMETER          ( TWENTY = 2.0E+01 )
00247       INTEGER            LDST
00248       PARAMETER          ( LDST = 4 )
00249       LOGICAL            WANDS
00250       PARAMETER          ( WANDS = .TRUE. )
00251 *     ..
00252 *     .. Local Scalars ..
00253       LOGICAL            STRONG, WEAK
00254       INTEGER            I, IDUM, LINFO, M
00255       REAL               BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
00256      $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
00257 *     ..
00258 *     .. Local Arrays ..
00259       INTEGER            IWORK( LDST )
00260       REAL               AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
00261      $                   IRCOP( LDST, LDST ), LI( LDST, LDST ),
00262      $                   LICOP( LDST, LDST ), S( LDST, LDST ),
00263      $                   SCPY( LDST, LDST ), T( LDST, LDST ),
00264      $                   TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
00265 *     ..
00266 *     .. External Functions ..
00267       REAL               SLAMCH
00268       EXTERNAL           SLAMCH
00269 *     ..
00270 *     .. External Subroutines ..
00271       EXTERNAL           SGEMM, SGEQR2, SGERQ2, SLACPY, SLAGV2, SLARTG,
00272      $                   SLASET, SLASSQ, SORG2R, SORGR2, SORM2R, SORMR2,
00273      $                   SROT, SSCAL, STGSY2
00274 *     ..
00275 *     .. Intrinsic Functions ..
00276       INTRINSIC          ABS, MAX, SQRT
00277 *     ..
00278 *     .. Executable Statements ..
00279 *
00280       INFO = 0
00281 *
00282 *     Quick return if possible
00283 *
00284       IF( N.LE.1 .OR. N1.LE.0 .OR. N2.LE.0 )
00285      $   RETURN
00286       IF( N1.GT.N .OR. ( J1+N1 ).GT.N )
00287      $   RETURN
00288       M = N1 + N2
00289       IF( LWORK.LT.MAX( N*M, M*M*2 ) ) THEN
00290          INFO = -16
00291          WORK( 1 ) = MAX( N*M, M*M*2 )
00292          RETURN
00293       END IF
00294 *
00295       WEAK = .FALSE.
00296       STRONG = .FALSE.
00297 *
00298 *     Make a local copy of selected block
00299 *
00300       CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, LI, LDST )
00301       CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, IR, LDST )
00302       CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
00303       CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
00304 *
00305 *     Compute threshold for testing acceptance of swapping.
00306 *
00307       EPS = SLAMCH( 'P' )
00308       SMLNUM = SLAMCH( 'S' ) / EPS
00309       DSCALE = ZERO
00310       DSUM = ONE
00311       CALL SLACPY( 'Full', M, M, S, LDST, WORK, M )
00312       CALL SLASSQ( M*M, WORK, 1, DSCALE, DSUM )
00313       CALL SLACPY( 'Full', M, M, T, LDST, WORK, M )
00314       CALL SLASSQ( M*M, WORK, 1, DSCALE, DSUM )
00315       DNORM = DSCALE*SQRT( DSUM )
00316 *
00317 *     THRES has been changed from 
00318 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
00319 *     to
00320 *        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
00321 *     on 04/01/10.
00322 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
00323 *     Jim Demmel and Guillaume Revy. See forum post 1783.
00324 *
00325       THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )
00326 *
00327       IF( M.EQ.2 ) THEN
00328 *
00329 *        CASE 1: Swap 1-by-1 and 1-by-1 blocks.
00330 *
00331 *        Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
00332 *        using Givens rotations and perform the swap tentatively.
00333 *
00334          F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
00335          G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
00336          SB = ABS( T( 2, 2 ) )
00337          SA = ABS( S( 2, 2 ) )
00338          CALL SLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
00339          IR( 2, 1 ) = -IR( 1, 2 )
00340          IR( 2, 2 ) = IR( 1, 1 )
00341          CALL SROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ),
00342      $              IR( 2, 1 ) )
00343          CALL SROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ),
00344      $              IR( 2, 1 ) )
00345          IF( SA.GE.SB ) THEN
00346             CALL SLARTG( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
00347      $                   DDUM )
00348          ELSE
00349             CALL SLARTG( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
00350      $                   DDUM )
00351          END IF
00352          CALL SROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ),
00353      $              LI( 2, 1 ) )
00354          CALL SROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, LI( 1, 1 ),
00355      $              LI( 2, 1 ) )
00356          LI( 2, 2 ) = LI( 1, 1 )
00357          LI( 1, 2 ) = -LI( 2, 1 )
00358 *
00359 *        Weak stability test:
00360 *           |S21| + |T21| <= O(EPS * F-norm((S, T)))
00361 *
00362          WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
00363          WEAK = WS.LE.THRESH
00364          IF( .NOT.WEAK )
00365      $      GO TO 70
00366 *
00367          IF( WANDS ) THEN
00368 *
00369 *           Strong stability test:
00370 *           F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A, B)))
00371 *
00372             CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
00373      $                   M )
00374             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
00375      $                  WORK, M )
00376             CALL SGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
00377      $                  WORK( M*M+1 ), M )
00378             DSCALE = ZERO
00379             DSUM = ONE
00380             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
00381 *
00382             CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
00383      $                   M )
00384             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
00385      $                  WORK, M )
00386             CALL SGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
00387      $                  WORK( M*M+1 ), M )
00388             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
00389             SS = DSCALE*SQRT( DSUM )
00390             STRONG = SS.LE.THRESH
00391             IF( .NOT.STRONG )
00392      $         GO TO 70
00393          END IF
00394 *
00395 *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
00396 *               (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
00397 *
00398          CALL SROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ),
00399      $              IR( 2, 1 ) )
00400          CALL SROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ),
00401      $              IR( 2, 1 ) )
00402          CALL SROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA,
00403      $              LI( 1, 1 ), LI( 2, 1 ) )
00404          CALL SROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB,
00405      $              LI( 1, 1 ), LI( 2, 1 ) )
00406 *
00407 *        Set  N1-by-N2 (2,1) - blocks to ZERO.
00408 *
00409          A( J1+1, J1 ) = ZERO
00410          B( J1+1, J1 ) = ZERO
00411 *
00412 *        Accumulate transformations into Q and Z if requested.
00413 *
00414          IF( WANTZ )
00415      $      CALL SROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
00416      $                 IR( 2, 1 ) )
00417          IF( WANTQ )
00418      $      CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, LI( 1, 1 ),
00419      $                 LI( 2, 1 ) )
00420 *
00421 *        Exit with INFO = 0 if swap was successfully performed.
00422 *
00423          RETURN
00424 *
00425       ELSE
00426 *
00427 *        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
00428 *                and 2-by-2 blocks.
00429 *
00430 *        Solve the generalized Sylvester equation
00431 *                 S11 * R - L * S22 = SCALE * S12
00432 *                 T11 * R - L * T22 = SCALE * T12
00433 *        for R and L. Solutions in LI and IR.
00434 *
00435          CALL SLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
00436          CALL SLACPY( 'Full', N1, N2, S( 1, N1+1 ), LDST,
00437      $                IR( N2+1, N1+1 ), LDST )
00438          CALL STGSY2( 'N', 0, N1, N2, S, LDST, S( N1+1, N1+1 ), LDST,
00439      $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
00440      $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
00441      $                LINFO )
00442 *
00443 *        Compute orthogonal matrix QL:
00444 *
00445 *                    QL**T * LI = [ TL ]
00446 *                                 [ 0  ]
00447 *        where
00448 *                    LI =  [      -L              ]
00449 *                          [ SCALE * identity(N2) ]
00450 *
00451          DO 10 I = 1, N2
00452             CALL SSCAL( N1, -ONE, LI( 1, I ), 1 )
00453             LI( N1+I, I ) = SCALE
00454    10    CONTINUE
00455          CALL SGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO )
00456          IF( LINFO.NE.0 )
00457      $      GO TO 70
00458          CALL SORG2R( M, M, N2, LI, LDST, TAUL, WORK, LINFO )
00459          IF( LINFO.NE.0 )
00460      $      GO TO 70
00461 *
00462 *        Compute orthogonal matrix RQ:
00463 *
00464 *                    IR * RQ**T =   [ 0  TR],
00465 *
00466 *         where IR = [ SCALE * identity(N1), R ]
00467 *
00468          DO 20 I = 1, N1
00469             IR( N2+I, I ) = SCALE
00470    20    CONTINUE
00471          CALL SGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
00472          IF( LINFO.NE.0 )
00473      $      GO TO 70
00474          CALL SORGR2( M, M, N1, IR, LDST, TAUR, WORK, LINFO )
00475          IF( LINFO.NE.0 )
00476      $      GO TO 70
00477 *
00478 *        Perform the swapping tentatively:
00479 *
00480          CALL SGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
00481      $               WORK, M )
00482          CALL SGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, S,
00483      $               LDST )
00484          CALL SGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
00485      $               WORK, M )
00486          CALL SGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, T,
00487      $               LDST )
00488          CALL SLACPY( 'F', M, M, S, LDST, SCPY, LDST )
00489          CALL SLACPY( 'F', M, M, T, LDST, TCPY, LDST )
00490          CALL SLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
00491          CALL SLACPY( 'F', M, M, LI, LDST, LICOP, LDST )
00492 *
00493 *        Triangularize the B-part by an RQ factorization.
00494 *        Apply transformation (from left) to A-part, giving S.
00495 *
00496          CALL SGERQ2( M, M, T, LDST, TAUR, WORK, LINFO )
00497          IF( LINFO.NE.0 )
00498      $      GO TO 70
00499          CALL SORMR2( 'R', 'T', M, M, M, T, LDST, TAUR, S, LDST, WORK,
00500      $                LINFO )
00501          IF( LINFO.NE.0 )
00502      $      GO TO 70
00503          CALL SORMR2( 'L', 'N', M, M, M, T, LDST, TAUR, IR, LDST, WORK,
00504      $                LINFO )
00505          IF( LINFO.NE.0 )
00506      $      GO TO 70
00507 *
00508 *        Compute F-norm(S21) in BRQA21. (T21 is 0.)
00509 *
00510          DSCALE = ZERO
00511          DSUM = ONE
00512          DO 30 I = 1, N2
00513             CALL SLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
00514    30    CONTINUE
00515          BRQA21 = DSCALE*SQRT( DSUM )
00516 *
00517 *        Triangularize the B-part by a QR factorization.
00518 *        Apply transformation (from right) to A-part, giving S.
00519 *
00520          CALL SGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO )
00521          IF( LINFO.NE.0 )
00522      $      GO TO 70
00523          CALL SORM2R( 'L', 'T', M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
00524      $                WORK, INFO )
00525          CALL SORM2R( 'R', 'N', M, M, M, TCPY, LDST, TAUL, LICOP, LDST,
00526      $                WORK, INFO )
00527          IF( LINFO.NE.0 )
00528      $      GO TO 70
00529 *
00530 *        Compute F-norm(S21) in BQRA21. (T21 is 0.)
00531 *
00532          DSCALE = ZERO
00533          DSUM = ONE
00534          DO 40 I = 1, N2
00535             CALL SLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
00536    40    CONTINUE
00537          BQRA21 = DSCALE*SQRT( DSUM )
00538 *
00539 *        Decide which method to use.
00540 *          Weak stability test:
00541 *             F-norm(S21) <= O(EPS * F-norm((S, T)))
00542 *
00543          IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
00544             CALL SLACPY( 'F', M, M, SCPY, LDST, S, LDST )
00545             CALL SLACPY( 'F', M, M, TCPY, LDST, T, LDST )
00546             CALL SLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
00547             CALL SLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
00548          ELSE IF( BRQA21.GE.THRESH ) THEN
00549             GO TO 70
00550          END IF
00551 *
00552 *        Set lower triangle of B-part to zero
00553 *
00554          CALL SLASET( 'Lower', M-1, M-1, ZERO, ZERO, T(2,1), LDST )
00555 *
00556          IF( WANDS ) THEN
00557 *
00558 *           Strong stability test:
00559 *              F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
00560 *
00561             CALL SLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
00562      $                   M )
00563             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
00564      $                  WORK, M )
00565             CALL SGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
00566      $                  WORK( M*M+1 ), M )
00567             DSCALE = ZERO
00568             DSUM = ONE
00569             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
00570 *
00571             CALL SLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
00572      $                   M )
00573             CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
00574      $                  WORK, M )
00575             CALL SGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
00576      $                  WORK( M*M+1 ), M )
00577             CALL SLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
00578             SS = DSCALE*SQRT( DSUM )
00579             STRONG = ( SS.LE.THRESH )
00580             IF( .NOT.STRONG )
00581      $         GO TO 70
00582 *
00583          END IF
00584 *
00585 *        If the swap is accepted ("weakly" and "strongly"), apply the
00586 *        transformations and set N1-by-N2 (2,1)-block to zero.
00587 *
00588          CALL SLASET( 'Full', N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
00589 *
00590 *        copy back M-by-M diagonal block starting at index J1 of (A, B)
00591 *
00592          CALL SLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA )
00593          CALL SLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB )
00594          CALL SLASET( 'Full', LDST, LDST, ZERO, ZERO, T, LDST )
00595 *
00596 *        Standardize existing 2-by-2 blocks.
00597 *
00598          DO 50 I = 1, M*M
00599             WORK(I) = ZERO
00600    50    CONTINUE
00601          WORK( 1 ) = ONE
00602          T( 1, 1 ) = ONE
00603          IDUM = LWORK - M*M - 2
00604          IF( N2.GT.1 ) THEN
00605             CALL SLAGV2( A( J1, J1 ), LDA, B( J1, J1 ), LDB, AR, AI, BE,
00606      $                   WORK( 1 ), WORK( 2 ), T( 1, 1 ), T( 2, 1 ) )
00607             WORK( M+1 ) = -WORK( 2 )
00608             WORK( M+2 ) = WORK( 1 )
00609             T( N2, N2 ) = T( 1, 1 )
00610             T( 1, 2 ) = -T( 2, 1 )
00611          END IF
00612          WORK( M*M ) = ONE
00613          T( M, M ) = ONE
00614 *
00615          IF( N1.GT.1 ) THEN
00616             CALL SLAGV2( A( J1+N2, J1+N2 ), LDA, B( J1+N2, J1+N2 ), LDB,
00617      $                   TAUR, TAUL, WORK( M*M+1 ), WORK( N2*M+N2+1 ),
00618      $                   WORK( N2*M+N2+2 ), T( N2+1, N2+1 ),
00619      $                   T( M, M-1 ) )
00620             WORK( M*M ) = WORK( N2*M+N2+1 )
00621             WORK( M*M-1 ) = -WORK( N2*M+N2+2 )
00622             T( M, M ) = T( N2+1, N2+1 )
00623             T( M-1, M ) = -T( M, M-1 )
00624          END IF
00625          CALL SGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, A( J1, J1+N2 ),
00626      $               LDA, ZERO, WORK( M*M+1 ), N2 )
00627          CALL SLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ),
00628      $                LDA )
00629          CALL SGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, B( J1, J1+N2 ),
00630      $               LDB, ZERO, WORK( M*M+1 ), N2 )
00631          CALL SLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ),
00632      $                LDB )
00633          CALL SGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, WORK, M, ZERO,
00634      $               WORK( M*M+1 ), M )
00635          CALL SLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST )
00636          CALL SGEMM( 'N', 'N', N2, N1, N1, ONE, A( J1, J1+N2 ), LDA,
00637      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
00638          CALL SLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA )
00639          CALL SGEMM( 'N', 'N', N2, N1, N1, ONE, B( J1, J1+N2 ), LDB,
00640      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
00641          CALL SLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB )
00642          CALL SGEMM( 'T', 'N', M, M, M, ONE, IR, LDST, T, LDST, ZERO,
00643      $               WORK, M )
00644          CALL SLACPY( 'Full', M, M, WORK, M, IR, LDST )
00645 *
00646 *        Accumulate transformations into Q and Z if requested.
00647 *
00648          IF( WANTQ ) THEN
00649             CALL SGEMM( 'N', 'N', N, M, M, ONE, Q( 1, J1 ), LDQ, LI,
00650      $                  LDST, ZERO, WORK, N )
00651             CALL SLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ )
00652 *
00653          END IF
00654 *
00655          IF( WANTZ ) THEN
00656             CALL SGEMM( 'N', 'N', N, M, M, ONE, Z( 1, J1 ), LDZ, IR,
00657      $                  LDST, ZERO, WORK, N )
00658             CALL SLACPY( 'Full', N, M, WORK, N, Z( 1, J1 ), LDZ )
00659 *
00660          END IF
00661 *
00662 *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
00663 *                (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
00664 *
00665          I = J1 + M
00666          IF( I.LE.N ) THEN
00667             CALL SGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
00668      $                  A( J1, I ), LDA, ZERO, WORK, M )
00669             CALL SLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA )
00670             CALL SGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
00671      $                  B( J1, I ), LDB, ZERO, WORK, M )
00672             CALL SLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB )
00673          END IF
00674          I = J1 - 1
00675          IF( I.GT.0 ) THEN
00676             CALL SGEMM( 'N', 'N', I, M, M, ONE, A( 1, J1 ), LDA, IR,
00677      $                  LDST, ZERO, WORK, I )
00678             CALL SLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA )
00679             CALL SGEMM( 'N', 'N', I, M, M, ONE, B( 1, J1 ), LDB, IR,
00680      $                  LDST, ZERO, WORK, I )
00681             CALL SLACPY( 'Full', I, M, WORK, I, B( 1, J1 ), LDB )
00682          END IF
00683 *
00684 *        Exit with INFO = 0 if swap was successfully performed.
00685 *
00686          RETURN
00687 *
00688       END IF
00689 *
00690 *     Exit with INFO = 1 if swap was rejected.
00691 *
00692    70 CONTINUE
00693 *
00694       INFO = 1
00695       RETURN
00696 *
00697 *     End of STGEX2
00698 *
00699       END
 All Files Functions