|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL WANTQ 00011 INTEGER INFO, J1, LDQ, LDT, N, N1, N2 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in 00021 * an upper quasi-triangular matrix T by an orthogonal similarity 00022 * transformation. 00023 * 00024 * T must be in Schur canonical form, that is, block upper triangular 00025 * with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block 00026 * has its diagonal elemnts equal and its off-diagonal elements of 00027 * opposite sign. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * WANTQ (input) LOGICAL 00033 * = .TRUE. : accumulate the transformation in the matrix Q; 00034 * = .FALSE.: do not accumulate the transformation. 00035 * 00036 * N (input) INTEGER 00037 * The order of the matrix T. N >= 0. 00038 * 00039 * T (input/output) DOUBLE PRECISION array, dimension (LDT,N) 00040 * On entry, the upper quasi-triangular matrix T, in Schur 00041 * canonical form. 00042 * On exit, the updated matrix T, again in Schur canonical form. 00043 * 00044 * LDT (input) INTEGER 00045 * The leading dimension of the array T. LDT >= max(1,N). 00046 * 00047 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 00048 * On entry, if WANTQ is .TRUE., the orthogonal matrix Q. 00049 * On exit, if WANTQ is .TRUE., the updated matrix Q. 00050 * If WANTQ is .FALSE., Q is not referenced. 00051 * 00052 * LDQ (input) INTEGER 00053 * The leading dimension of the array Q. 00054 * LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. 00055 * 00056 * J1 (input) INTEGER 00057 * The index of the first row of the first block T11. 00058 * 00059 * N1 (input) INTEGER 00060 * The order of the first block T11. N1 = 0, 1 or 2. 00061 * 00062 * N2 (input) INTEGER 00063 * The order of the second block T22. N2 = 0, 1 or 2. 00064 * 00065 * WORK (workspace) DOUBLE PRECISION array, dimension (N) 00066 * 00067 * INFO (output) INTEGER 00068 * = 0: successful exit 00069 * = 1: the transformed matrix T would be too far from Schur 00070 * form; the blocks are not swapped and T and Q are 00071 * unchanged. 00072 * 00073 * ===================================================================== 00074 * 00075 * .. Parameters .. 00076 DOUBLE PRECISION ZERO, ONE 00077 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00078 DOUBLE PRECISION TEN 00079 PARAMETER ( TEN = 1.0D+1 ) 00080 INTEGER LDD, LDX 00081 PARAMETER ( LDD = 4, LDX = 2 ) 00082 * .. 00083 * .. Local Scalars .. 00084 INTEGER IERR, J2, J3, J4, K, ND 00085 DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22, 00086 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2, 00087 $ WR1, WR2, XNORM 00088 * .. 00089 * .. Local Arrays .. 00090 DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ), 00091 $ X( LDX, 2 ) 00092 * .. 00093 * .. External Functions .. 00094 DOUBLE PRECISION DLAMCH, DLANGE 00095 EXTERNAL DLAMCH, DLANGE 00096 * .. 00097 * .. External Subroutines .. 00098 EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2, 00099 $ DROT 00100 * .. 00101 * .. Intrinsic Functions .. 00102 INTRINSIC ABS, MAX 00103 * .. 00104 * .. Executable Statements .. 00105 * 00106 INFO = 0 00107 * 00108 * Quick return if possible 00109 * 00110 IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 ) 00111 $ RETURN 00112 IF( J1+N1.GT.N ) 00113 $ RETURN 00114 * 00115 J2 = J1 + 1 00116 J3 = J1 + 2 00117 J4 = J1 + 3 00118 * 00119 IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN 00120 * 00121 * Swap two 1-by-1 blocks. 00122 * 00123 T11 = T( J1, J1 ) 00124 T22 = T( J2, J2 ) 00125 * 00126 * Determine the transformation to perform the interchange. 00127 * 00128 CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP ) 00129 * 00130 * Apply transformation to the matrix T. 00131 * 00132 IF( J3.LE.N ) 00133 $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS, 00134 $ SN ) 00135 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 00136 * 00137 T( J1, J1 ) = T22 00138 T( J2, J2 ) = T11 00139 * 00140 IF( WANTQ ) THEN 00141 * 00142 * Accumulate transformation in the matrix Q. 00143 * 00144 CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 00145 END IF 00146 * 00147 ELSE 00148 * 00149 * Swapping involves at least one 2-by-2 block. 00150 * 00151 * Copy the diagonal block of order N1+N2 to the local array D 00152 * and compute its norm. 00153 * 00154 ND = N1 + N2 00155 CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD ) 00156 DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK ) 00157 * 00158 * Compute machine-dependent threshold for test for accepting 00159 * swap. 00160 * 00161 EPS = DLAMCH( 'P' ) 00162 SMLNUM = DLAMCH( 'S' ) / EPS 00163 THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) 00164 * 00165 * Solve T11*X - X*T22 = scale*T12 for X. 00166 * 00167 CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD, 00168 $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X, 00169 $ LDX, XNORM, IERR ) 00170 * 00171 * Swap the adjacent diagonal blocks. 00172 * 00173 K = N1 + N1 + N2 - 3 00174 GO TO ( 10, 20, 30 )K 00175 * 00176 10 CONTINUE 00177 * 00178 * N1 = 1, N2 = 2: generate elementary reflector H so that: 00179 * 00180 * ( scale, X11, X12 ) H = ( 0, 0, * ) 00181 * 00182 U( 1 ) = SCALE 00183 U( 2 ) = X( 1, 1 ) 00184 U( 3 ) = X( 1, 2 ) 00185 CALL DLARFG( 3, U( 3 ), U, 1, TAU ) 00186 U( 3 ) = ONE 00187 T11 = T( J1, J1 ) 00188 * 00189 * Perform swap provisionally on diagonal block in D. 00190 * 00191 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 00192 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 00193 * 00194 * Test whether to reject swap. 00195 * 00196 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3, 00197 $ 3 )-T11 ) ).GT.THRESH )GO TO 50 00198 * 00199 * Accept swap: apply transformation to the entire matrix T. 00200 * 00201 CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK ) 00202 CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 00203 * 00204 T( J3, J1 ) = ZERO 00205 T( J3, J2 ) = ZERO 00206 T( J3, J3 ) = T11 00207 * 00208 IF( WANTQ ) THEN 00209 * 00210 * Accumulate transformation in the matrix Q. 00211 * 00212 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 00213 END IF 00214 GO TO 40 00215 * 00216 20 CONTINUE 00217 * 00218 * N1 = 2, N2 = 1: generate elementary reflector H so that: 00219 * 00220 * H ( -X11 ) = ( * ) 00221 * ( -X21 ) = ( 0 ) 00222 * ( scale ) = ( 0 ) 00223 * 00224 U( 1 ) = -X( 1, 1 ) 00225 U( 2 ) = -X( 2, 1 ) 00226 U( 3 ) = SCALE 00227 CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU ) 00228 U( 1 ) = ONE 00229 T33 = T( J3, J3 ) 00230 * 00231 * Perform swap provisionally on diagonal block in D. 00232 * 00233 CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 00234 CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 00235 * 00236 * Test whether to reject swap. 00237 * 00238 IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1, 00239 $ 1 )-T33 ) ).GT.THRESH )GO TO 50 00240 * 00241 * Accept swap: apply transformation to the entire matrix T. 00242 * 00243 CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 00244 CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK ) 00245 * 00246 T( J1, J1 ) = T33 00247 T( J2, J1 ) = ZERO 00248 T( J3, J1 ) = ZERO 00249 * 00250 IF( WANTQ ) THEN 00251 * 00252 * Accumulate transformation in the matrix Q. 00253 * 00254 CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 00255 END IF 00256 GO TO 40 00257 * 00258 30 CONTINUE 00259 * 00260 * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so 00261 * that: 00262 * 00263 * H(2) H(1) ( -X11 -X12 ) = ( * * ) 00264 * ( -X21 -X22 ) ( 0 * ) 00265 * ( scale 0 ) ( 0 0 ) 00266 * ( 0 scale ) ( 0 0 ) 00267 * 00268 U1( 1 ) = -X( 1, 1 ) 00269 U1( 2 ) = -X( 2, 1 ) 00270 U1( 3 ) = SCALE 00271 CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 ) 00272 U1( 1 ) = ONE 00273 * 00274 TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) ) 00275 U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 ) 00276 U2( 2 ) = -TEMP*U1( 3 ) 00277 U2( 3 ) = SCALE 00278 CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 ) 00279 U2( 1 ) = ONE 00280 * 00281 * Perform swap provisionally on diagonal block in D. 00282 * 00283 CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK ) 00284 CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK ) 00285 CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK ) 00286 CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK ) 00287 * 00288 * Test whether to reject swap. 00289 * 00290 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ), 00291 $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50 00292 * 00293 * Accept swap: apply transformation to the entire matrix T. 00294 * 00295 CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK ) 00296 CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK ) 00297 CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK ) 00298 CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK ) 00299 * 00300 T( J3, J1 ) = ZERO 00301 T( J3, J2 ) = ZERO 00302 T( J4, J1 ) = ZERO 00303 T( J4, J2 ) = ZERO 00304 * 00305 IF( WANTQ ) THEN 00306 * 00307 * Accumulate transformation in the matrix Q. 00308 * 00309 CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK ) 00310 CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK ) 00311 END IF 00312 * 00313 40 CONTINUE 00314 * 00315 IF( N2.EQ.2 ) THEN 00316 * 00317 * Standardize new 2-by-2 block T11 00318 * 00319 CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ), 00320 $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN ) 00321 CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT, 00322 $ CS, SN ) 00323 CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 00324 IF( WANTQ ) 00325 $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 00326 END IF 00327 * 00328 IF( N1.EQ.2 ) THEN 00329 * 00330 * Standardize new 2-by-2 block T22 00331 * 00332 J3 = J1 + N2 00333 J4 = J3 + 1 00334 CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ), 00335 $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN ) 00336 IF( J3+2.LE.N ) 00337 $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ), 00338 $ LDT, CS, SN ) 00339 CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN ) 00340 IF( WANTQ ) 00341 $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN ) 00342 END IF 00343 * 00344 END IF 00345 RETURN 00346 * 00347 * Exit with INFO = 1 if swap was rejected. 00348 * 00349 50 CONTINUE 00350 INFO = 1 00351 RETURN 00352 * 00353 * End of DLAEXC 00354 * 00355 END
1.7.3 
