![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAEXC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAEXC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaexc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaexc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaexc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL WANTQ 00026 * INTEGER INFO, J1, LDQ, LDT, N, N1, N2 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL Q( LDQ, * ), T( LDT, * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in 00039 *> an upper quasi-triangular matrix T by an orthogonal similarity 00040 *> transformation. 00041 *> 00042 *> T must be in Schur canonical form, that is, block upper triangular 00043 *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block 00044 *> has its diagonal elemnts equal and its off-diagonal elements of 00045 *> opposite sign. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] WANTQ 00052 *> \verbatim 00053 *> WANTQ is LOGICAL 00054 *> = .TRUE. : accumulate the transformation in the matrix Q; 00055 *> = .FALSE.: do not accumulate the transformation. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The order of the matrix T. N >= 0. 00062 *> \endverbatim 00063 *> 00064 *> \param[in,out] T 00065 *> \verbatim 00066 *> T is REAL array, dimension (LDT,N) 00067 *> On entry, the upper quasi-triangular matrix T, in Schur 00068 *> canonical form. 00069 *> On exit, the updated matrix T, again in Schur canonical form. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] LDT 00073 *> \verbatim 00074 *> LDT is INTEGER 00075 *> The leading dimension of the array T. LDT >= max(1,N). 00076 *> \endverbatim 00077 *> 00078 *> \param[in,out] Q 00079 *> \verbatim 00080 *> Q is REAL array, dimension (LDQ,N) 00081 *> On entry, if WANTQ is .TRUE., the orthogonal matrix Q. 00082 *> On exit, if WANTQ is .TRUE., the updated matrix Q. 00083 *> If WANTQ is .FALSE., Q is not referenced. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDQ 00087 *> \verbatim 00088 *> LDQ is INTEGER 00089 *> The leading dimension of the array Q. 00090 *> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. 00091 *> \endverbatim 00092 *> 00093 *> \param[in] J1 00094 *> \verbatim 00095 *> J1 is INTEGER 00096 *> The index of the first row of the first block T11. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] N1 00100 *> \verbatim 00101 *> N1 is INTEGER 00102 *> The order of the first block T11. N1 = 0, 1 or 2. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] N2 00106 *> \verbatim 00107 *> N2 is INTEGER 00108 *> The order of the second block T22. N2 = 0, 1 or 2. 00109 *> \endverbatim 00110 *> 00111 *> \param[out] WORK 00112 *> \verbatim 00113 *> WORK is REAL array, dimension (N) 00114 *> \endverbatim 00115 *> 00116 *> \param[out] INFO 00117 *> \verbatim 00118 *> INFO is INTEGER 00119 *> = 0: successful exit 00120 *> = 1: the transformed matrix T would be too far from Schur 00121 *> form; the blocks are not swapped and T and Q are 00122 *> unchanged. 00123 *> \endverbatim 00124 * 00125 * Authors: 00126 * ======== 00127 * 00128 *> \author Univ. of Tennessee 00129 *> \author Univ. of California Berkeley 00130 *> \author Univ. of Colorado Denver 00131 *> \author NAG Ltd. 00132 * 00133 *> \date November 2011 00134 * 00135 *> \ingroup realOTHERauxiliary 00136 * 00137 * ===================================================================== 00138 SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, 00139 $ INFO ) 00140 * 00141 * -- LAPACK auxiliary routine (version 3.4.0) -- 00142 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00144 * November 2011 00145 * 00146 * .. Scalar Arguments .. 00147 LOGICAL WANTQ 00148 INTEGER INFO, J1, LDQ, LDT, N, N1, N2 00149 * .. 00150 * .. Array Arguments .. 00151 REAL Q( LDQ, * ), T( LDT, * ), WORK( * ) 00152 * .. 00153 * 00154 * ===================================================================== 00155 * 00156 * .. Parameters .. 00157 REAL ZERO, ONE 00158 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00159 REAL TEN 00160 PARAMETER ( TEN = 1.0E+1 ) 00161 INTEGER LDD, LDX 00162 PARAMETER ( LDD = 4, LDX = 2 ) 00163 * .. 00164 * .. Local Scalars .. 00165 INTEGER IERR, J2, J3, J4, K, ND 00166 REAL CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22, 00167 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2, 00168 $ WR1, WR2, XNORM 00169 * .. 00170 * .. Local Arrays .. 00171 REAL D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ), 00172 $ X( LDX, 2 ) 00173 * .. 00174 * .. External Functions .. 00175 REAL SLAMCH, SLANGE 00176 EXTERNAL SLAMCH, SLANGE 00177 * .. 00178 * .. External Subroutines .. 00179 EXTERNAL SLACPY, SLANV2, SLARFG, SLARFX, SLARTG, SLASY2, 00180 $ SROT 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC ABS, MAX 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 INFO = 0 00188 * 00189 * Quick return if possible 00190 * 00191 IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 ) 00192 $ RETURN 00193 IF( J1+N1.GT.N ) 00194 $ RETURN 00195 * 00196 J2 = J1 + 1 00197 J3 = J1 + 2 00198 J4 = J1 + 3 00199 * 00200 IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN 00201 * 00202 * Swap two 1-by-1 blocks. 00203 * 00204 T11 = T( J1, J1 ) 00205 T22 = T( J2, J2 ) 00206 * 00207 * Determine the transformation to perform the interchange. 00208 * 00209 CALL SLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP ) 00210 * 00211 * Apply transformation to the matrix T. 00212 * 00213 IF( J3.LE.N ) 00214 $ CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS, 00215 $ SN ) 00216 CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 00217 * 00218 T( J1, J1 ) = T22 00219 T( J2, J2 ) = T11 00220 * 00221 IF( WANTQ ) THEN 00222 * 00223 * Accumulate transformation in the matrix Q. 00224 * 00225 CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 00226 END IF 00227 * 00228 ELSE 00229 * 00230 * Swapping involves at least one 2-by-2 block. 00231 * 00232 * Copy the diagonal block of order N1+N2 to the local array D 00233 * and compute its norm. 00234 * 00235 ND = N1 + N2 00236 CALL SLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD ) 00237 DNORM = SLANGE( 'Max', ND, ND, D, LDD, WORK ) 00238 * 00239 * Compute machine-dependent threshold for test for accepting 00240 * swap. 00241 * 00242 EPS = SLAMCH( 'P' ) 00243 SMLNUM = SLAMCH( 'S' ) / EPS 00244 THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) 00245 * 00246 * Solve T11*X - X*T22 = scale*T12 for X. 00247 * 00248 CALL SLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD, 00249 $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X, 00250 $ LDX, XNORM, IERR ) 00251 * 00252 * Swap the adjacent diagonal blocks. 00253 * 00254 K = N1 + N1 + N2 - 3 00255 GO TO ( 10, 20, 30 )K 00256 * 00257 10 CONTINUE 00258 * 00259 * N1 = 1, N2 = 2: generate elementary reflector H so that: 00260 * 00261 * ( scale, X11, X12 ) H = ( 0, 0, * ) 00262 * 00263 U( 1 ) = SCALE 00264 U( 2 ) = X( 1, 1 ) 00265 U( 3 ) = X( 1, 2 ) 00266 CALL SLARFG( 3, U( 3 ), U, 1, TAU ) 00267 U( 3 ) = ONE 00268 T11 = T( J1, J1 ) 00269 * 00270 * Perform swap provisionally on diagonal block in D. 00271 * 00272 CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 00273 CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 00274 * 00275 * Test whether to reject swap. 00276 * 00277 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3, 00278 $ 3 )-T11 ) ).GT.THRESH )GO TO 50 00279 * 00280 * Accept swap: apply transformation to the entire matrix T. 00281 * 00282 CALL SLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK ) 00283 CALL SLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 00284 * 00285 T( J3, J1 ) = ZERO 00286 T( J3, J2 ) = ZERO 00287 T( J3, J3 ) = T11 00288 * 00289 IF( WANTQ ) THEN 00290 * 00291 * Accumulate transformation in the matrix Q. 00292 * 00293 CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 00294 END IF 00295 GO TO 40 00296 * 00297 20 CONTINUE 00298 * 00299 * N1 = 2, N2 = 1: generate elementary reflector H so that: 00300 * 00301 * H ( -X11 ) = ( * ) 00302 * ( -X21 ) = ( 0 ) 00303 * ( scale ) = ( 0 ) 00304 * 00305 U( 1 ) = -X( 1, 1 ) 00306 U( 2 ) = -X( 2, 1 ) 00307 U( 3 ) = SCALE 00308 CALL SLARFG( 3, U( 1 ), U( 2 ), 1, TAU ) 00309 U( 1 ) = ONE 00310 T33 = T( J3, J3 ) 00311 * 00312 * Perform swap provisionally on diagonal block in D. 00313 * 00314 CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) 00315 CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) 00316 * 00317 * Test whether to reject swap. 00318 * 00319 IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1, 00320 $ 1 )-T33 ) ).GT.THRESH )GO TO 50 00321 * 00322 * Accept swap: apply transformation to the entire matrix T. 00323 * 00324 CALL SLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK ) 00325 CALL SLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK ) 00326 * 00327 T( J1, J1 ) = T33 00328 T( J2, J1 ) = ZERO 00329 T( J3, J1 ) = ZERO 00330 * 00331 IF( WANTQ ) THEN 00332 * 00333 * Accumulate transformation in the matrix Q. 00334 * 00335 CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) 00336 END IF 00337 GO TO 40 00338 * 00339 30 CONTINUE 00340 * 00341 * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so 00342 * that: 00343 * 00344 * H(2) H(1) ( -X11 -X12 ) = ( * * ) 00345 * ( -X21 -X22 ) ( 0 * ) 00346 * ( scale 0 ) ( 0 0 ) 00347 * ( 0 scale ) ( 0 0 ) 00348 * 00349 U1( 1 ) = -X( 1, 1 ) 00350 U1( 2 ) = -X( 2, 1 ) 00351 U1( 3 ) = SCALE 00352 CALL SLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 ) 00353 U1( 1 ) = ONE 00354 * 00355 TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) ) 00356 U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 ) 00357 U2( 2 ) = -TEMP*U1( 3 ) 00358 U2( 3 ) = SCALE 00359 CALL SLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 ) 00360 U2( 1 ) = ONE 00361 * 00362 * Perform swap provisionally on diagonal block in D. 00363 * 00364 CALL SLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK ) 00365 CALL SLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK ) 00366 CALL SLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK ) 00367 CALL SLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK ) 00368 * 00369 * Test whether to reject swap. 00370 * 00371 IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ), 00372 $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50 00373 * 00374 * Accept swap: apply transformation to the entire matrix T. 00375 * 00376 CALL SLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK ) 00377 CALL SLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK ) 00378 CALL SLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK ) 00379 CALL SLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK ) 00380 * 00381 T( J3, J1 ) = ZERO 00382 T( J3, J2 ) = ZERO 00383 T( J4, J1 ) = ZERO 00384 T( J4, J2 ) = ZERO 00385 * 00386 IF( WANTQ ) THEN 00387 * 00388 * Accumulate transformation in the matrix Q. 00389 * 00390 CALL SLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK ) 00391 CALL SLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK ) 00392 END IF 00393 * 00394 40 CONTINUE 00395 * 00396 IF( N2.EQ.2 ) THEN 00397 * 00398 * Standardize new 2-by-2 block T11 00399 * 00400 CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ), 00401 $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN ) 00402 CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT, 00403 $ CS, SN ) 00404 CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) 00405 IF( WANTQ ) 00406 $ CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) 00407 END IF 00408 * 00409 IF( N1.EQ.2 ) THEN 00410 * 00411 * Standardize new 2-by-2 block T22 00412 * 00413 J3 = J1 + N2 00414 J4 = J3 + 1 00415 CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ), 00416 $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN ) 00417 IF( J3+2.LE.N ) 00418 $ CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ), 00419 $ LDT, CS, SN ) 00420 CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN ) 00421 IF( WANTQ ) 00422 $ CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN ) 00423 END IF 00424 * 00425 END IF 00426 RETURN 00427 * 00428 * Exit with INFO = 1 if swap was rejected. 00429 * 00430 50 INFO = 1 00431 RETURN 00432 * 00433 * End of SLAEXC 00434 * 00435 END