![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DSYSVX computes the solution to system of linear equations A * X = B for SY matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSYSVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsysvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsysvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsysvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00022 * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00023 * IWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER FACT, UPLO 00027 * INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00028 * DOUBLE PRECISION RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IPIV( * ), IWORK( * ) 00032 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00033 * $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> DSYSVX uses the diagonal pivoting factorization to compute the 00043 *> solution to a real system of linear equations A * X = B, 00044 *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS 00045 *> matrices. 00046 *> 00047 *> Error bounds on the solution and a condition estimate are also 00048 *> provided. 00049 *> \endverbatim 00050 * 00051 *> \par Description: 00052 * ================= 00053 *> 00054 *> \verbatim 00055 *> 00056 *> The following steps are performed: 00057 *> 00058 *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A. 00059 *> The form of the factorization is 00060 *> A = U * D * U**T, if UPLO = 'U', or 00061 *> A = L * D * L**T, if UPLO = 'L', 00062 *> where U (or L) is a product of permutation and unit upper (lower) 00063 *> triangular matrices, and D is symmetric and block diagonal with 00064 *> 1-by-1 and 2-by-2 diagonal blocks. 00065 *> 00066 *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine 00067 *> returns with INFO = i. Otherwise, the factored form of A is used 00068 *> to estimate the condition number of the matrix A. If the 00069 *> reciprocal of the condition number is less than machine precision, 00070 *> INFO = N+1 is returned as a warning, but the routine still goes on 00071 *> to solve for X and compute error bounds as described below. 00072 *> 00073 *> 3. The system of equations is solved for X using the factored form 00074 *> of A. 00075 *> 00076 *> 4. Iterative refinement is applied to improve the computed solution 00077 *> matrix and calculate error bounds and backward error estimates 00078 *> for it. 00079 *> \endverbatim 00080 * 00081 * Arguments: 00082 * ========== 00083 * 00084 *> \param[in] FACT 00085 *> \verbatim 00086 *> FACT is CHARACTER*1 00087 *> Specifies whether or not the factored form of A has been 00088 *> supplied on entry. 00089 *> = 'F': On entry, AF and IPIV contain the factored form of 00090 *> A. AF and IPIV will not be modified. 00091 *> = 'N': The matrix A will be copied to AF and factored. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] UPLO 00095 *> \verbatim 00096 *> UPLO is CHARACTER*1 00097 *> = 'U': Upper triangle of A is stored; 00098 *> = 'L': Lower triangle of A is stored. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] N 00102 *> \verbatim 00103 *> N is INTEGER 00104 *> The number of linear equations, i.e., the order of the 00105 *> matrix A. N >= 0. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] NRHS 00109 *> \verbatim 00110 *> NRHS is INTEGER 00111 *> The number of right hand sides, i.e., the number of columns 00112 *> of the matrices B and X. NRHS >= 0. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] A 00116 *> \verbatim 00117 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00118 *> The symmetric matrix A. If UPLO = 'U', the leading N-by-N 00119 *> upper triangular part of A contains the upper triangular part 00120 *> of the matrix A, and the strictly lower triangular part of A 00121 *> is not referenced. If UPLO = 'L', the leading N-by-N lower 00122 *> triangular part of A contains the lower triangular part of 00123 *> the matrix A, and the strictly upper triangular part of A is 00124 *> not referenced. 00125 *> \endverbatim 00126 *> 00127 *> \param[in] LDA 00128 *> \verbatim 00129 *> LDA is INTEGER 00130 *> The leading dimension of the array A. LDA >= max(1,N). 00131 *> \endverbatim 00132 *> 00133 *> \param[in,out] AF 00134 *> \verbatim 00135 *> AF is DOUBLE PRECISION array, dimension (LDAF,N) 00136 *> If FACT = 'F', then AF is an input argument and on entry 00137 *> contains the block diagonal matrix D and the multipliers used 00138 *> to obtain the factor U or L from the factorization 00139 *> A = U*D*U**T or A = L*D*L**T as computed by DSYTRF. 00140 *> 00141 *> If FACT = 'N', then AF is an output argument and on exit 00142 *> returns the block diagonal matrix D and the multipliers used 00143 *> to obtain the factor U or L from the factorization 00144 *> A = U*D*U**T or A = L*D*L**T. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] LDAF 00148 *> \verbatim 00149 *> LDAF is INTEGER 00150 *> The leading dimension of the array AF. LDAF >= max(1,N). 00151 *> \endverbatim 00152 *> 00153 *> \param[in,out] IPIV 00154 *> \verbatim 00155 *> IPIV is INTEGER array, dimension (N) 00156 *> If FACT = 'F', then IPIV is an input argument and on entry 00157 *> contains details of the interchanges and the block structure 00158 *> of D, as determined by DSYTRF. 00159 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00160 *> interchanged and D(k,k) is a 1-by-1 diagonal block. 00161 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00162 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00163 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00164 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00165 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00166 *> 00167 *> If FACT = 'N', then IPIV is an output argument and on exit 00168 *> contains details of the interchanges and the block structure 00169 *> of D, as determined by DSYTRF. 00170 *> \endverbatim 00171 *> 00172 *> \param[in] B 00173 *> \verbatim 00174 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00175 *> The N-by-NRHS right hand side matrix B. 00176 *> \endverbatim 00177 *> 00178 *> \param[in] LDB 00179 *> \verbatim 00180 *> LDB is INTEGER 00181 *> The leading dimension of the array B. LDB >= max(1,N). 00182 *> \endverbatim 00183 *> 00184 *> \param[out] X 00185 *> \verbatim 00186 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00187 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 00188 *> \endverbatim 00189 *> 00190 *> \param[in] LDX 00191 *> \verbatim 00192 *> LDX is INTEGER 00193 *> The leading dimension of the array X. LDX >= max(1,N). 00194 *> \endverbatim 00195 *> 00196 *> \param[out] RCOND 00197 *> \verbatim 00198 *> RCOND is DOUBLE PRECISION 00199 *> The estimate of the reciprocal condition number of the matrix 00200 *> A. If RCOND is less than the machine precision (in 00201 *> particular, if RCOND = 0), the matrix is singular to working 00202 *> precision. This condition is indicated by a return code of 00203 *> INFO > 0. 00204 *> \endverbatim 00205 *> 00206 *> \param[out] FERR 00207 *> \verbatim 00208 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00209 *> The estimated forward error bound for each solution vector 00210 *> X(j) (the j-th column of the solution matrix X). 00211 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00212 *> is an estimated upper bound for the magnitude of the largest 00213 *> element in (X(j) - XTRUE) divided by the magnitude of the 00214 *> largest element in X(j). The estimate is as reliable as 00215 *> the estimate for RCOND, and is almost always a slight 00216 *> overestimate of the true error. 00217 *> \endverbatim 00218 *> 00219 *> \param[out] BERR 00220 *> \verbatim 00221 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00222 *> The componentwise relative backward error of each solution 00223 *> vector X(j) (i.e., the smallest relative change in 00224 *> any element of A or B that makes X(j) an exact solution). 00225 *> \endverbatim 00226 *> 00227 *> \param[out] WORK 00228 *> \verbatim 00229 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00230 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00231 *> \endverbatim 00232 *> 00233 *> \param[in] LWORK 00234 *> \verbatim 00235 *> LWORK is INTEGER 00236 *> The length of WORK. LWORK >= max(1,3*N), and for best 00237 *> performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where 00238 *> NB is the optimal blocksize for DSYTRF. 00239 *> 00240 *> If LWORK = -1, then a workspace query is assumed; the routine 00241 *> only calculates the optimal size of the WORK array, returns 00242 *> this value as the first entry of the WORK array, and no error 00243 *> message related to LWORK is issued by XERBLA. 00244 *> \endverbatim 00245 *> 00246 *> \param[out] IWORK 00247 *> \verbatim 00248 *> IWORK is INTEGER array, dimension (N) 00249 *> \endverbatim 00250 *> 00251 *> \param[out] INFO 00252 *> \verbatim 00253 *> INFO is INTEGER 00254 *> = 0: successful exit 00255 *> < 0: if INFO = -i, the i-th argument had an illegal value 00256 *> > 0: if INFO = i, and i is 00257 *> <= N: D(i,i) is exactly zero. The factorization 00258 *> has been completed but the factor D is exactly 00259 *> singular, so the solution and error bounds could 00260 *> not be computed. RCOND = 0 is returned. 00261 *> = N+1: D is nonsingular, but RCOND is less than machine 00262 *> precision, meaning that the matrix is singular 00263 *> to working precision. Nevertheless, the 00264 *> solution and error bounds are computed because 00265 *> there are a number of situations where the 00266 *> computed solution can be more accurate than the 00267 *> value of RCOND would suggest. 00268 *> \endverbatim 00269 * 00270 * Authors: 00271 * ======== 00272 * 00273 *> \author Univ. of Tennessee 00274 *> \author Univ. of California Berkeley 00275 *> \author Univ. of Colorado Denver 00276 *> \author NAG Ltd. 00277 * 00278 *> \date April 2012 00279 * 00280 *> \ingroup doubleSYsolve 00281 * 00282 * ===================================================================== 00283 SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00284 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00285 $ IWORK, INFO ) 00286 * 00287 * -- LAPACK driver routine (version 3.4.1) -- 00288 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00290 * April 2012 00291 * 00292 * .. Scalar Arguments .. 00293 CHARACTER FACT, UPLO 00294 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00295 DOUBLE PRECISION RCOND 00296 * .. 00297 * .. Array Arguments .. 00298 INTEGER IPIV( * ), IWORK( * ) 00299 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00300 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 00301 * .. 00302 * 00303 * ===================================================================== 00304 * 00305 * .. Parameters .. 00306 DOUBLE PRECISION ZERO 00307 PARAMETER ( ZERO = 0.0D+0 ) 00308 * .. 00309 * .. Local Scalars .. 00310 LOGICAL LQUERY, NOFACT 00311 INTEGER LWKOPT, NB 00312 DOUBLE PRECISION ANORM 00313 * .. 00314 * .. External Functions .. 00315 LOGICAL LSAME 00316 INTEGER ILAENV 00317 DOUBLE PRECISION DLAMCH, DLANSY 00318 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY 00319 * .. 00320 * .. External Subroutines .. 00321 EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA 00322 * .. 00323 * .. Intrinsic Functions .. 00324 INTRINSIC MAX 00325 * .. 00326 * .. Executable Statements .. 00327 * 00328 * Test the input parameters. 00329 * 00330 INFO = 0 00331 NOFACT = LSAME( FACT, 'N' ) 00332 LQUERY = ( LWORK.EQ.-1 ) 00333 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00334 INFO = -1 00335 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00336 $ THEN 00337 INFO = -2 00338 ELSE IF( N.LT.0 ) THEN 00339 INFO = -3 00340 ELSE IF( NRHS.LT.0 ) THEN 00341 INFO = -4 00342 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00343 INFO = -6 00344 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00345 INFO = -8 00346 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00347 INFO = -11 00348 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00349 INFO = -13 00350 ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN 00351 INFO = -18 00352 END IF 00353 * 00354 IF( INFO.EQ.0 ) THEN 00355 LWKOPT = MAX( 1, 3*N ) 00356 IF( NOFACT ) THEN 00357 NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) 00358 LWKOPT = MAX( LWKOPT, N*NB ) 00359 END IF 00360 WORK( 1 ) = LWKOPT 00361 END IF 00362 * 00363 IF( INFO.NE.0 ) THEN 00364 CALL XERBLA( 'DSYSVX', -INFO ) 00365 RETURN 00366 ELSE IF( LQUERY ) THEN 00367 RETURN 00368 END IF 00369 * 00370 IF( NOFACT ) THEN 00371 * 00372 * Compute the factorization A = U*D*U**T or A = L*D*L**T. 00373 * 00374 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00375 CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 00376 * 00377 * Return if INFO is non-zero. 00378 * 00379 IF( INFO.GT.0 )THEN 00380 RCOND = ZERO 00381 RETURN 00382 END IF 00383 END IF 00384 * 00385 * Compute the norm of the matrix A. 00386 * 00387 ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) 00388 * 00389 * Compute the reciprocal of the condition number of A. 00390 * 00391 CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, 00392 $ INFO ) 00393 * 00394 * Compute the solution vectors X. 00395 * 00396 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00397 CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00398 * 00399 * Use iterative refinement to improve the computed solutions and 00400 * compute error bounds and backward error estimates for them. 00401 * 00402 CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00403 $ LDX, FERR, BERR, WORK, IWORK, INFO ) 00404 * 00405 * Set INFO = N+1 if the matrix is singular to working precision. 00406 * 00407 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00408 $ INFO = N + 1 00409 * 00410 WORK( 1 ) = LWKOPT 00411 * 00412 RETURN 00413 * 00414 * End of DSYSVX 00415 * 00416 END