|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, LDA, LDB, N, NRHS 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ) 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DSYTRS solves a system of linear equations A*X = B with a real 00021 * symmetric matrix A using the factorization A = U*D*U**T or 00022 * A = L*D*L**T computed by DSYTRF. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * Specifies whether the details of the factorization are stored 00029 * as an upper or lower triangular matrix. 00030 * = 'U': Upper triangular, form is A = U*D*U**T; 00031 * = 'L': Lower triangular, form is A = L*D*L**T. 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. N >= 0. 00035 * 00036 * NRHS (input) INTEGER 00037 * The number of right hand sides, i.e., the number of columns 00038 * of the matrix B. NRHS >= 0. 00039 * 00040 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00041 * The block diagonal matrix D and the multipliers used to 00042 * obtain the factor U or L as computed by DSYTRF. 00043 * 00044 * LDA (input) INTEGER 00045 * The leading dimension of the array A. LDA >= max(1,N). 00046 * 00047 * IPIV (input) INTEGER array, dimension (N) 00048 * Details of the interchanges and the block structure of D 00049 * as determined by DSYTRF. 00050 * 00051 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00052 * On entry, the right hand side matrix B. 00053 * On exit, the solution matrix X. 00054 * 00055 * LDB (input) INTEGER 00056 * The leading dimension of the array B. LDB >= max(1,N). 00057 * 00058 * INFO (output) INTEGER 00059 * = 0: successful exit 00060 * < 0: if INFO = -i, the i-th argument had an illegal value 00061 * 00062 * ===================================================================== 00063 * 00064 * .. Parameters .. 00065 DOUBLE PRECISION ONE 00066 PARAMETER ( ONE = 1.0D+0 ) 00067 * .. 00068 * .. Local Scalars .. 00069 LOGICAL UPPER 00070 INTEGER J, K, KP 00071 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM 00072 * .. 00073 * .. External Functions .. 00074 LOGICAL LSAME 00075 EXTERNAL LSAME 00076 * .. 00077 * .. External Subroutines .. 00078 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA 00079 * .. 00080 * .. Intrinsic Functions .. 00081 INTRINSIC MAX 00082 * .. 00083 * .. Executable Statements .. 00084 * 00085 INFO = 0 00086 UPPER = LSAME( UPLO, 'U' ) 00087 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00088 INFO = -1 00089 ELSE IF( N.LT.0 ) THEN 00090 INFO = -2 00091 ELSE IF( NRHS.LT.0 ) THEN 00092 INFO = -3 00093 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00094 INFO = -5 00095 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00096 INFO = -8 00097 END IF 00098 IF( INFO.NE.0 ) THEN 00099 CALL XERBLA( 'DSYTRS', -INFO ) 00100 RETURN 00101 END IF 00102 * 00103 * Quick return if possible 00104 * 00105 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00106 $ RETURN 00107 * 00108 IF( UPPER ) THEN 00109 * 00110 * Solve A*X = B, where A = U*D*U'. 00111 * 00112 * First solve U*D*X = B, overwriting B with X. 00113 * 00114 * K is the main loop index, decreasing from N to 1 in steps of 00115 * 1 or 2, depending on the size of the diagonal blocks. 00116 * 00117 K = N 00118 10 CONTINUE 00119 * 00120 * If K < 1, exit from loop. 00121 * 00122 IF( K.LT.1 ) 00123 $ GO TO 30 00124 * 00125 IF( IPIV( K ).GT.0 ) THEN 00126 * 00127 * 1 x 1 diagonal block 00128 * 00129 * Interchange rows K and IPIV(K). 00130 * 00131 KP = IPIV( K ) 00132 IF( KP.NE.K ) 00133 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00134 * 00135 * Multiply by inv(U(K)), where U(K) is the transformation 00136 * stored in column K of A. 00137 * 00138 CALL DGER( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, 00139 $ B( 1, 1 ), LDB ) 00140 * 00141 * Multiply by the inverse of the diagonal block. 00142 * 00143 CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) 00144 K = K - 1 00145 ELSE 00146 * 00147 * 2 x 2 diagonal block 00148 * 00149 * Interchange rows K-1 and -IPIV(K). 00150 * 00151 KP = -IPIV( K ) 00152 IF( KP.NE.K-1 ) 00153 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00154 * 00155 * Multiply by inv(U(K)), where U(K) is the transformation 00156 * stored in columns K-1 and K of A. 00157 * 00158 CALL DGER( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, 00159 $ B( 1, 1 ), LDB ) 00160 CALL DGER( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ), 00161 $ LDB, B( 1, 1 ), LDB ) 00162 * 00163 * Multiply by the inverse of the diagonal block. 00164 * 00165 AKM1K = A( K-1, K ) 00166 AKM1 = A( K-1, K-1 ) / AKM1K 00167 AK = A( K, K ) / AKM1K 00168 DENOM = AKM1*AK - ONE 00169 DO 20 J = 1, NRHS 00170 BKM1 = B( K-1, J ) / AKM1K 00171 BK = B( K, J ) / AKM1K 00172 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM 00173 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM 00174 20 CONTINUE 00175 K = K - 2 00176 END IF 00177 * 00178 GO TO 10 00179 30 CONTINUE 00180 * 00181 * Next solve U'*X = B, overwriting B with X. 00182 * 00183 * K is the main loop index, increasing from 1 to N in steps of 00184 * 1 or 2, depending on the size of the diagonal blocks. 00185 * 00186 K = 1 00187 40 CONTINUE 00188 * 00189 * If K > N, exit from loop. 00190 * 00191 IF( K.GT.N ) 00192 $ GO TO 50 00193 * 00194 IF( IPIV( K ).GT.0 ) THEN 00195 * 00196 * 1 x 1 diagonal block 00197 * 00198 * Multiply by inv(U'(K)), where U(K) is the transformation 00199 * stored in column K of A. 00200 * 00201 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), 00202 $ 1, ONE, B( K, 1 ), LDB ) 00203 * 00204 * Interchange rows K and IPIV(K). 00205 * 00206 KP = IPIV( K ) 00207 IF( KP.NE.K ) 00208 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00209 K = K + 1 00210 ELSE 00211 * 00212 * 2 x 2 diagonal block 00213 * 00214 * Multiply by inv(U'(K+1)), where U(K+1) is the transformation 00215 * stored in columns K and K+1 of A. 00216 * 00217 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), 00218 $ 1, ONE, B( K, 1 ), LDB ) 00219 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, 00220 $ A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB ) 00221 * 00222 * Interchange rows K and -IPIV(K). 00223 * 00224 KP = -IPIV( K ) 00225 IF( KP.NE.K ) 00226 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00227 K = K + 2 00228 END IF 00229 * 00230 GO TO 40 00231 50 CONTINUE 00232 * 00233 ELSE 00234 * 00235 * Solve A*X = B, where A = L*D*L'. 00236 * 00237 * First solve L*D*X = B, overwriting B with X. 00238 * 00239 * K is the main loop index, increasing from 1 to N in steps of 00240 * 1 or 2, depending on the size of the diagonal blocks. 00241 * 00242 K = 1 00243 60 CONTINUE 00244 * 00245 * If K > N, exit from loop. 00246 * 00247 IF( K.GT.N ) 00248 $ GO TO 80 00249 * 00250 IF( IPIV( K ).GT.0 ) THEN 00251 * 00252 * 1 x 1 diagonal block 00253 * 00254 * Interchange rows K and IPIV(K). 00255 * 00256 KP = IPIV( K ) 00257 IF( KP.NE.K ) 00258 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00259 * 00260 * Multiply by inv(L(K)), where L(K) is the transformation 00261 * stored in column K of A. 00262 * 00263 IF( K.LT.N ) 00264 $ CALL DGER( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ), 00265 $ LDB, B( K+1, 1 ), LDB ) 00266 * 00267 * Multiply by the inverse of the diagonal block. 00268 * 00269 CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) 00270 K = K + 1 00271 ELSE 00272 * 00273 * 2 x 2 diagonal block 00274 * 00275 * Interchange rows K+1 and -IPIV(K). 00276 * 00277 KP = -IPIV( K ) 00278 IF( KP.NE.K+1 ) 00279 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00280 * 00281 * Multiply by inv(L(K)), where L(K) is the transformation 00282 * stored in columns K and K+1 of A. 00283 * 00284 IF( K.LT.N-1 ) THEN 00285 CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ), 00286 $ LDB, B( K+2, 1 ), LDB ) 00287 CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1, 00288 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) 00289 END IF 00290 * 00291 * Multiply by the inverse of the diagonal block. 00292 * 00293 AKM1K = A( K+1, K ) 00294 AKM1 = A( K, K ) / AKM1K 00295 AK = A( K+1, K+1 ) / AKM1K 00296 DENOM = AKM1*AK - ONE 00297 DO 70 J = 1, NRHS 00298 BKM1 = B( K, J ) / AKM1K 00299 BK = B( K+1, J ) / AKM1K 00300 B( K, J ) = ( AK*BKM1-BK ) / DENOM 00301 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00302 70 CONTINUE 00303 K = K + 2 00304 END IF 00305 * 00306 GO TO 60 00307 80 CONTINUE 00308 * 00309 * Next solve L'*X = B, overwriting B with X. 00310 * 00311 * K is the main loop index, decreasing from N to 1 in steps of 00312 * 1 or 2, depending on the size of the diagonal blocks. 00313 * 00314 K = N 00315 90 CONTINUE 00316 * 00317 * If K < 1, exit from loop. 00318 * 00319 IF( K.LT.1 ) 00320 $ GO TO 100 00321 * 00322 IF( IPIV( K ).GT.0 ) THEN 00323 * 00324 * 1 x 1 diagonal block 00325 * 00326 * Multiply by inv(L'(K)), where L(K) is the transformation 00327 * stored in column K of A. 00328 * 00329 IF( K.LT.N ) 00330 $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 00331 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 00332 * 00333 * Interchange rows K and IPIV(K). 00334 * 00335 KP = IPIV( K ) 00336 IF( KP.NE.K ) 00337 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00338 K = K - 1 00339 ELSE 00340 * 00341 * 2 x 2 diagonal block 00342 * 00343 * Multiply by inv(L'(K-1)), where L(K-1) is the transformation 00344 * stored in columns K-1 and K of A. 00345 * 00346 IF( K.LT.N ) THEN 00347 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 00348 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 00349 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), 00350 $ LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ), 00351 $ LDB ) 00352 END IF 00353 * 00354 * Interchange rows K and -IPIV(K). 00355 * 00356 KP = -IPIV( K ) 00357 IF( KP.NE.K ) 00358 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00359 K = K - 2 00360 END IF 00361 * 00362 GO TO 90 00363 100 CONTINUE 00364 END IF 00365 * 00366 RETURN 00367 * 00368 * End of DSYTRS 00369 * 00370 END
1.7.3 
