![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASY2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASY2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasy2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasy2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasy2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, 00022 * LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL LTRANL, LTRANR 00026 * INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2 00027 * REAL SCALE, XNORM 00028 * .. 00029 * .. Array Arguments .. 00030 * REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ), 00031 * $ X( LDX, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in 00041 *> 00042 *> op(TL)*X + ISGN*X*op(TR) = SCALE*B, 00043 *> 00044 *> where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or 00045 *> -1. op(T) = T or T**T, where T**T denotes the transpose of T. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] LTRANL 00052 *> \verbatim 00053 *> LTRANL is LOGICAL 00054 *> On entry, LTRANL specifies the op(TL): 00055 *> = .FALSE., op(TL) = TL, 00056 *> = .TRUE., op(TL) = TL**T. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] LTRANR 00060 *> \verbatim 00061 *> LTRANR is LOGICAL 00062 *> On entry, LTRANR specifies the op(TR): 00063 *> = .FALSE., op(TR) = TR, 00064 *> = .TRUE., op(TR) = TR**T. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] ISGN 00068 *> \verbatim 00069 *> ISGN is INTEGER 00070 *> On entry, ISGN specifies the sign of the equation 00071 *> as described before. ISGN may only be 1 or -1. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] N1 00075 *> \verbatim 00076 *> N1 is INTEGER 00077 *> On entry, N1 specifies the order of matrix TL. 00078 *> N1 may only be 0, 1 or 2. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] N2 00082 *> \verbatim 00083 *> N2 is INTEGER 00084 *> On entry, N2 specifies the order of matrix TR. 00085 *> N2 may only be 0, 1 or 2. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] TL 00089 *> \verbatim 00090 *> TL is REAL array, dimension (LDTL,2) 00091 *> On entry, TL contains an N1 by N1 matrix. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] LDTL 00095 *> \verbatim 00096 *> LDTL is INTEGER 00097 *> The leading dimension of the matrix TL. LDTL >= max(1,N1). 00098 *> \endverbatim 00099 *> 00100 *> \param[in] TR 00101 *> \verbatim 00102 *> TR is REAL array, dimension (LDTR,2) 00103 *> On entry, TR contains an N2 by N2 matrix. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDTR 00107 *> \verbatim 00108 *> LDTR is INTEGER 00109 *> The leading dimension of the matrix TR. LDTR >= max(1,N2). 00110 *> \endverbatim 00111 *> 00112 *> \param[in] B 00113 *> \verbatim 00114 *> B is REAL array, dimension (LDB,2) 00115 *> On entry, the N1 by N2 matrix B contains the right-hand 00116 *> side of the equation. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] LDB 00120 *> \verbatim 00121 *> LDB is INTEGER 00122 *> The leading dimension of the matrix B. LDB >= max(1,N1). 00123 *> \endverbatim 00124 *> 00125 *> \param[out] SCALE 00126 *> \verbatim 00127 *> SCALE is REAL 00128 *> On exit, SCALE contains the scale factor. SCALE is chosen 00129 *> less than or equal to 1 to prevent the solution overflowing. 00130 *> \endverbatim 00131 *> 00132 *> \param[out] X 00133 *> \verbatim 00134 *> X is REAL array, dimension (LDX,2) 00135 *> On exit, X contains the N1 by N2 solution. 00136 *> \endverbatim 00137 *> 00138 *> \param[in] LDX 00139 *> \verbatim 00140 *> LDX is INTEGER 00141 *> The leading dimension of the matrix X. LDX >= max(1,N1). 00142 *> \endverbatim 00143 *> 00144 *> \param[out] XNORM 00145 *> \verbatim 00146 *> XNORM is REAL 00147 *> On exit, XNORM is the infinity-norm of the solution. 00148 *> \endverbatim 00149 *> 00150 *> \param[out] INFO 00151 *> \verbatim 00152 *> INFO is INTEGER 00153 *> On exit, INFO is set to 00154 *> 0: successful exit. 00155 *> 1: TL and TR have too close eigenvalues, so TL or 00156 *> TR is perturbed to get a nonsingular equation. 00157 *> NOTE: In the interests of speed, this routine does not 00158 *> check the inputs for errors. 00159 *> \endverbatim 00160 * 00161 * Authors: 00162 * ======== 00163 * 00164 *> \author Univ. of Tennessee 00165 *> \author Univ. of California Berkeley 00166 *> \author Univ. of Colorado Denver 00167 *> \author NAG Ltd. 00168 * 00169 *> \date November 2011 00170 * 00171 *> \ingroup realSYauxiliary 00172 * 00173 * ===================================================================== 00174 SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, 00175 $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO ) 00176 * 00177 * -- LAPACK auxiliary routine (version 3.4.0) -- 00178 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00180 * November 2011 00181 * 00182 * .. Scalar Arguments .. 00183 LOGICAL LTRANL, LTRANR 00184 INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2 00185 REAL SCALE, XNORM 00186 * .. 00187 * .. Array Arguments .. 00188 REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ), 00189 $ X( LDX, * ) 00190 * .. 00191 * 00192 * ===================================================================== 00193 * 00194 * .. Parameters .. 00195 REAL ZERO, ONE 00196 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00197 REAL TWO, HALF, EIGHT 00198 PARAMETER ( TWO = 2.0E+0, HALF = 0.5E+0, EIGHT = 8.0E+0 ) 00199 * .. 00200 * .. Local Scalars .. 00201 LOGICAL BSWAP, XSWAP 00202 INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K 00203 REAL BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1, 00204 $ TEMP, U11, U12, U22, XMAX 00205 * .. 00206 * .. Local Arrays .. 00207 LOGICAL BSWPIV( 4 ), XSWPIV( 4 ) 00208 INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ), 00209 $ LOCU22( 4 ) 00210 REAL BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 ) 00211 * .. 00212 * .. External Functions .. 00213 INTEGER ISAMAX 00214 REAL SLAMCH 00215 EXTERNAL ISAMAX, SLAMCH 00216 * .. 00217 * .. External Subroutines .. 00218 EXTERNAL SCOPY, SSWAP 00219 * .. 00220 * .. Intrinsic Functions .. 00221 INTRINSIC ABS, MAX 00222 * .. 00223 * .. Data statements .. 00224 DATA LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / , 00225 $ LOCU22 / 4, 3, 2, 1 / 00226 DATA XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. / 00227 DATA BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. / 00228 * .. 00229 * .. Executable Statements .. 00230 * 00231 * Do not check the input parameters for errors 00232 * 00233 INFO = 0 00234 * 00235 * Quick return if possible 00236 * 00237 IF( N1.EQ.0 .OR. N2.EQ.0 ) 00238 $ RETURN 00239 * 00240 * Set constants to control overflow 00241 * 00242 EPS = SLAMCH( 'P' ) 00243 SMLNUM = SLAMCH( 'S' ) / EPS 00244 SGN = ISGN 00245 * 00246 K = N1 + N1 + N2 - 2 00247 GO TO ( 10, 20, 30, 50 )K 00248 * 00249 * 1 by 1: TL11*X + SGN*X*TR11 = B11 00250 * 00251 10 CONTINUE 00252 TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 ) 00253 BET = ABS( TAU1 ) 00254 IF( BET.LE.SMLNUM ) THEN 00255 TAU1 = SMLNUM 00256 BET = SMLNUM 00257 INFO = 1 00258 END IF 00259 * 00260 SCALE = ONE 00261 GAM = ABS( B( 1, 1 ) ) 00262 IF( SMLNUM*GAM.GT.BET ) 00263 $ SCALE = ONE / GAM 00264 * 00265 X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1 00266 XNORM = ABS( X( 1, 1 ) ) 00267 RETURN 00268 * 00269 * 1 by 2: 00270 * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12] 00271 * [TR21 TR22] 00272 * 00273 20 CONTINUE 00274 * 00275 SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ), 00276 $ ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ), 00277 $ SMLNUM ) 00278 TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 ) 00279 TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 ) 00280 IF( LTRANR ) THEN 00281 TMP( 2 ) = SGN*TR( 2, 1 ) 00282 TMP( 3 ) = SGN*TR( 1, 2 ) 00283 ELSE 00284 TMP( 2 ) = SGN*TR( 1, 2 ) 00285 TMP( 3 ) = SGN*TR( 2, 1 ) 00286 END IF 00287 BTMP( 1 ) = B( 1, 1 ) 00288 BTMP( 2 ) = B( 1, 2 ) 00289 GO TO 40 00290 * 00291 * 2 by 1: 00292 * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11] 00293 * [TL21 TL22] [X21] [X21] [B21] 00294 * 00295 30 CONTINUE 00296 SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ), 00297 $ ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ), 00298 $ SMLNUM ) 00299 TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 ) 00300 TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 ) 00301 IF( LTRANL ) THEN 00302 TMP( 2 ) = TL( 1, 2 ) 00303 TMP( 3 ) = TL( 2, 1 ) 00304 ELSE 00305 TMP( 2 ) = TL( 2, 1 ) 00306 TMP( 3 ) = TL( 1, 2 ) 00307 END IF 00308 BTMP( 1 ) = B( 1, 1 ) 00309 BTMP( 2 ) = B( 2, 1 ) 00310 40 CONTINUE 00311 * 00312 * Solve 2 by 2 system using complete pivoting. 00313 * Set pivots less than SMIN to SMIN. 00314 * 00315 IPIV = ISAMAX( 4, TMP, 1 ) 00316 U11 = TMP( IPIV ) 00317 IF( ABS( U11 ).LE.SMIN ) THEN 00318 INFO = 1 00319 U11 = SMIN 00320 END IF 00321 U12 = TMP( LOCU12( IPIV ) ) 00322 L21 = TMP( LOCL21( IPIV ) ) / U11 00323 U22 = TMP( LOCU22( IPIV ) ) - U12*L21 00324 XSWAP = XSWPIV( IPIV ) 00325 BSWAP = BSWPIV( IPIV ) 00326 IF( ABS( U22 ).LE.SMIN ) THEN 00327 INFO = 1 00328 U22 = SMIN 00329 END IF 00330 IF( BSWAP ) THEN 00331 TEMP = BTMP( 2 ) 00332 BTMP( 2 ) = BTMP( 1 ) - L21*TEMP 00333 BTMP( 1 ) = TEMP 00334 ELSE 00335 BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 ) 00336 END IF 00337 SCALE = ONE 00338 IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR. 00339 $ ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN 00340 SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) ) 00341 BTMP( 1 ) = BTMP( 1 )*SCALE 00342 BTMP( 2 ) = BTMP( 2 )*SCALE 00343 END IF 00344 X2( 2 ) = BTMP( 2 ) / U22 00345 X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 ) 00346 IF( XSWAP ) THEN 00347 TEMP = X2( 2 ) 00348 X2( 2 ) = X2( 1 ) 00349 X2( 1 ) = TEMP 00350 END IF 00351 X( 1, 1 ) = X2( 1 ) 00352 IF( N1.EQ.1 ) THEN 00353 X( 1, 2 ) = X2( 2 ) 00354 XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) 00355 ELSE 00356 X( 2, 1 ) = X2( 2 ) 00357 XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) ) 00358 END IF 00359 RETURN 00360 * 00361 * 2 by 2: 00362 * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12] 00363 * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22] 00364 * 00365 * Solve equivalent 4 by 4 system using complete pivoting. 00366 * Set pivots less than SMIN to SMIN. 00367 * 00368 50 CONTINUE 00369 SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ), 00370 $ ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ) 00371 SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ), 00372 $ ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ) 00373 SMIN = MAX( EPS*SMIN, SMLNUM ) 00374 BTMP( 1 ) = ZERO 00375 CALL SCOPY( 16, BTMP, 0, T16, 1 ) 00376 T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 ) 00377 T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 ) 00378 T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 ) 00379 T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 ) 00380 IF( LTRANL ) THEN 00381 T16( 1, 2 ) = TL( 2, 1 ) 00382 T16( 2, 1 ) = TL( 1, 2 ) 00383 T16( 3, 4 ) = TL( 2, 1 ) 00384 T16( 4, 3 ) = TL( 1, 2 ) 00385 ELSE 00386 T16( 1, 2 ) = TL( 1, 2 ) 00387 T16( 2, 1 ) = TL( 2, 1 ) 00388 T16( 3, 4 ) = TL( 1, 2 ) 00389 T16( 4, 3 ) = TL( 2, 1 ) 00390 END IF 00391 IF( LTRANR ) THEN 00392 T16( 1, 3 ) = SGN*TR( 1, 2 ) 00393 T16( 2, 4 ) = SGN*TR( 1, 2 ) 00394 T16( 3, 1 ) = SGN*TR( 2, 1 ) 00395 T16( 4, 2 ) = SGN*TR( 2, 1 ) 00396 ELSE 00397 T16( 1, 3 ) = SGN*TR( 2, 1 ) 00398 T16( 2, 4 ) = SGN*TR( 2, 1 ) 00399 T16( 3, 1 ) = SGN*TR( 1, 2 ) 00400 T16( 4, 2 ) = SGN*TR( 1, 2 ) 00401 END IF 00402 BTMP( 1 ) = B( 1, 1 ) 00403 BTMP( 2 ) = B( 2, 1 ) 00404 BTMP( 3 ) = B( 1, 2 ) 00405 BTMP( 4 ) = B( 2, 2 ) 00406 * 00407 * Perform elimination 00408 * 00409 DO 100 I = 1, 3 00410 XMAX = ZERO 00411 DO 70 IP = I, 4 00412 DO 60 JP = I, 4 00413 IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN 00414 XMAX = ABS( T16( IP, JP ) ) 00415 IPSV = IP 00416 JPSV = JP 00417 END IF 00418 60 CONTINUE 00419 70 CONTINUE 00420 IF( IPSV.NE.I ) THEN 00421 CALL SSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 ) 00422 TEMP = BTMP( I ) 00423 BTMP( I ) = BTMP( IPSV ) 00424 BTMP( IPSV ) = TEMP 00425 END IF 00426 IF( JPSV.NE.I ) 00427 $ CALL SSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 ) 00428 JPIV( I ) = JPSV 00429 IF( ABS( T16( I, I ) ).LT.SMIN ) THEN 00430 INFO = 1 00431 T16( I, I ) = SMIN 00432 END IF 00433 DO 90 J = I + 1, 4 00434 T16( J, I ) = T16( J, I ) / T16( I, I ) 00435 BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I ) 00436 DO 80 K = I + 1, 4 00437 T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K ) 00438 80 CONTINUE 00439 90 CONTINUE 00440 100 CONTINUE 00441 IF( ABS( T16( 4, 4 ) ).LT.SMIN ) 00442 $ T16( 4, 4 ) = SMIN 00443 SCALE = ONE 00444 IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR. 00445 $ ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR. 00446 $ ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR. 00447 $ ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN 00448 SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ), 00449 $ ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) ) 00450 BTMP( 1 ) = BTMP( 1 )*SCALE 00451 BTMP( 2 ) = BTMP( 2 )*SCALE 00452 BTMP( 3 ) = BTMP( 3 )*SCALE 00453 BTMP( 4 ) = BTMP( 4 )*SCALE 00454 END IF 00455 DO 120 I = 1, 4 00456 K = 5 - I 00457 TEMP = ONE / T16( K, K ) 00458 TMP( K ) = BTMP( K )*TEMP 00459 DO 110 J = K + 1, 4 00460 TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J ) 00461 110 CONTINUE 00462 120 CONTINUE 00463 DO 130 I = 1, 3 00464 IF( JPIV( 4-I ).NE.4-I ) THEN 00465 TEMP = TMP( 4-I ) 00466 TMP( 4-I ) = TMP( JPIV( 4-I ) ) 00467 TMP( JPIV( 4-I ) ) = TEMP 00468 END IF 00469 130 CONTINUE 00470 X( 1, 1 ) = TMP( 1 ) 00471 X( 2, 1 ) = TMP( 2 ) 00472 X( 1, 2 ) = TMP( 3 ) 00473 X( 2, 2 ) = TMP( 4 ) 00474 XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ), 00475 $ ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) ) 00476 RETURN 00477 * 00478 * End of SLASY2 00479 * 00480 END