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 December 2016 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.7.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* December 2016 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 698