![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTGEX2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTGEX2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTGEX2( 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 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00030 * $ WORK( * ), Z( LDZ, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DTGEX2 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 DGGES), 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 DOUBLE PRECISION array, 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 DOUBLE PRECISION array, 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 DOUBLE PRECISION array, dimension (LDQ,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 DOUBLE PRECISION 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 DOUBLE PRECISION 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( 1, 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 doubleGEauxiliary 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 DTGEX2( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00235 $ WORK( * ), Z( LDZ, * ) 00236 * .. 00237 * 00238 * ===================================================================== 00239 * Replaced various illegal calls to DCOPY by calls to DLASET, or by DO 00240 * loops. Sven Hammarling, 1/5/02. 00241 * 00242 * .. Parameters .. 00243 DOUBLE PRECISION ZERO, ONE 00244 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00245 DOUBLE PRECISION TWENTY 00246 PARAMETER ( TWENTY = 2.0D+01 ) 00247 INTEGER LDST 00248 PARAMETER ( LDST = 4 ) 00249 LOGICAL WANDS 00250 PARAMETER ( WANDS = .TRUE. ) 00251 * .. 00252 * .. Local Scalars .. 00253 LOGICAL DTRONG, WEAK 00254 INTEGER I, IDUM, LINFO, M 00255 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00268 EXTERNAL DLAMCH 00269 * .. 00270 * .. External Subroutines .. 00271 EXTERNAL DGEMM, DGEQR2, DGERQ2, DLACPY, DLAGV2, DLARTG, 00272 $ DLASET, DLASSQ, DORG2R, DORGR2, DORM2R, DORMR2, 00273 $ DROT, DSCAL, DTGSY2 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( 1, N*M, M*M*2 ) ) THEN 00290 INFO = -16 00291 WORK( 1 ) = MAX( 1, N*M, M*M*2 ) 00292 RETURN 00293 END IF 00294 * 00295 WEAK = .FALSE. 00296 DTRONG = .FALSE. 00297 * 00298 * Make a local copy of selected block 00299 * 00300 CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, LI, LDST ) 00301 CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, IR, LDST ) 00302 CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST ) 00303 CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST ) 00304 * 00305 * Compute threshold for testing acceptance of swapping. 00306 * 00307 EPS = DLAMCH( 'P' ) 00308 SMLNUM = DLAMCH( 'S' ) / EPS 00309 DSCALE = ZERO 00310 DSUM = ONE 00311 CALL DLACPY( 'Full', M, M, S, LDST, WORK, M ) 00312 CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) 00313 CALL DLACPY( 'Full', M, M, T, LDST, WORK, M ) 00314 CALL DLASSQ( 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 DLARTG( 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 DROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ), 00342 $ IR( 2, 1 ) ) 00343 CALL DROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ), 00344 $ IR( 2, 1 ) ) 00345 IF( SA.GE.SB ) THEN 00346 CALL DLARTG( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ), 00347 $ DDUM ) 00348 ELSE 00349 CALL DLARTG( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ), 00350 $ DDUM ) 00351 END IF 00352 CALL DROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ), 00353 $ LI( 2, 1 ) ) 00354 CALL DROT( 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 DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), 00373 $ M ) 00374 CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO, 00375 $ WORK, M ) 00376 CALL DGEMM( '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 DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) 00381 * 00382 CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), 00383 $ M ) 00384 CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO, 00385 $ WORK, M ) 00386 CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE, 00387 $ WORK( M*M+1 ), M ) 00388 CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) 00389 SS = DSCALE*SQRT( DSUM ) 00390 DTRONG = SS.LE.THRESH 00391 IF( .NOT.DTRONG ) 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 DROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ), 00399 $ IR( 2, 1 ) ) 00400 CALL DROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ), 00401 $ IR( 2, 1 ) ) 00402 CALL DROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, 00403 $ LI( 1, 1 ), LI( 2, 1 ) ) 00404 CALL DROT( 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 DROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ), 00416 $ IR( 2, 1 ) ) 00417 IF( WANTQ ) 00418 $ CALL DROT( 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 DLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST ) 00436 CALL DLACPY( 'Full', N1, N2, S( 1, N1+1 ), LDST, 00437 $ IR( N2+1, N1+1 ), LDST ) 00438 CALL DTGSY2( '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 DSCAL( N1, -ONE, LI( 1, I ), 1 ) 00453 LI( N1+I, I ) = SCALE 00454 10 CONTINUE 00455 CALL DGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO ) 00456 IF( LINFO.NE.0 ) 00457 $ GO TO 70 00458 CALL DORG2R( 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 DGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO ) 00472 IF( LINFO.NE.0 ) 00473 $ GO TO 70 00474 CALL DORGR2( 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 DGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO, 00481 $ WORK, M ) 00482 CALL DGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, S, 00483 $ LDST ) 00484 CALL DGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO, 00485 $ WORK, M ) 00486 CALL DGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO, T, 00487 $ LDST ) 00488 CALL DLACPY( 'F', M, M, S, LDST, SCPY, LDST ) 00489 CALL DLACPY( 'F', M, M, T, LDST, TCPY, LDST ) 00490 CALL DLACPY( 'F', M, M, IR, LDST, IRCOP, LDST ) 00491 CALL DLACPY( '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 DGERQ2( M, M, T, LDST, TAUR, WORK, LINFO ) 00497 IF( LINFO.NE.0 ) 00498 $ GO TO 70 00499 CALL DORMR2( 'R', 'T', M, M, M, T, LDST, TAUR, S, LDST, WORK, 00500 $ LINFO ) 00501 IF( LINFO.NE.0 ) 00502 $ GO TO 70 00503 CALL DORMR2( '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 DLASSQ( 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 DGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO ) 00521 IF( LINFO.NE.0 ) 00522 $ GO TO 70 00523 CALL DORM2R( 'L', 'T', M, M, M, TCPY, LDST, TAUL, SCPY, LDST, 00524 $ WORK, INFO ) 00525 CALL DORM2R( '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 DLASSQ( 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 DLACPY( 'F', M, M, SCPY, LDST, S, LDST ) 00545 CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST ) 00546 CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST ) 00547 CALL DLACPY( '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 DLASET( '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 DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), 00562 $ M ) 00563 CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO, 00564 $ WORK, M ) 00565 CALL DGEMM( '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 DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) 00570 * 00571 CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), 00572 $ M ) 00573 CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO, 00574 $ WORK, M ) 00575 CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE, 00576 $ WORK( M*M+1 ), M ) 00577 CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) 00578 SS = DSCALE*SQRT( DSUM ) 00579 DTRONG = ( SS.LE.THRESH ) 00580 IF( .NOT.DTRONG ) 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 DLASET( '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 DLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA ) 00593 CALL DLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB ) 00594 CALL DLASET( '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 DLAGV2( 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 DLAGV2( 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 DGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, A( J1, J1+N2 ), 00626 $ LDA, ZERO, WORK( M*M+1 ), N2 ) 00627 CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ), 00628 $ LDA ) 00629 CALL DGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, B( J1, J1+N2 ), 00630 $ LDB, ZERO, WORK( M*M+1 ), N2 ) 00631 CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ), 00632 $ LDB ) 00633 CALL DGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, WORK, M, ZERO, 00634 $ WORK( M*M+1 ), M ) 00635 CALL DLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST ) 00636 CALL DGEMM( 'N', 'N', N2, N1, N1, ONE, A( J1, J1+N2 ), LDA, 00637 $ T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 ) 00638 CALL DLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA ) 00639 CALL DGEMM( 'N', 'N', N2, N1, N1, ONE, B( J1, J1+N2 ), LDB, 00640 $ T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 ) 00641 CALL DLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB ) 00642 CALL DGEMM( 'T', 'N', M, M, M, ONE, IR, LDST, T, LDST, ZERO, 00643 $ WORK, M ) 00644 CALL DLACPY( 'Full', M, M, WORK, M, IR, LDST ) 00645 * 00646 * Accumulate transformations into Q and Z if requested. 00647 * 00648 IF( WANTQ ) THEN 00649 CALL DGEMM( 'N', 'N', N, M, M, ONE, Q( 1, J1 ), LDQ, LI, 00650 $ LDST, ZERO, WORK, N ) 00651 CALL DLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ ) 00652 * 00653 END IF 00654 * 00655 IF( WANTZ ) THEN 00656 CALL DGEMM( 'N', 'N', N, M, M, ONE, Z( 1, J1 ), LDZ, IR, 00657 $ LDST, ZERO, WORK, N ) 00658 CALL DLACPY( '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 DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST, 00668 $ A( J1, I ), LDA, ZERO, WORK, M ) 00669 CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA ) 00670 CALL DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST, 00671 $ B( J1, I ), LDA, ZERO, WORK, M ) 00672 CALL DLACPY( '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 DGEMM( 'N', 'N', I, M, M, ONE, A( 1, J1 ), LDA, IR, 00677 $ LDST, ZERO, WORK, I ) 00678 CALL DLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA ) 00679 CALL DGEMM( 'N', 'N', I, M, M, ONE, B( 1, J1 ), LDB, IR, 00680 $ LDST, ZERO, WORK, I ) 00681 CALL DLACPY( '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 DTGEX2 00698 * 00699 END