1*> \brief \b DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAEXC + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 22* INFO ) 23* 24* .. Scalar Arguments .. 25* LOGICAL WANTQ 26* INTEGER INFO, J1, LDQ, LDT, N, N1, N2 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in 39*> an upper quasi-triangular matrix T by an orthogonal similarity 40*> transformation. 41*> 42*> T must be in Schur canonical form, that is, block upper triangular 43*> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block 44*> has its diagonal elemnts equal and its off-diagonal elements of 45*> opposite sign. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] WANTQ 52*> \verbatim 53*> WANTQ is LOGICAL 54*> = .TRUE. : accumulate the transformation in the matrix Q; 55*> = .FALSE.: do not accumulate the transformation. 56*> \endverbatim 57*> 58*> \param[in] N 59*> \verbatim 60*> N is INTEGER 61*> The order of the matrix T. N >= 0. 62*> \endverbatim 63*> 64*> \param[in,out] T 65*> \verbatim 66*> T is DOUBLE PRECISION array, dimension (LDT,N) 67*> On entry, the upper quasi-triangular matrix T, in Schur 68*> canonical form. 69*> On exit, the updated matrix T, again in Schur canonical form. 70*> \endverbatim 71*> 72*> \param[in] LDT 73*> \verbatim 74*> LDT is INTEGER 75*> The leading dimension of the array T. LDT >= max(1,N). 76*> \endverbatim 77*> 78*> \param[in,out] Q 79*> \verbatim 80*> Q is DOUBLE PRECISION array, dimension (LDQ,N) 81*> On entry, if WANTQ is .TRUE., the orthogonal matrix Q. 82*> On exit, if WANTQ is .TRUE., the updated matrix Q. 83*> If WANTQ is .FALSE., Q is not referenced. 84*> \endverbatim 85*> 86*> \param[in] LDQ 87*> \verbatim 88*> LDQ is INTEGER 89*> The leading dimension of the array Q. 90*> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. 91*> \endverbatim 92*> 93*> \param[in] J1 94*> \verbatim 95*> J1 is INTEGER 96*> The index of the first row of the first block T11. 97*> \endverbatim 98*> 99*> \param[in] N1 100*> \verbatim 101*> N1 is INTEGER 102*> The order of the first block T11. N1 = 0, 1 or 2. 103*> \endverbatim 104*> 105*> \param[in] N2 106*> \verbatim 107*> N2 is INTEGER 108*> The order of the second block T22. N2 = 0, 1 or 2. 109*> \endverbatim 110*> 111*> \param[out] WORK 112*> \verbatim 113*> WORK is DOUBLE PRECISION array, dimension (N) 114*> \endverbatim 115*> 116*> \param[out] INFO 117*> \verbatim 118*> INFO is INTEGER 119*> = 0: successful exit 120*> = 1: the transformed matrix T would be too far from Schur 121*> form; the blocks are not swapped and T and Q are 122*> unchanged. 123*> \endverbatim 124* 125* Authors: 126* ======== 127* 128*> \author Univ. of Tennessee 129*> \author Univ. of California Berkeley 130*> \author Univ. of Colorado Denver 131*> \author NAG Ltd. 132* 133*> \date September 2012 134* 135*> \ingroup doubleOTHERauxiliary 136* 137* ===================================================================== 138 SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 139 $ INFO ) 140* 141* -- LAPACK auxiliary routine (version 3.4.2) -- 142* -- LAPACK is a software package provided by Univ. of Tennessee, -- 143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 144* September 2012 145* 146* .. Scalar Arguments .. 147 LOGICAL WANTQ 148 INTEGER INFO, J1, LDQ, LDT, N, N1, N2 149* .. 150* .. Array Arguments .. 151 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 152* .. 153* 154* ===================================================================== 155* 156* .. Parameters .. 157 DOUBLE PRECISION ZERO, ONE 158 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 159 DOUBLE PRECISION TEN 160 PARAMETER ( TEN = 1.0D+1 ) 161 INTEGER LDD, LDX 162 PARAMETER ( LDD = 4, LDX = 2 ) 163* .. 164* .. Local Scalars .. 165 INTEGER IERR, J2, J3, J4, K, ND 166 DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22, 167 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2, 168 $ WR1, WR2, XNORM 169* .. 170* .. Local Arrays .. 171 DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ), 172 $ X( LDX, 2 ) 173* .. 174* .. External Functions .. 175 DOUBLE PRECISION DLAMCH, DLANGE 176 EXTERNAL DLAMCH, DLANGE 177* .. 178* .. External Subroutines .. 179 EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2, 180 $ DROT 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC ABS, MAX 184* .. 185* .. Executable Statements .. 186* 187 INFO = 0 188* 189* Quick return if possible 190* 191 IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 ) 192 $ RETURN 193 IF( J1+N1.GT.N ) 194 $ RETURN 195* 196 J2 = J1 + 1 197 J3 = J1 + 2 198 J4 = J1 + 3 199* 200 IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN 201* 202* Swap two 1-by-1 blocks. 203* 204 T11 = T( J1, J1 ) 205 T22 = T( J2, J2 ) 206* 207* Determine the transformation to perform the interchange. 208* 209 CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP ) 210* 211* Apply transformation to the matrix T. 212* 213 IF( J3.LE.N ) 214 $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS, 215 $ SN ) 216 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 217* 218 T( J1, J1 ) = T22 219 T( J2, J2 ) = T11 220* 221 IF( WANTQ ) THEN 222* 223* Accumulate transformation in the matrix Q. 224* 225 CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 226 END IF 227* 228 ELSE 229* 230* Swapping involves at least one 2-by-2 block. 231* 232* Copy the diagonal block of order N1+N2 to the local array D 233* and compute its norm. 234* 235 ND = N1 + N2 236 CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD ) 237 DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK ) 238* 239* Compute machine-dependent threshold for test for accepting 240* swap. 241* 242 EPS = DLAMCH( 'P' ) 243 SMLNUM = DLAMCH( 'S' ) / EPS 244 THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) 245* 246* Solve T11*X - X*T22 = scale*T12 for X. 247* 248 CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD, 249 $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X, 250 $ LDX, XNORM, IERR ) 251* 252* Swap the adjacent diagonal blocks. 253* 254 K = N1 + N1 + N2 - 3 255 GO TO ( 10, 20, 30 )K 256* 257 10 CONTINUE 258* 259* N1 = 1, N2 = 2: generate elementary reflector H so that: 260* 261* ( scale, X11, X12 ) H = ( 0, 0, * ) 262* 263 U( 1 ) = SCALE 264 U( 2 ) = X( 1, 1 ) 265 U( 3 ) = X( 1, 2 ) 266 CALL DLARFG( 3, U( 3 ), U, 1, TAU ) 267 U( 3 ) = ONE 268 T11 = T( J1, J1 ) 269* 270* Perform swap provisionally on diagonal block in D. 271* 272 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 273 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 274* 275* Test whether to reject swap. 276* 277 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3, 278 $ 3 )-T11 ) ).GT.THRESH )GO TO 50 279* 280* Accept swap: apply transformation to the entire matrix T. 281* 282 CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK ) 283 CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 284* 285 T( J3, J1 ) = ZERO 286 T( J3, J2 ) = ZERO 287 T( J3, J3 ) = T11 288* 289 IF( WANTQ ) THEN 290* 291* Accumulate transformation in the matrix Q. 292* 293 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 294 END IF 295 GO TO 40 296* 297 20 CONTINUE 298* 299* N1 = 2, N2 = 1: generate elementary reflector H so that: 300* 301* H ( -X11 ) = ( * ) 302* ( -X21 ) = ( 0 ) 303* ( scale ) = ( 0 ) 304* 305 U( 1 ) = -X( 1, 1 ) 306 U( 2 ) = -X( 2, 1 ) 307 U( 3 ) = SCALE 308 CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU ) 309 U( 1 ) = ONE 310 T33 = T( J3, J3 ) 311* 312* Perform swap provisionally on diagonal block in D. 313* 314 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 315 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 316* 317* Test whether to reject swap. 318* 319 IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1, 320 $ 1 )-T33 ) ).GT.THRESH )GO TO 50 321* 322* Accept swap: apply transformation to the entire matrix T. 323* 324 CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 325 CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK ) 326* 327 T( J1, J1 ) = T33 328 T( J2, J1 ) = ZERO 329 T( J3, J1 ) = ZERO 330* 331 IF( WANTQ ) THEN 332* 333* Accumulate transformation in the matrix Q. 334* 335 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 336 END IF 337 GO TO 40 338* 339 30 CONTINUE 340* 341* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so 342* that: 343* 344* H(2) H(1) ( -X11 -X12 ) = ( * * ) 345* ( -X21 -X22 ) ( 0 * ) 346* ( scale 0 ) ( 0 0 ) 347* ( 0 scale ) ( 0 0 ) 348* 349 U1( 1 ) = -X( 1, 1 ) 350 U1( 2 ) = -X( 2, 1 ) 351 U1( 3 ) = SCALE 352 CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 ) 353 U1( 1 ) = ONE 354* 355 TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) ) 356 U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 ) 357 U2( 2 ) = -TEMP*U1( 3 ) 358 U2( 3 ) = SCALE 359 CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 ) 360 U2( 1 ) = ONE 361* 362* Perform swap provisionally on diagonal block in D. 363* 364 CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK ) 365 CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK ) 366 CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK ) 367 CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK ) 368* 369* Test whether to reject swap. 370* 371 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ), 372 $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50 373* 374* Accept swap: apply transformation to the entire matrix T. 375* 376 CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK ) 377 CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK ) 378 CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK ) 379 CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK ) 380* 381 T( J3, J1 ) = ZERO 382 T( J3, J2 ) = ZERO 383 T( J4, J1 ) = ZERO 384 T( J4, J2 ) = ZERO 385* 386 IF( WANTQ ) THEN 387* 388* Accumulate transformation in the matrix Q. 389* 390 CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK ) 391 CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK ) 392 END IF 393* 394 40 CONTINUE 395* 396 IF( N2.EQ.2 ) THEN 397* 398* Standardize new 2-by-2 block T11 399* 400 CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ), 401 $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN ) 402 CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT, 403 $ CS, SN ) 404 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 405 IF( WANTQ ) 406 $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 407 END IF 408* 409 IF( N1.EQ.2 ) THEN 410* 411* Standardize new 2-by-2 block T22 412* 413 J3 = J1 + N2 414 J4 = J3 + 1 415 CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ), 416 $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN ) 417 IF( J3+2.LE.N ) 418 $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ), 419 $ LDT, CS, SN ) 420 CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN ) 421 IF( WANTQ ) 422 $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN ) 423 END IF 424* 425 END IF 426 RETURN 427* 428* Exit with INFO = 1 if swap was rejected. 429* 430 50 CONTINUE 431 INFO = 1 432 RETURN 433* 434* End of DLAEXC 435* 436 END 437c $Id$ 438