1 SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 2 $ INFO ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* February 29, 1992 8* 9* .. Scalar Arguments .. 10 LOGICAL WANTQ 11 INTEGER INFO, J1, LDQ, LDT, N, N1, N2 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 15* .. 16* 17* Purpose 18* ======= 19* 20* DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in 21* an upper quasi-triangular matrix T by an orthogonal similarity 22* transformation. 23* 24* T must be in Schur canonical form, that is, block upper triangular 25* with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block 26* has its diagonal elemnts equal and its off-diagonal elements of 27* opposite sign. 28* 29* Arguments 30* ========= 31* 32* WANTQ (input) LOGICAL 33* = .TRUE. : accumulate the transformation in the matrix Q; 34* = .FALSE.: do not accumulate the transformation. 35* 36* N (input) INTEGER 37* The order of the matrix T. N >= 0. 38* 39* T (input/output) DOUBLE PRECISION array, dimension (LDT,N) 40* On entry, the upper quasi-triangular matrix T, in Schur 41* canonical form. 42* On exit, the updated matrix T, again in Schur canonical form. 43* 44* LDT (input) INTEGER 45* The leading dimension of the array T. LDT >= max(1,N). 46* 47* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 48* On entry, if WANTQ is .TRUE., the orthogonal matrix Q. 49* On exit, if WANTQ is .TRUE., the updated matrix Q. 50* If WANTQ is .FALSE., Q is not referenced. 51* 52* LDQ (input) INTEGER 53* The leading dimension of the array Q. 54* LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. 55* 56* J1 (input) INTEGER 57* The index of the first row of the first block T11. 58* 59* N1 (input) INTEGER 60* The order of the first block T11. N1 = 0, 1 or 2. 61* 62* N2 (input) INTEGER 63* The order of the second block T22. N2 = 0, 1 or 2. 64* 65* WORK (workspace) DOUBLE PRECISION array, dimension (N) 66* 67* INFO (output) INTEGER 68* = 0: successful exit 69* = 1: the transformed matrix T would be too far from Schur 70* form; the blocks are not swapped and T and Q are 71* unchanged. 72* 73* ===================================================================== 74* 75* .. Parameters .. 76 DOUBLE PRECISION ZERO, ONE 77 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 78 DOUBLE PRECISION TEN 79 PARAMETER ( TEN = 1.0D+1 ) 80 INTEGER LDD, LDX 81 PARAMETER ( LDD = 4, LDX = 2 ) 82* .. 83* .. Local Scalars .. 84 INTEGER IERR, J2, J3, J4, K, ND 85 DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22, 86 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2, 87 $ WR1, WR2, XNORM 88* .. 89* .. Local Arrays .. 90 DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ), 91 $ X( LDX, 2 ) 92* .. 93* .. External Functions .. 94 DOUBLE PRECISION DLAMCH, DLANGE 95 EXTERNAL DLAMCH, DLANGE 96* .. 97* .. External Subroutines .. 98 EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2, 99 $ DROT 100* .. 101* .. Intrinsic Functions .. 102 INTRINSIC ABS, MAX 103* .. 104* .. Executable Statements .. 105* 106 INFO = 0 107* 108* Quick return if possible 109* 110 IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 ) 111 $ RETURN 112 IF( J1+N1.GT.N ) 113 $ RETURN 114* 115 J2 = J1 + 1 116 J3 = J1 + 2 117 J4 = J1 + 3 118* 119 IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN 120* 121* Swap two 1-by-1 blocks. 122* 123 T11 = T( J1, J1 ) 124 T22 = T( J2, J2 ) 125* 126* Determine the transformation to perform the interchange. 127* 128 CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP ) 129* 130* Apply transformation to the matrix T. 131* 132 IF( J3.LE.N ) 133 $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS, 134 $ SN ) 135 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 136* 137 T( J1, J1 ) = T22 138 T( J2, J2 ) = T11 139* 140 IF( WANTQ ) THEN 141* 142* Accumulate transformation in the matrix Q. 143* 144 CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 145 END IF 146* 147 ELSE 148* 149* Swapping involves at least one 2-by-2 block. 150* 151* Copy the diagonal block of order N1+N2 to the local array D 152* and compute its norm. 153* 154 ND = N1 + N2 155 CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD ) 156 DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK ) 157* 158* Compute machine-dependent threshold for test for accepting 159* swap. 160* 161 EPS = DLAMCH( 'P' ) 162 SMLNUM = DLAMCH( 'S' ) / EPS 163 THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) 164* 165* Solve T11*X - X*T22 = scale*T12 for X. 166* 167 CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD, 168 $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X, 169 $ LDX, XNORM, IERR ) 170* 171* Swap the adjacent diagonal blocks. 172* 173 K = N1 + N1 + N2 - 3 174 GO TO ( 10, 20, 30 )K 175* 176 10 CONTINUE 177* 178* N1 = 1, N2 = 2: generate elementary reflector H so that: 179* 180* ( scale, X11, X12 ) H = ( 0, 0, * ) 181* 182 U( 1 ) = SCALE 183 U( 2 ) = X( 1, 1 ) 184 U( 3 ) = X( 1, 2 ) 185 CALL DLARFG( 3, U( 3 ), U, 1, TAU ) 186 U( 3 ) = ONE 187 T11 = T( J1, J1 ) 188* 189* Perform swap provisionally on diagonal block in D. 190* 191 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 192 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 193* 194* Test whether to reject swap. 195* 196 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3, 197 $ 3 )-T11 ) ).GT.THRESH )GO TO 50 198* 199* Accept swap: apply transformation to the entire matrix T. 200* 201 CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK ) 202 CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 203* 204 T( J3, J1 ) = ZERO 205 T( J3, J2 ) = ZERO 206 T( J3, J3 ) = T11 207* 208 IF( WANTQ ) THEN 209* 210* Accumulate transformation in the matrix Q. 211* 212 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 213 END IF 214 GO TO 40 215* 216 20 CONTINUE 217* 218* N1 = 2, N2 = 1: generate elementary reflector H so that: 219* 220* H ( -X11 ) = ( * ) 221* ( -X21 ) = ( 0 ) 222* ( scale ) = ( 0 ) 223* 224 U( 1 ) = -X( 1, 1 ) 225 U( 2 ) = -X( 2, 1 ) 226 U( 3 ) = SCALE 227 CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU ) 228 U( 1 ) = ONE 229 T33 = T( J3, J3 ) 230* 231* Perform swap provisionally on diagonal block in D. 232* 233 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 234 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 235* 236* Test whether to reject swap. 237* 238 IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1, 239 $ 1 )-T33 ) ).GT.THRESH )GO TO 50 240* 241* Accept swap: apply transformation to the entire matrix T. 242* 243 CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 244 CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK ) 245* 246 T( J1, J1 ) = T33 247 T( J2, J1 ) = ZERO 248 T( J3, J1 ) = ZERO 249* 250 IF( WANTQ ) THEN 251* 252* Accumulate transformation in the matrix Q. 253* 254 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 255 END IF 256 GO TO 40 257* 258 30 CONTINUE 259* 260* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so 261* that: 262* 263* H(2) H(1) ( -X11 -X12 ) = ( * * ) 264* ( -X21 -X22 ) ( 0 * ) 265* ( scale 0 ) ( 0 0 ) 266* ( 0 scale ) ( 0 0 ) 267* 268 U1( 1 ) = -X( 1, 1 ) 269 U1( 2 ) = -X( 2, 1 ) 270 U1( 3 ) = SCALE 271 CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 ) 272 U1( 1 ) = ONE 273* 274 TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) ) 275 U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 ) 276 U2( 2 ) = -TEMP*U1( 3 ) 277 U2( 3 ) = SCALE 278 CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 ) 279 U2( 1 ) = ONE 280* 281* Perform swap provisionally on diagonal block in D. 282* 283 CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK ) 284 CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK ) 285 CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK ) 286 CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK ) 287* 288* Test whether to reject swap. 289* 290 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ), 291 $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50 292* 293* Accept swap: apply transformation to the entire matrix T. 294* 295 CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK ) 296 CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK ) 297 CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK ) 298 CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK ) 299* 300 T( J3, J1 ) = ZERO 301 T( J3, J2 ) = ZERO 302 T( J4, J1 ) = ZERO 303 T( J4, J2 ) = ZERO 304* 305 IF( WANTQ ) THEN 306* 307* Accumulate transformation in the matrix Q. 308* 309 CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK ) 310 CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK ) 311 END IF 312* 313 40 CONTINUE 314* 315 IF( N2.EQ.2 ) THEN 316* 317* Standardize new 2-by-2 block T11 318* 319 CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ), 320 $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN ) 321 CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT, 322 $ CS, SN ) 323 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 324 IF( WANTQ ) 325 $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 326 END IF 327* 328 IF( N1.EQ.2 ) THEN 329* 330* Standardize new 2-by-2 block T22 331* 332 J3 = J1 + N2 333 J4 = J3 + 1 334 CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ), 335 $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN ) 336 IF( J3+2.LE.N ) 337 $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ), 338 $ LDT, CS, SN ) 339 CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN ) 340 IF( WANTQ ) 341 $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN ) 342 END IF 343* 344 END IF 345 RETURN 346* 347* Exit with INFO = 1 if swap was rejected. 348* 349 50 CONTINUE 350 INFO = 1 351 RETURN 352* 353* End of DLAEXC 354* 355 END 356