SRC/dtgex2.f

   1 *> \brief \b DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal 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 DTGEX2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
  22 *                          LDZ, J1, N1, N2, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            WANTQ, WANTZ
  26 *       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  30 *      $                   WORK( * ), Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
  40 *> of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
  41 *> (A, B) by an orthogonal equivalence transformation.
  42 *>
  43 *> (A, B) must be in generalized real Schur canonical form (as returned
  44 *> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
  45 *> diagonal blocks. B is upper triangular.
  46 *>
  47 *> Optionally, the matrices Q and Z of generalized Schur vectors are
  48 *> updated.
  49 *>
  50 *>        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
  51 *>        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
  52 *>
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] WANTQ
  59 *> \verbatim
  60 *>          WANTQ is LOGICAL
  61 *>          .TRUE. : update the left transformation matrix Q;
  62 *>          .FALSE.: do not update Q.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] WANTZ
  66 *> \verbatim
  67 *>          WANTZ is LOGICAL
  68 *>          .TRUE. : update the right transformation matrix Z;
  69 *>          .FALSE.: do not update Z.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] N
  73 *> \verbatim
  74 *>          N is INTEGER
  75 *>          The order of the matrices A and B. N >= 0.
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] A
  79 *> \verbatim
  80 *>          A is DOUBLE PRECISION array, dimensions (LDA,N)
  81 *>          On entry, the matrix A in the pair (A, B).
  82 *>          On exit, the updated matrix A.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] LDA
  86 *> \verbatim
  87 *>          LDA is INTEGER
  88 *>          The leading dimension of the array A. LDA >= max(1,N).
  89 *> \endverbatim
  90 *>
  91 *> \param[in,out] B
  92 *> \verbatim
  93 *>          B is DOUBLE PRECISION array, dimensions (LDB,N)
  94 *>          On entry, the matrix B in the pair (A, B).
  95 *>          On exit, the updated matrix B.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] LDB
  99 *> \verbatim
 100 *>          LDB is INTEGER
 101 *>          The leading dimension of the array B. LDB >= max(1,N).
 102 *> \endverbatim
 103 *>
 104 *> \param[in,out] Q
 105 *> \verbatim
 106 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
 107 *>          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
 108 *>          On exit, the updated matrix Q.
 109 *>          Not referenced if WANTQ = .FALSE..
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDQ
 113 *> \verbatim
 114 *>          LDQ is INTEGER
 115 *>          The leading dimension of the array Q. LDQ >= 1.
 116 *>          If WANTQ = .TRUE., LDQ >= N.
 117 *> \endverbatim
 118 *>
 119 *> \param[in,out] Z
 120 *> \verbatim
 121 *>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
 122 *>          On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
 123 *>          On exit, the updated matrix Z.
 124 *>          Not referenced if WANTZ = .FALSE..
 125 *> \endverbatim
 126 *>
 127 *> \param[in] LDZ
 128 *> \verbatim
 129 *>          LDZ is INTEGER
 130 *>          The leading dimension of the array Z. LDZ >= 1.
 131 *>          If WANTZ = .TRUE., LDZ >= N.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] J1
 135 *> \verbatim
 136 *>          J1 is INTEGER
 137 *>          The index to the first block (A11, B11). 1 <= J1 <= N.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] N1
 141 *> \verbatim
 142 *>          N1 is INTEGER
 143 *>          The order of the first block (A11, B11). N1 = 0, 1 or 2.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] N2
 147 *> \verbatim
 148 *>          N2 is INTEGER
 149 *>          The order of the second block (A22, B22). N2 = 0, 1 or 2.
 150 *> \endverbatim
 151 *>
 152 *> \param[out] WORK
 153 *> \verbatim
 154 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
 155 *> \endverbatim
 156 *>
 157 *> \param[in] LWORK
 158 *> \verbatim
 159 *>          LWORK is INTEGER
 160 *>          The dimension of the array WORK.
 161 *>          LWORK >=  MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
 162 *> \endverbatim
 163 *>
 164 *> \param[out] INFO
 165 *> \verbatim
 166 *>          INFO is INTEGER
 167 *>            =0: Successful exit
 168 *>            >0: If INFO = 1, the transformed matrix (A, B) would be
 169 *>                too far from generalized Schur form; the blocks are
 170 *>                not swapped and (A, B) and (Q, Z) are unchanged.
 171 *>                The problem of swapping is too ill-conditioned.
 172 *>            <0: If INFO = -16: LWORK is too small. Appropriate value
 173 *>                for LWORK is returned in WORK(1).
 174 *> \endverbatim
 175 *
 176 *  Authors:
 177 *  ========
 178 *
 179 *> \author Univ. of Tennessee
 180 *> \author Univ. of California Berkeley
 181 *> \author Univ. of Colorado Denver
 182 *> \author NAG Ltd.
 183 *
 184 *> \date November 2015
 185 *
 186 *> \ingroup doubleGEauxiliary
 187 *
 188 *> \par Further Details:
 189 *  =====================
 190 *>
 191 *>  In the current code both weak and strong stability tests are
 192 *>  performed. The user can omit the strong stability test by changing
 193 *>  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
 194 *>  details.
 195 *
 196 *> \par Contributors:
 197 *  ==================
 198 *>
 199 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 200 *>     Umea University, S-901 87 Umea, Sweden.
 201 *
 202 *> \par References:
 203 *  ================
 204 *>
 205 *> \verbatim
 206 *>
 207 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
 208 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
 209 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
 210 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
 211 *>
 212 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
 213 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
 214 *>      Estimation: Theory, Algorithms and Software,
 215 *>      Report UMINF - 94.04, Department of Computing Science, Umea
 216 *>      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
 217 *>      Note 87. To appear in Numerical Algorithms, 1996.
 218 *> \endverbatim
 219 *>
 220 *  =====================================================================
 221       SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 222      $                   LDZ, J1, N1, N2, WORK, LWORK, INFO )
 223 *
 224 *  -- LAPACK auxiliary routine (version 3.6.0) --
 225 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 226 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 227 *     November 2015
 228 *
 229 *     .. Scalar Arguments ..
 230       LOGICAL            WANTQ, WANTZ
 231       INTEGER            INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
 232 *     ..
 233 *     .. Array Arguments ..
 234       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 235      $                   WORK( * ), Z( LDZ, * )
 236 *     ..
 237 *
 238 *  =====================================================================
 239 *  Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
 240 *  loops. Sven Hammarling, 1/5/02.
 241 *
 242 *     .. Parameters ..
 243       DOUBLE PRECISION   ZERO, ONE
 244       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 245       DOUBLE PRECISION   TWENTY
 246       PARAMETER          ( TWENTY = 2.0D+01 )
 247       INTEGER            LDST
 248       PARAMETER          ( LDST = 4 )
 249       LOGICAL            WANDS
 250       PARAMETER          ( WANDS = .TRUE. )
 251 *     ..
 252 *     .. Local Scalars ..
 253       LOGICAL            DTRONG, WEAK
 254       INTEGER            I, IDUM, LINFO, M
 255       DOUBLE PRECISION   BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
 256      $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
 257 *     ..
 258 *     .. Local Arrays ..
 259       INTEGER            IWORK( LDST )
 260       DOUBLE PRECISION   AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
 261      $                   IRCOP( LDST, LDST ), LI( LDST, LDST ),
 262      $                   LICOP( LDST, LDST ), S( LDST, LDST ),
 263      $                   SCPY( LDST, LDST ), T( LDST, LDST ),
 264      $                   TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
 265 *     ..
 266 *     .. External Functions ..
 267       DOUBLE PRECISION   DLAMCH
 268       EXTERNAL           DLAMCH
 269 *     ..
 270 *     .. External Subroutines ..
 271       EXTERNAL           DGEMM, DGEQR2, DGERQ2, DLACPY, DLAGV2, DLARTG,
 272      $                   DLASET, DLASSQ, DORG2R, DORGR2, DORM2R, DORMR2,
 273      $                   DROT, DSCAL, DTGSY2
 274 *     ..
 275 *     .. Intrinsic Functions ..
 276       INTRINSIC          ABS, MAX, SQRT
 277 *     ..
 278 *     .. Executable Statements ..
 279 *
 280       INFO = 0
 281 *
 282 *     Quick return if possible
 283 *
 284       IF( N.LE.1 .OR. N1.LE.0 .OR. N2.LE.0 )
 285      $   RETURN
 286       IF( N1.GT.N .OR. ( J1+N1 ).GT.N )
 287      $   RETURN
 288       M = N1 + N2
 289       IF( LWORK.LT.MAX( 1, N*M, M*M*2 ) ) THEN
 290          INFO = -16
 291          WORK( 1 ) = MAX( 1, N*M, M*M*2 )
 292          RETURN
 293       END IF
 294 *
 295       WEAK = .FALSE.
 296       DTRONG = .FALSE.
 297 *
 298 *     Make a local copy of selected block
 299 *
 300       CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, LI, LDST )
 301       CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, IR, LDST )
 302       CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
 303       CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
 304 *
 305 *     Compute threshold for testing acceptance of swapping.
 306 *
 307       EPS = DLAMCH( 'P' )
 308       SMLNUM = DLAMCH( 'S' ) / EPS
 309       DSCALE = ZERO
 310       DSUM = ONE
 311       CALL DLACPY( 'Full', M, M, S, LDST, WORK, M )
 312       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
 313       CALL DLACPY( 'Full', M, M, T, LDST, WORK, M )
 314       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
 315       DNORM = DSCALE*SQRT( DSUM )
 316 *
 317 *     THRES has been changed from
 318 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
 319 *     to
 320 *        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
 321 *     on 04/01/10.
 322 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
 323 *     Jim Demmel and Guillaume Revy. See forum post 1783.
 324 *
 325       THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )
 326 *
 327       IF( M.EQ.2 ) THEN
 328 *
 329 *        CASE 1: Swap 1-by-1 and 1-by-1 blocks.
 330 *
 331 *        Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
 332 *        using Givens rotations and perform the swap tentatively.
 333 *
 334          F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
 335          G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
 336          SB = ABS( T( 2, 2 ) )
 337          SA = ABS( S( 2, 2 ) )
 338          CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
 339          IR( 2, 1 ) = -IR( 1, 2 )
 340          IR( 2, 2 ) = IR( 1, 1 )
 341          CALL DROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ),
 342      $              IR( 2, 1 ) )
 343          CALL DROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ),
 344      $              IR( 2, 1 ) )
 345          IF( SA.GE.SB ) THEN
 346             CALL DLARTG( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
 347      $                   DDUM )
 348          ELSE
 349             CALL DLARTG( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
 350      $                   DDUM )
 351          END IF
 352          CALL DROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ),
 353      $              LI( 2, 1 ) )
 354          CALL DROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, LI( 1, 1 ),
 355      $              LI( 2, 1 ) )
 356          LI( 2, 2 ) = LI( 1, 1 )
 357          LI( 1, 2 ) = -LI( 2, 1 )
 358 *
 359 *        Weak stability test:
 360 *           |S21| + |T21| <= O(EPS * F-norm((S, T)))
 361 *
 362          WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
 363          WEAK = WS.LE.THRESH
 364          IF( .NOT.WEAK )
 365      $      GO TO 70
 366 *
 367          IF( WANDS ) THEN
 368 *
 369 *           Strong stability test:
 370 *             F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))
 371 *
 372             CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
 373      $                   M )
 374             CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
 375      $                  WORK, M )
 376             CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
 377      $                  WORK( M*M+1 ), M )
 378             DSCALE = ZERO
 379             DSUM = ONE
 380             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
 381 *
 382             CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
 383      $                   M )
 384             CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
 385      $                  WORK, M )
 386             CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
 387      $                  WORK( M*M+1 ), M )
 388             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
 389             SS = DSCALE*SQRT( DSUM )
 390             DTRONG = SS.LE.THRESH
 391             IF( .NOT.DTRONG )
 392      $         GO TO 70
 393          END IF
 394 *
 395 *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
 396 *               (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
 397 *
 398          CALL DROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ),
 399      $              IR( 2, 1 ) )
 400          CALL DROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ),
 401      $              IR( 2, 1 ) )
 402          CALL DROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA,
 403      $              LI( 1, 1 ), LI( 2, 1 ) )
 404          CALL DROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB,
 405      $              LI( 1, 1 ), LI( 2, 1 ) )
 406 *
 407 *        Set  N1-by-N2 (2,1) - blocks to ZERO.
 408 *
 409          A( J1+1, J1 ) = ZERO
 410          B( J1+1, J1 ) = ZERO
 411 *
 412 *        Accumulate transformations into Q and Z if requested.
 413 *
 414          IF( WANTZ )
 415      $      CALL DROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
 416      $                 IR( 2, 1 ) )
 417          IF( WANTQ )
 418      $      CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, LI( 1, 1 ),
 419      $                 LI( 2, 1 ) )
 420 *
 421 *        Exit with INFO = 0 if swap was successfully performed.
 422 *
 423          RETURN
 424 *
 425       ELSE
 426 *
 427 *        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
 428 *                and 2-by-2 blocks.
 429 *
 430 *        Solve the generalized Sylvester equation
 431 *                 S11 * R - L * S22 = SCALE * S12
 432 *                 T11 * R - L * T22 = SCALE * T12
 433 *        for R and L. Solutions in LI and IR.
 434 *
 435          CALL DLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
 436          CALL DLACPY( 'Full', N1, N2, S( 1, N1+1 ), LDST,
 437      $                IR( N2+1, N1+1 ), LDST )
 438          CALL DTGSY2( 'N', 0, N1, N2, S, LDST, S( N1+1, N1+1 ), LDST,
 439      $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
 440      $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
 441      $                LINFO )
 442 *
 443 *        Compute orthogonal matrix QL:
 444 *
 445 *                    QL**T * LI = [ TL ]
 446 *                                 [ 0  ]
 447 *        where
 448 *                    LI =  [      -L              ]
 449 *                          [ SCALE * identity(N2) ]
 450 *
 451          DO 10 I = 1, N2
 452             CALL DSCAL( N1, -ONE, LI( 1, I ), 1 )
 453             LI( N1+I, I ) = SCALE
 454    10    CONTINUE
 455          CALL DGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO )
 456          IF( LINFO.NE.0 )
 457      $      GO TO 70
 458          CALL DORG2R( M, M, N2, LI, LDST, TAUL, WORK, LINFO )
 459          IF( LINFO.NE.0 )
 460      $      GO TO 70
 461 *
 462 *        Compute orthogonal matrix RQ:
 463 *
 464 *                    IR * RQ**T =   [ 0  TR],
 465 *
 466 *         where IR = [ SCALE * identity(N1), R ]
 467 *
 468          DO 20 I = 1, N1
 469             IR( N2+I, I ) = SCALE
 470    20    CONTINUE
 471          CALL DGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
 472          IF( LINFO.NE.0 )
 473      $      GO TO 70
 474          CALL DORGR2( M, M, N1, IR, LDST, TAUR, WORK, LINFO )
 475          IF( LINFO.NE.0 )
 476      $      GO TO 70
 477 *
 478 *        Perform the swapping tentatively:
 479 *
 480          CALL DGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
 481      $               WORK, M )
 482          CALL DGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, S,
 483      $               LDST )
 484          CALL DGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
 485      $               WORK, M )
 486          CALL DGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, T,
 487      $               LDST )
 488          CALL DLACPY( 'F', M, M, S, LDST, SCPY, LDST )
 489          CALL DLACPY( 'F', M, M, T, LDST, TCPY, LDST )
 490          CALL DLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
 491          CALL DLACPY( 'F', M, M, LI, LDST, LICOP, LDST )
 492 *
 493 *        Triangularize the B-part by an RQ factorization.
 494 *        Apply transformation (from left) to A-part, giving S.
 495 *
 496          CALL DGERQ2( M, M, T, LDST, TAUR, WORK, LINFO )
 497          IF( LINFO.NE.0 )
 498      $      GO TO 70
 499          CALL DORMR2( 'R', 'T', M, M, M, T, LDST, TAUR, S, LDST, WORK,
 500      $                LINFO )
 501          IF( LINFO.NE.0 )
 502      $      GO TO 70
 503          CALL DORMR2( 'L', 'N', M, M, M, T, LDST, TAUR, IR, LDST, WORK,
 504      $                LINFO )
 505          IF( LINFO.NE.0 )
 506      $      GO TO 70
 507 *
 508 *        Compute F-norm(S21) in BRQA21. (T21 is 0.)
 509 *
 510          DSCALE = ZERO
 511          DSUM = ONE
 512          DO 30 I = 1, N2
 513             CALL DLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
 514    30    CONTINUE
 515          BRQA21 = DSCALE*SQRT( DSUM )
 516 *
 517 *        Triangularize the B-part by a QR factorization.
 518 *        Apply transformation (from right) to A-part, giving S.
 519 *
 520          CALL DGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO )
 521          IF( LINFO.NE.0 )
 522      $      GO TO 70
 523          CALL DORM2R( 'L', 'T', M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
 524      $                WORK, INFO )
 525          CALL DORM2R( 'R', 'N', M, M, M, TCPY, LDST, TAUL, LICOP, LDST,
 526      $                WORK, INFO )
 527          IF( LINFO.NE.0 )
 528      $      GO TO 70
 529 *
 530 *        Compute F-norm(S21) in BQRA21. (T21 is 0.)
 531 *
 532          DSCALE = ZERO
 533          DSUM = ONE
 534          DO 40 I = 1, N2
 535             CALL DLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
 536    40    CONTINUE
 537          BQRA21 = DSCALE*SQRT( DSUM )
 538 *
 539 *        Decide which method to use.
 540 *          Weak stability test:
 541 *             F-norm(S21) <= O(EPS * F-norm((S, T)))
 542 *
 543          IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
 544             CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST )
 545             CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST )
 546             CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
 547             CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
 548          ELSE IF( BRQA21.GE.THRESH ) THEN
 549             GO TO 70
 550          END IF
 551 *
 552 *        Set lower triangle of B-part to zero
 553 *
 554          CALL DLASET( 'Lower', M-1, M-1, ZERO, ZERO, T(2,1), LDST )
 555 *
 556          IF( WANDS ) THEN
 557 *
 558 *           Strong stability test:
 559 *              F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
 560 *
 561             CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
 562      $                   M )
 563             CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
 564      $                  WORK, M )
 565             CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
 566      $                  WORK( M*M+1 ), M )
 567             DSCALE = ZERO
 568             DSUM = ONE
 569             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
 570 *
 571             CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
 572      $                   M )
 573             CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
 574      $                  WORK, M )
 575             CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
 576      $                  WORK( M*M+1 ), M )
 577             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
 578             SS = DSCALE*SQRT( DSUM )
 579             DTRONG = ( SS.LE.THRESH )
 580             IF( .NOT.DTRONG )
 581      $         GO TO 70
 582 *
 583          END IF
 584 *
 585 *        If the swap is accepted ("weakly" and "strongly"), apply the
 586 *        transformations and set N1-by-N2 (2,1)-block to zero.
 587 *
 588          CALL DLASET( 'Full', N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
 589 *
 590 *        copy back M-by-M diagonal block starting at index J1 of (A, B)
 591 *
 592          CALL DLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA )
 593          CALL DLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB )
 594          CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, T, LDST )
 595 *
 596 *        Standardize existing 2-by-2 blocks.
 597 *
 598          CALL DLASET( 'Full', M, M, ZERO, ZERO, WORK, M )
 599          WORK( 1 ) = ONE
 600          T( 1, 1 ) = ONE
 601          IDUM = LWORK - M*M - 2
 602          IF( N2.GT.1 ) THEN
 603             CALL DLAGV2( A( J1, J1 ), LDA, B( J1, J1 ), LDB, AR, AI, BE,
 604      $                   WORK( 1 ), WORK( 2 ), T( 1, 1 ), T( 2, 1 ) )
 605             WORK( M+1 ) = -WORK( 2 )
 606             WORK( M+2 ) = WORK( 1 )
 607             T( N2, N2 ) = T( 1, 1 )
 608             T( 1, 2 ) = -T( 2, 1 )
 609          END IF
 610          WORK( M*M ) = ONE
 611          T( M, M ) = ONE
 612 *
 613          IF( N1.GT.1 ) THEN
 614             CALL DLAGV2( A( J1+N2, J1+N2 ), LDA, B( J1+N2, J1+N2 ), LDB,
 615      $                   TAUR, TAUL, WORK( M*M+1 ), WORK( N2*M+N2+1 ),
 616      $                   WORK( N2*M+N2+2 ), T( N2+1, N2+1 ),
 617      $                   T( M, M-1 ) )
 618             WORK( M*M ) = WORK( N2*M+N2+1 )
 619             WORK( M*M-1 ) = -WORK( N2*M+N2+2 )
 620             T( M, M ) = T( N2+1, N2+1 )
 621             T( M-1, M ) = -T( M, M-1 )
 622          END IF
 623          CALL DGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, A( J1, J1+N2 ),
 624      $               LDA, ZERO, WORK( M*M+1 ), N2 )
 625          CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ),
 626      $                LDA )
 627          CALL DGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, B( J1, J1+N2 ),
 628      $               LDB, ZERO, WORK( M*M+1 ), N2 )
 629          CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ),
 630      $                LDB )
 631          CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, WORK, M, ZERO,
 632      $               WORK( M*M+1 ), M )
 633          CALL DLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST )
 634          CALL DGEMM( 'N', 'N', N2, N1, N1, ONE, A( J1, J1+N2 ), LDA,
 635      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
 636          CALL DLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA )
 637          CALL DGEMM( 'N', 'N', N2, N1, N1, ONE, B( J1, J1+N2 ), LDB,
 638      $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
 639          CALL DLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB )
 640          CALL DGEMM( 'T', 'N', M, M, M, ONE, IR, LDST, T, LDST, ZERO,
 641      $               WORK, M )
 642          CALL DLACPY( 'Full', M, M, WORK, M, IR, LDST )
 643 *
 644 *        Accumulate transformations into Q and Z if requested.
 645 *
 646          IF( WANTQ ) THEN
 647             CALL DGEMM( 'N', 'N', N, M, M, ONE, Q( 1, J1 ), LDQ, LI,
 648      $                  LDST, ZERO, WORK, N )
 649             CALL DLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ )
 650 *
 651          END IF
 652 *
 653          IF( WANTZ ) THEN
 654             CALL DGEMM( 'N', 'N', N, M, M, ONE, Z( 1, J1 ), LDZ, IR,
 655      $                  LDST, ZERO, WORK, N )
 656             CALL DLACPY( 'Full', N, M, WORK, N, Z( 1, J1 ), LDZ )
 657 *
 658          END IF
 659 *
 660 *        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
 661 *                (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
 662 *
 663          I = J1 + M
 664          IF( I.LE.N ) THEN
 665             CALL DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
 666      $                  A( J1, I ), LDA, ZERO, WORK, M )
 667             CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA )
 668             CALL DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
 669      $                  B( J1, I ), LDB, ZERO, WORK, M )
 670             CALL DLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB )
 671          END IF
 672          I = J1 - 1
 673          IF( I.GT.0 ) THEN
 674             CALL DGEMM( 'N', 'N', I, M, M, ONE, A( 1, J1 ), LDA, IR,
 675      $                  LDST, ZERO, WORK, I )
 676             CALL DLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA )
 677             CALL DGEMM( 'N', 'N', I, M, M, ONE, B( 1, J1 ), LDB, IR,
 678      $                  LDST, ZERO, WORK, I )
 679             CALL DLACPY( 'Full', I, M, WORK, I, B( 1, J1 ), LDB )
 680          END IF
 681 *
 682 *        Exit with INFO = 0 if swap was successfully performed.
 683 *
 684          RETURN
 685 *
 686       END IF
 687 *
 688 *     Exit with INFO = 1 if swap was rejected.
 689 *
 690    70 CONTINUE
 691 *
 692       INFO = 1
 693       RETURN
 694 *
 695 *     End of DTGEX2
 696 *
 697       END