![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSYTRS2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSYTRS2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrs2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrs2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrs2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00022 * WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, LDA, LDB, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IPIV( * ) 00030 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DSYTRS2 solves a system of linear equations A*X = B with a real 00040 *> symmetric matrix A using the factorization A = U*D*U**T or 00041 *> A = L*D*L**T computed by DSYTRF and converted by DSYCONV. 00042 *> \endverbatim 00043 * 00044 * Arguments: 00045 * ========== 00046 * 00047 *> \param[in] UPLO 00048 *> \verbatim 00049 *> UPLO is CHARACTER*1 00050 *> Specifies whether the details of the factorization are stored 00051 *> as an upper or lower triangular matrix. 00052 *> = 'U': Upper triangular, form is A = U*D*U**T; 00053 *> = 'L': Lower triangular, form is A = L*D*L**T. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] N 00057 *> \verbatim 00058 *> N is INTEGER 00059 *> The order of the matrix A. N >= 0. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] NRHS 00063 *> \verbatim 00064 *> NRHS is INTEGER 00065 *> The number of right hand sides, i.e., the number of columns 00066 *> of the matrix B. NRHS >= 0. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] A 00070 *> \verbatim 00071 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00072 *> The block diagonal matrix D and the multipliers used to 00073 *> obtain the factor U or L as computed by DSYTRF. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] LDA 00077 *> \verbatim 00078 *> LDA is INTEGER 00079 *> The leading dimension of the array A. LDA >= max(1,N). 00080 *> \endverbatim 00081 *> 00082 *> \param[in] IPIV 00083 *> \verbatim 00084 *> IPIV is INTEGER array, dimension (N) 00085 *> Details of the interchanges and the block structure of D 00086 *> as determined by DSYTRF. 00087 *> \endverbatim 00088 *> 00089 *> \param[in,out] B 00090 *> \verbatim 00091 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00092 *> On entry, the right hand side matrix B. 00093 *> On exit, the solution matrix X. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDB 00097 *> \verbatim 00098 *> LDB is INTEGER 00099 *> The leading dimension of the array B. LDB >= max(1,N). 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is REAL array, dimension (N) 00105 *> \endverbatim 00106 *> 00107 *> \param[out] INFO 00108 *> \verbatim 00109 *> INFO is INTEGER 00110 *> = 0: successful exit 00111 *> < 0: if INFO = -i, the i-th argument had an illegal value 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date November 2011 00123 * 00124 *> \ingroup doubleSYcomputational 00125 * 00126 * ===================================================================== 00127 SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00128 $ WORK, INFO ) 00129 * 00130 * -- LAPACK computational routine (version 3.4.0) -- 00131 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00133 * November 2011 00134 * 00135 * .. Scalar Arguments .. 00136 CHARACTER UPLO 00137 INTEGER INFO, LDA, LDB, N, NRHS 00138 * .. 00139 * .. Array Arguments .. 00140 INTEGER IPIV( * ) 00141 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00142 * .. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. Parameters .. 00147 DOUBLE PRECISION ONE 00148 PARAMETER ( ONE = 1.0D+0 ) 00149 * .. 00150 * .. Local Scalars .. 00151 LOGICAL UPPER 00152 INTEGER I, IINFO, J, K, KP 00153 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM 00154 * .. 00155 * .. External Functions .. 00156 LOGICAL LSAME 00157 EXTERNAL LSAME 00158 * .. 00159 * .. External Subroutines .. 00160 EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA 00161 * .. 00162 * .. Intrinsic Functions .. 00163 INTRINSIC MAX 00164 * .. 00165 * .. Executable Statements .. 00166 * 00167 INFO = 0 00168 UPPER = LSAME( UPLO, 'U' ) 00169 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00170 INFO = -1 00171 ELSE IF( N.LT.0 ) THEN 00172 INFO = -2 00173 ELSE IF( NRHS.LT.0 ) THEN 00174 INFO = -3 00175 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00176 INFO = -5 00177 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00178 INFO = -8 00179 END IF 00180 IF( INFO.NE.0 ) THEN 00181 CALL XERBLA( 'DSYTRS2', -INFO ) 00182 RETURN 00183 END IF 00184 * 00185 * Quick return if possible 00186 * 00187 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00188 $ RETURN 00189 * 00190 * Convert A 00191 * 00192 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00193 * 00194 IF( UPPER ) THEN 00195 * 00196 * Solve A*X = B, where A = U*D*U**T. 00197 * 00198 * P**T * B 00199 K=N 00200 DO WHILE ( K .GE. 1 ) 00201 IF( IPIV( K ).GT.0 ) THEN 00202 * 1 x 1 diagonal block 00203 * Interchange rows K and IPIV(K). 00204 KP = IPIV( K ) 00205 IF( KP.NE.K ) 00206 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00207 K=K-1 00208 ELSE 00209 * 2 x 2 diagonal block 00210 * Interchange rows K-1 and -IPIV(K). 00211 KP = -IPIV( K ) 00212 IF( KP.EQ.-IPIV( K-1 ) ) 00213 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00214 K=K-2 00215 END IF 00216 END DO 00217 * 00218 * Compute (U \P**T * B) -> B [ (U \P**T * B) ] 00219 * 00220 CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00221 * 00222 * Compute D \ B -> B [ D \ (U \P**T * B) ] 00223 * 00224 I=N 00225 DO WHILE ( I .GE. 1 ) 00226 IF( IPIV(I) .GT. 0 ) THEN 00227 CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB ) 00228 ELSEIF ( I .GT. 1) THEN 00229 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 00230 AKM1K = WORK(I) 00231 AKM1 = A( I-1, I-1 ) / AKM1K 00232 AK = A( I, I ) / AKM1K 00233 DENOM = AKM1*AK - ONE 00234 DO 15 J = 1, NRHS 00235 BKM1 = B( I-1, J ) / AKM1K 00236 BK = B( I, J ) / AKM1K 00237 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 00238 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 00239 15 CONTINUE 00240 I = I - 1 00241 ENDIF 00242 ENDIF 00243 I = I - 1 00244 END DO 00245 * 00246 * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ] 00247 * 00248 CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB) 00249 * 00250 * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ] 00251 * 00252 K=1 00253 DO WHILE ( K .LE. N ) 00254 IF( IPIV( K ).GT.0 ) THEN 00255 * 1 x 1 diagonal block 00256 * Interchange rows K and IPIV(K). 00257 KP = IPIV( K ) 00258 IF( KP.NE.K ) 00259 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00260 K=K+1 00261 ELSE 00262 * 2 x 2 diagonal block 00263 * Interchange rows K-1 and -IPIV(K). 00264 KP = -IPIV( K ) 00265 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 00266 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00267 K=K+2 00268 ENDIF 00269 END DO 00270 * 00271 ELSE 00272 * 00273 * Solve A*X = B, where A = L*D*L**T. 00274 * 00275 * P**T * B 00276 K=1 00277 DO WHILE ( K .LE. N ) 00278 IF( IPIV( K ).GT.0 ) THEN 00279 * 1 x 1 diagonal block 00280 * Interchange rows K and IPIV(K). 00281 KP = IPIV( K ) 00282 IF( KP.NE.K ) 00283 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00284 K=K+1 00285 ELSE 00286 * 2 x 2 diagonal block 00287 * Interchange rows K and -IPIV(K+1). 00288 KP = -IPIV( K+1 ) 00289 IF( KP.EQ.-IPIV( K ) ) 00290 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00291 K=K+2 00292 ENDIF 00293 END DO 00294 * 00295 * Compute (L \P**T * B) -> B [ (L \P**T * B) ] 00296 * 00297 CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00298 * 00299 * Compute D \ B -> B [ D \ (L \P**T * B) ] 00300 * 00301 I=1 00302 DO WHILE ( I .LE. N ) 00303 IF( IPIV(I) .GT. 0 ) THEN 00304 CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB ) 00305 ELSE 00306 AKM1K = WORK(I) 00307 AKM1 = A( I, I ) / AKM1K 00308 AK = A( I+1, I+1 ) / AKM1K 00309 DENOM = AKM1*AK - ONE 00310 DO 25 J = 1, NRHS 00311 BKM1 = B( I, J ) / AKM1K 00312 BK = B( I+1, J ) / AKM1K 00313 B( I, J ) = ( AK*BKM1-BK ) / DENOM 00314 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00315 25 CONTINUE 00316 I = I + 1 00317 ENDIF 00318 I = I + 1 00319 END DO 00320 * 00321 * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ] 00322 * 00323 CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB) 00324 * 00325 * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ] 00326 * 00327 K=N 00328 DO WHILE ( K .GE. 1 ) 00329 IF( IPIV( K ).GT.0 ) THEN 00330 * 1 x 1 diagonal block 00331 * Interchange rows K and IPIV(K). 00332 KP = IPIV( K ) 00333 IF( KP.NE.K ) 00334 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00335 K=K-1 00336 ELSE 00337 * 2 x 2 diagonal block 00338 * Interchange rows K-1 and -IPIV(K). 00339 KP = -IPIV( K ) 00340 IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) ) 00341 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00342 K=K-2 00343 ENDIF 00344 END DO 00345 * 00346 END IF 00347 * 00348 * Revert A 00349 * 00350 CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 00351 * 00352 RETURN 00353 * 00354 * End of DSYTRS2 00355 * 00356 END