![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DPOSVX computes the solution to system of linear equations A * X = B for PO 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 DPOSVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dposvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dposvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dposvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 00022 * S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 00023 * IWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER EQUED, FACT, UPLO 00027 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00028 * DOUBLE PRECISION RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IWORK( * ) 00032 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00033 * $ BERR( * ), FERR( * ), S( * ), WORK( * ), 00034 * $ X( LDX, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to 00044 *> compute the solution to a real system of linear equations 00045 *> A * X = B, 00046 *> where A is an N-by-N symmetric positive definite matrix and X and B 00047 *> are N-by-NRHS matrices. 00048 *> 00049 *> Error bounds on the solution and a condition estimate are also 00050 *> provided. 00051 *> \endverbatim 00052 * 00053 *> \par Description: 00054 * ================= 00055 *> 00056 *> \verbatim 00057 *> 00058 *> The following steps are performed: 00059 *> 00060 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate 00061 *> the system: 00062 *> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B 00063 *> Whether or not the system will be equilibrated depends on the 00064 *> scaling of the matrix A, but if equilibration is used, A is 00065 *> overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 00066 *> 00067 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to 00068 *> factor the matrix A (after equilibration if FACT = 'E') as 00069 *> A = U**T* U, if UPLO = 'U', or 00070 *> A = L * L**T, if UPLO = 'L', 00071 *> where U is an upper triangular matrix and L is a lower triangular 00072 *> matrix. 00073 *> 00074 *> 3. If the leading i-by-i principal minor is not positive definite, 00075 *> then the routine returns with INFO = i. Otherwise, the factored 00076 *> form of A is used to estimate the condition number of the matrix 00077 *> A. If the reciprocal of the condition number is less than machine 00078 *> precision, INFO = N+1 is returned as a warning, but the routine 00079 *> still goes on to solve for X and compute error bounds as 00080 *> described below. 00081 *> 00082 *> 4. The system of equations is solved for X using the factored form 00083 *> of A. 00084 *> 00085 *> 5. Iterative refinement is applied to improve the computed solution 00086 *> matrix and calculate error bounds and backward error estimates 00087 *> for it. 00088 *> 00089 *> 6. If equilibration was used, the matrix X is premultiplied by 00090 *> diag(S) so that it solves the original system before 00091 *> equilibration. 00092 *> \endverbatim 00093 * 00094 * Arguments: 00095 * ========== 00096 * 00097 *> \param[in] FACT 00098 *> \verbatim 00099 *> FACT is CHARACTER*1 00100 *> Specifies whether or not the factored form of the matrix A is 00101 *> supplied on entry, and if not, whether the matrix A should be 00102 *> equilibrated before it is factored. 00103 *> = 'F': On entry, AF contains the factored form of A. 00104 *> If EQUED = 'Y', the matrix A has been equilibrated 00105 *> with scaling factors given by S. A and AF will not 00106 *> be modified. 00107 *> = 'N': The matrix A will be copied to AF and factored. 00108 *> = 'E': The matrix A will be equilibrated if necessary, then 00109 *> copied to AF and factored. 00110 *> \endverbatim 00111 *> 00112 *> \param[in] UPLO 00113 *> \verbatim 00114 *> UPLO is CHARACTER*1 00115 *> = 'U': Upper triangle of A is stored; 00116 *> = 'L': Lower triangle of A is stored. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] N 00120 *> \verbatim 00121 *> N is INTEGER 00122 *> The number of linear equations, i.e., the order of the 00123 *> matrix A. N >= 0. 00124 *> \endverbatim 00125 *> 00126 *> \param[in] NRHS 00127 *> \verbatim 00128 *> NRHS is INTEGER 00129 *> The number of right hand sides, i.e., the number of columns 00130 *> of the matrices B and X. NRHS >= 0. 00131 *> \endverbatim 00132 *> 00133 *> \param[in,out] A 00134 *> \verbatim 00135 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00136 *> On entry, the symmetric matrix A, except if FACT = 'F' and 00137 *> EQUED = 'Y', then A must contain the equilibrated matrix 00138 *> diag(S)*A*diag(S). If UPLO = 'U', the leading 00139 *> N-by-N upper triangular part of A contains the upper 00140 *> triangular part of the matrix A, and the strictly lower 00141 *> triangular part of A is not referenced. If UPLO = 'L', the 00142 *> leading N-by-N lower triangular part of A contains the lower 00143 *> triangular part of the matrix A, and the strictly upper 00144 *> triangular part of A is not referenced. A is not modified if 00145 *> FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. 00146 *> 00147 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 00148 *> diag(S)*A*diag(S). 00149 *> \endverbatim 00150 *> 00151 *> \param[in] LDA 00152 *> \verbatim 00153 *> LDA is INTEGER 00154 *> The leading dimension of the array A. LDA >= max(1,N). 00155 *> \endverbatim 00156 *> 00157 *> \param[in,out] AF 00158 *> \verbatim 00159 *> AF is DOUBLE PRECISION array, dimension (LDAF,N) 00160 *> If FACT = 'F', then AF is an input argument and on entry 00161 *> contains the triangular factor U or L from the Cholesky 00162 *> factorization A = U**T*U or A = L*L**T, in the same storage 00163 *> format as A. If EQUED .ne. 'N', then AF is the factored form 00164 *> of the equilibrated matrix diag(S)*A*diag(S). 00165 *> 00166 *> If FACT = 'N', then AF is an output argument and on exit 00167 *> returns the triangular factor U or L from the Cholesky 00168 *> factorization A = U**T*U or A = L*L**T of the original 00169 *> matrix A. 00170 *> 00171 *> If FACT = 'E', then AF is an output argument and on exit 00172 *> returns the triangular factor U or L from the Cholesky 00173 *> factorization A = U**T*U or A = L*L**T of the equilibrated 00174 *> matrix A (see the description of A for the form of the 00175 *> equilibrated matrix). 00176 *> \endverbatim 00177 *> 00178 *> \param[in] LDAF 00179 *> \verbatim 00180 *> LDAF is INTEGER 00181 *> The leading dimension of the array AF. LDAF >= max(1,N). 00182 *> \endverbatim 00183 *> 00184 *> \param[in,out] EQUED 00185 *> \verbatim 00186 *> EQUED is CHARACTER*1 00187 *> Specifies the form of equilibration that was done. 00188 *> = 'N': No equilibration (always true if FACT = 'N'). 00189 *> = 'Y': Equilibration was done, i.e., A has been replaced by 00190 *> diag(S) * A * diag(S). 00191 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an 00192 *> output argument. 00193 *> \endverbatim 00194 *> 00195 *> \param[in,out] S 00196 *> \verbatim 00197 *> S is DOUBLE PRECISION array, dimension (N) 00198 *> The scale factors for A; not accessed if EQUED = 'N'. S is 00199 *> an input argument if FACT = 'F'; otherwise, S is an output 00200 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S 00201 *> must be positive. 00202 *> \endverbatim 00203 *> 00204 *> \param[in,out] B 00205 *> \verbatim 00206 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00207 *> On entry, the N-by-NRHS right hand side matrix B. 00208 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', 00209 *> B is overwritten by diag(S) * B. 00210 *> \endverbatim 00211 *> 00212 *> \param[in] LDB 00213 *> \verbatim 00214 *> LDB is INTEGER 00215 *> The leading dimension of the array B. LDB >= max(1,N). 00216 *> \endverbatim 00217 *> 00218 *> \param[out] X 00219 *> \verbatim 00220 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00221 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to 00222 *> the original system of equations. Note that if EQUED = 'Y', 00223 *> A and B are modified on exit, and the solution to the 00224 *> equilibrated system is inv(diag(S))*X. 00225 *> \endverbatim 00226 *> 00227 *> \param[in] LDX 00228 *> \verbatim 00229 *> LDX is INTEGER 00230 *> The leading dimension of the array X. LDX >= max(1,N). 00231 *> \endverbatim 00232 *> 00233 *> \param[out] RCOND 00234 *> \verbatim 00235 *> RCOND is DOUBLE PRECISION 00236 *> The estimate of the reciprocal condition number of the matrix 00237 *> A after equilibration (if done). If RCOND is less than the 00238 *> machine precision (in particular, if RCOND = 0), the matrix 00239 *> is singular to working precision. This condition is 00240 *> indicated by a return code of INFO > 0. 00241 *> \endverbatim 00242 *> 00243 *> \param[out] FERR 00244 *> \verbatim 00245 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00246 *> The estimated forward error bound for each solution vector 00247 *> X(j) (the j-th column of the solution matrix X). 00248 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00249 *> is an estimated upper bound for the magnitude of the largest 00250 *> element in (X(j) - XTRUE) divided by the magnitude of the 00251 *> largest element in X(j). The estimate is as reliable as 00252 *> the estimate for RCOND, and is almost always a slight 00253 *> overestimate of the true error. 00254 *> \endverbatim 00255 *> 00256 *> \param[out] BERR 00257 *> \verbatim 00258 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00259 *> The componentwise relative backward error of each solution 00260 *> vector X(j) (i.e., the smallest relative change in 00261 *> any element of A or B that makes X(j) an exact solution). 00262 *> \endverbatim 00263 *> 00264 *> \param[out] WORK 00265 *> \verbatim 00266 *> WORK is DOUBLE PRECISION array, dimension (3*N) 00267 *> \endverbatim 00268 *> 00269 *> \param[out] IWORK 00270 *> \verbatim 00271 *> IWORK is INTEGER array, dimension (N) 00272 *> \endverbatim 00273 *> 00274 *> \param[out] INFO 00275 *> \verbatim 00276 *> INFO is INTEGER 00277 *> = 0: successful exit 00278 *> < 0: if INFO = -i, the i-th argument had an illegal value 00279 *> > 0: if INFO = i, and i is 00280 *> <= N: the leading minor of order i of A is 00281 *> not positive definite, so the factorization 00282 *> could not be completed, and the solution has not 00283 *> been computed. RCOND = 0 is returned. 00284 *> = N+1: U is nonsingular, but RCOND is less than machine 00285 *> precision, meaning that the matrix is singular 00286 *> to working precision. Nevertheless, the 00287 *> solution and error bounds are computed because 00288 *> there are a number of situations where the 00289 *> computed solution can be more accurate than the 00290 *> value of RCOND would suggest. 00291 *> \endverbatim 00292 * 00293 * Authors: 00294 * ======== 00295 * 00296 *> \author Univ. of Tennessee 00297 *> \author Univ. of California Berkeley 00298 *> \author Univ. of Colorado Denver 00299 *> \author NAG Ltd. 00300 * 00301 *> \date April 2012 00302 * 00303 *> \ingroup doublePOsolve 00304 * 00305 * ===================================================================== 00306 SUBROUTINE DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 00307 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 00308 $ IWORK, INFO ) 00309 * 00310 * -- LAPACK driver routine (version 3.4.1) -- 00311 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00312 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00313 * April 2012 00314 * 00315 * .. Scalar Arguments .. 00316 CHARACTER EQUED, FACT, UPLO 00317 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00318 DOUBLE PRECISION RCOND 00319 * .. 00320 * .. Array Arguments .. 00321 INTEGER IWORK( * ) 00322 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00323 $ BERR( * ), FERR( * ), S( * ), WORK( * ), 00324 $ X( LDX, * ) 00325 * .. 00326 * 00327 * ===================================================================== 00328 * 00329 * .. Parameters .. 00330 DOUBLE PRECISION ZERO, ONE 00331 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00332 * .. 00333 * .. Local Scalars .. 00334 LOGICAL EQUIL, NOFACT, RCEQU 00335 INTEGER I, INFEQU, J 00336 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 00337 * .. 00338 * .. External Functions .. 00339 LOGICAL LSAME 00340 DOUBLE PRECISION DLAMCH, DLANSY 00341 EXTERNAL LSAME, DLAMCH, DLANSY 00342 * .. 00343 * .. External Subroutines .. 00344 EXTERNAL DLACPY, DLAQSY, DPOCON, DPOEQU, DPORFS, DPOTRF, 00345 $ DPOTRS, XERBLA 00346 * .. 00347 * .. Intrinsic Functions .. 00348 INTRINSIC MAX, MIN 00349 * .. 00350 * .. Executable Statements .. 00351 * 00352 INFO = 0 00353 NOFACT = LSAME( FACT, 'N' ) 00354 EQUIL = LSAME( FACT, 'E' ) 00355 IF( NOFACT .OR. EQUIL ) THEN 00356 EQUED = 'N' 00357 RCEQU = .FALSE. 00358 ELSE 00359 RCEQU = LSAME( EQUED, 'Y' ) 00360 SMLNUM = DLAMCH( 'Safe minimum' ) 00361 BIGNUM = ONE / SMLNUM 00362 END IF 00363 * 00364 * Test the input parameters. 00365 * 00366 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 00367 $ THEN 00368 INFO = -1 00369 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00370 $ THEN 00371 INFO = -2 00372 ELSE IF( N.LT.0 ) THEN 00373 INFO = -3 00374 ELSE IF( NRHS.LT.0 ) THEN 00375 INFO = -4 00376 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00377 INFO = -6 00378 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00379 INFO = -8 00380 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 00381 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 00382 INFO = -9 00383 ELSE 00384 IF( RCEQU ) THEN 00385 SMIN = BIGNUM 00386 SMAX = ZERO 00387 DO 10 J = 1, N 00388 SMIN = MIN( SMIN, S( J ) ) 00389 SMAX = MAX( SMAX, S( J ) ) 00390 10 CONTINUE 00391 IF( SMIN.LE.ZERO ) THEN 00392 INFO = -10 00393 ELSE IF( N.GT.0 ) THEN 00394 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 00395 ELSE 00396 SCOND = ONE 00397 END IF 00398 END IF 00399 IF( INFO.EQ.0 ) THEN 00400 IF( LDB.LT.MAX( 1, N ) ) THEN 00401 INFO = -12 00402 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00403 INFO = -14 00404 END IF 00405 END IF 00406 END IF 00407 * 00408 IF( INFO.NE.0 ) THEN 00409 CALL XERBLA( 'DPOSVX', -INFO ) 00410 RETURN 00411 END IF 00412 * 00413 IF( EQUIL ) THEN 00414 * 00415 * Compute row and column scalings to equilibrate the matrix A. 00416 * 00417 CALL DPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 00418 IF( INFEQU.EQ.0 ) THEN 00419 * 00420 * Equilibrate the matrix. 00421 * 00422 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 00423 RCEQU = LSAME( EQUED, 'Y' ) 00424 END IF 00425 END IF 00426 * 00427 * Scale the right hand side. 00428 * 00429 IF( RCEQU ) THEN 00430 DO 30 J = 1, NRHS 00431 DO 20 I = 1, N 00432 B( I, J ) = S( I )*B( I, J ) 00433 20 CONTINUE 00434 30 CONTINUE 00435 END IF 00436 * 00437 IF( NOFACT .OR. EQUIL ) THEN 00438 * 00439 * Compute the Cholesky factorization A = U**T *U or A = L*L**T. 00440 * 00441 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00442 CALL DPOTRF( UPLO, N, AF, LDAF, INFO ) 00443 * 00444 * Return if INFO is non-zero. 00445 * 00446 IF( INFO.GT.0 )THEN 00447 RCOND = ZERO 00448 RETURN 00449 END IF 00450 END IF 00451 * 00452 * Compute the norm of the matrix A. 00453 * 00454 ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK ) 00455 * 00456 * Compute the reciprocal of the condition number of A. 00457 * 00458 CALL DPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO ) 00459 * 00460 * Compute the solution matrix X. 00461 * 00462 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00463 CALL DPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 00464 * 00465 * Use iterative refinement to improve the computed solution and 00466 * compute error bounds and backward error estimates for it. 00467 * 00468 CALL DPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 00469 $ FERR, BERR, WORK, IWORK, INFO ) 00470 * 00471 * Transform the solution matrix X to a solution of the original 00472 * system. 00473 * 00474 IF( RCEQU ) THEN 00475 DO 50 J = 1, NRHS 00476 DO 40 I = 1, N 00477 X( I, J ) = S( I )*X( I, J ) 00478 40 CONTINUE 00479 50 CONTINUE 00480 DO 60 J = 1, NRHS 00481 FERR( J ) = FERR( J ) / SCOND 00482 60 CONTINUE 00483 END IF 00484 * 00485 * Set INFO = N+1 if the matrix is singular to working precision. 00486 * 00487 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00488 $ INFO = N + 1 00489 * 00490 RETURN 00491 * 00492 * End of DPOSVX 00493 * 00494 END