 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00002 $ WORK, INFO ) 00003 * 00004 * -- LAPACK PROTOTYPE routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * -- Written by Julie Langou of the Univ. of TN -- 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER UPLO 00013 INTEGER INFO, LDA, LDB, N, NRHS 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IPIV( * ) 00017 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DSYTRS2 solves a system of linear equations A*X = B with a real 00024 * symmetric matrix A using the factorization A = U*D*U**T or 00025 * A = L*D*L**T computed by DSYTRF and converted by DSYCONV. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * Specifies whether the details of the factorization are stored 00032 * as an upper or lower triangular matrix. 00033 * = 'U': Upper triangular, form is A = U*D*U**T; 00034 * = 'L': Lower triangular, form is A = L*D*L**T. 00035 * 00036 * N (input) INTEGER 00037 * The order of the matrix A. N >= 0. 00038 * 00039 * NRHS (input) INTEGER 00040 * The number of right hand sides, i.e., the number of columns 00041 * of the matrix B. NRHS >= 0. 00042 * 00043 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00044 * The block diagonal matrix D and the multipliers used to 00045 * obtain the factor U or L as computed by DSYTRF. 00046 * 00047 * LDA (input) INTEGER 00048 * The leading dimension of the array A. LDA >= max(1,N). 00049 * 00050 * IPIV (input) INTEGER array, dimension (N) 00051 * Details of the interchanges and the block structure of D 00052 * as determined by DSYTRF. 00053 * 00054 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00055 * On entry, the right hand side matrix B. 00056 * On exit, the solution matrix X. 00057 * 00058 * LDB (input) INTEGER 00059 * The leading dimension of the array B. LDB >= max(1,N). 00060 * 00061 * WORK (workspace) REAL array, dimension (N) 00062 * 00063 * INFO (output) INTEGER 00064 * = 0: successful exit 00065 * < 0: if INFO = -i, the i-th argument had an illegal value 00066 * 00067 * ===================================================================== 00068 * 00069 * .. Parameters .. 00070 DOUBLE PRECISION ONE 00071 PARAMETER ( ONE = 1.0D+0 ) 00072 * .. 00073 * .. Local Scalars .. 00074 LOGICAL UPPER 00075 INTEGER I, IINFO, J, K, KP 00076 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM 00077 * .. 00078 * .. External Functions .. 00079 LOGICAL LSAME 00080 EXTERNAL LSAME 00081 * .. 00082 * .. External Subroutines .. 00083 EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA 00084 * .. 00085 * .. Intrinsic Functions .. 00086 INTRINSIC MAX 00087 * .. 00088 * .. Executable Statements .. 00089 * 00090 INFO = 0 00091 UPPER = LSAME( UPLO, 'U' ) 00092 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00093 INFO = -1 00094 ELSE IF( N.LT.0 ) THEN 00095 INFO = -2 00096 ELSE IF( NRHS.LT.0 ) THEN 00097 INFO = -3 00098 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00099 INFO = -5 00100 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00101 INFO = -8 00102 END IF 00103 IF( INFO.NE.0 ) THEN 00104 CALL XERBLA( 'DSYTRS2', -INFO ) 00105 RETURN 00106 END IF 00107 * 00108 * Quick return if possible 00109 * 00110 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00111 $ RETURN 00112 * 00113 * Convert A 00114 * 00115 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00116 * 00117 IF( UPPER ) THEN 00118 * 00119 * Solve A*X = B, where A = U*D*U**T. 00120 * 00121 * P**T * B 00122 K=N 00123 DO WHILE ( K .GE. 1 ) 00124 IF( IPIV( K ).GT.0 ) THEN 00125 * 1 x 1 diagonal block 00126 * Interchange rows K and IPIV(K). 00127 KP = IPIV( K ) 00128 IF( KP.NE.K ) 00129 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00130 K=K-1 00131 ELSE 00132 * 2 x 2 diagonal block 00133 * Interchange rows K-1 and -IPIV(K). 00134 KP = -IPIV( K ) 00135 IF( KP.EQ.-IPIV( K-1 ) ) 00136 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00137 K=K-2 00138 END IF 00139 END DO 00140 * 00141 * Compute (U \P**T * B) -> B [ (U \P**T * B) ] 00142 * 00143 CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00144 * 00145 * Compute D \ B -> B [ D \ (U \P**T * B) ] 00146 * 00147 I=N 00148 DO WHILE ( I .GE. 1 ) 00149 IF( IPIV(I) .GT. 0 ) THEN 00150 CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB ) 00151 ELSEIF ( I .GT. 1) THEN 00152 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 00153 AKM1K = WORK(I) 00154 AKM1 = A( I-1, I-1 ) / AKM1K 00155 AK = A( I, I ) / AKM1K 00156 DENOM = AKM1*AK - ONE 00157 DO 15 J = 1, NRHS 00158 BKM1 = B( I-1, J ) / AKM1K 00159 BK = B( I, J ) / AKM1K 00160 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 00161 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 00162 15 CONTINUE 00163 I = I - 1 00164 ENDIF 00165 ENDIF 00166 I = I - 1 00167 END DO 00168 * 00169 * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ] 00170 * 00171 CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB) 00172 * 00173 * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ] 00174 * 00175 K=1 00176 DO WHILE ( K .LE. N ) 00177 IF( IPIV( K ).GT.0 ) THEN 00178 * 1 x 1 diagonal block 00179 * Interchange rows K and IPIV(K). 00180 KP = IPIV( K ) 00181 IF( KP.NE.K ) 00182 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00183 K=K+1 00184 ELSE 00185 * 2 x 2 diagonal block 00186 * Interchange rows K-1 and -IPIV(K). 00187 KP = -IPIV( K ) 00188 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 00189 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00190 K=K+2 00191 ENDIF 00192 END DO 00193 * 00194 ELSE 00195 * 00196 * Solve A*X = B, where A = L*D*L**T. 00197 * 00198 * P**T * B 00199 K=1 00200 DO WHILE ( K .LE. N ) 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 and -IPIV(K+1). 00211 KP = -IPIV( K+1 ) 00212 IF( KP.EQ.-IPIV( K ) ) 00213 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00214 K=K+2 00215 ENDIF 00216 END DO 00217 * 00218 * Compute (L \P**T * B) -> B [ (L \P**T * B) ] 00219 * 00220 CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00221 * 00222 * Compute D \ B -> B [ D \ (L \P**T * B) ] 00223 * 00224 I=1 00225 DO WHILE ( I .LE. N ) 00226 IF( IPIV(I) .GT. 0 ) THEN 00227 CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB ) 00228 ELSE 00229 AKM1K = WORK(I) 00230 AKM1 = A( I, I ) / AKM1K 00231 AK = A( I+1, I+1 ) / AKM1K 00232 DENOM = AKM1*AK - ONE 00233 DO 25 J = 1, NRHS 00234 BKM1 = B( I, J ) / AKM1K 00235 BK = B( I+1, J ) / AKM1K 00236 B( I, J ) = ( AK*BKM1-BK ) / DENOM 00237 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00238 25 CONTINUE 00239 I = I + 1 00240 ENDIF 00241 I = I + 1 00242 END DO 00243 * 00244 * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ] 00245 * 00246 CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB) 00247 * 00248 * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ] 00249 * 00250 K=N 00251 DO WHILE ( K .GE. 1 ) 00252 IF( IPIV( K ).GT.0 ) THEN 00253 * 1 x 1 diagonal block 00254 * Interchange rows K and IPIV(K). 00255 KP = IPIV( K ) 00256 IF( KP.NE.K ) 00257 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00258 K=K-1 00259 ELSE 00260 * 2 x 2 diagonal block 00261 * Interchange rows K-1 and -IPIV(K). 00262 KP = -IPIV( K ) 00263 IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) ) 00264 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00265 K=K-2 00266 ENDIF 00267 END DO 00268 * 00269 END IF 00270 * 00271 * Revert A 00272 * 00273 CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 00274 * 00275 RETURN 00276 * 00277 * End of DSYTRS2 00278 * 00279 END
1.7.3 
