![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DPORFSX 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DPORFSX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dporfsx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dporfsx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dporfsx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B, 00022 * LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, 00023 * ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, 00024 * WORK, IWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * CHARACTER UPLO, EQUED 00028 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 00029 * $ N_ERR_BNDS 00030 * DOUBLE PRECISION RCOND 00031 * .. 00032 * .. Array Arguments .. 00033 * INTEGER IWORK( * ) 00034 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00035 * $ X( LDX, * ), WORK( * ) 00036 * DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), 00037 * $ ERR_BNDS_NORM( NRHS, * ), 00038 * $ ERR_BNDS_COMP( NRHS, * ) 00039 * .. 00040 * 00041 * 00042 *> \par Purpose: 00043 * ============= 00044 *> 00045 *> \verbatim 00046 *> 00047 *> DPORFSX improves the computed solution to a system of linear 00048 *> equations when the coefficient matrix is symmetric positive 00049 *> definite, and provides error bounds and backward error estimates 00050 *> for the solution. In addition to normwise error bound, the code 00051 *> provides maximum componentwise error bound if possible. See 00052 *> comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the 00053 *> error bounds. 00054 *> 00055 *> The original system of linear equations may have been equilibrated 00056 *> before calling this routine, as described by arguments EQUED and S 00057 *> below. In this case, the solution and error bounds returned are 00058 *> for the original unequilibrated system. 00059 *> \endverbatim 00060 * 00061 * Arguments: 00062 * ========== 00063 * 00064 *> \verbatim 00065 *> Some optional parameters are bundled in the PARAMS array. These 00066 *> settings determine how refinement is performed, but often the 00067 *> defaults are acceptable. If the defaults are acceptable, users 00068 *> can pass NPARAMS = 0 which prevents the source code from accessing 00069 *> the PARAMS argument. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] UPLO 00073 *> \verbatim 00074 *> UPLO is CHARACTER*1 00075 *> = 'U': Upper triangle of A is stored; 00076 *> = 'L': Lower triangle of A is stored. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] EQUED 00080 *> \verbatim 00081 *> EQUED is CHARACTER*1 00082 *> Specifies the form of equilibration that was done to A 00083 *> before calling this routine. This is needed to compute 00084 *> the solution and error bounds correctly. 00085 *> = 'N': No equilibration 00086 *> = 'Y': Both row and column equilibration, i.e., A has been 00087 *> replaced by diag(S) * A * diag(S). 00088 *> The right hand side B has been changed accordingly. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] N 00092 *> \verbatim 00093 *> N is INTEGER 00094 *> The order of the matrix A. N >= 0. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] NRHS 00098 *> \verbatim 00099 *> NRHS is INTEGER 00100 *> The number of right hand sides, i.e., the number of columns 00101 *> of the matrices B and X. NRHS >= 0. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] A 00105 *> \verbatim 00106 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00107 *> The symmetric matrix A. If UPLO = 'U', the leading N-by-N 00108 *> upper triangular part of A contains the upper triangular part 00109 *> of the matrix A, and the strictly lower triangular part of A 00110 *> is not referenced. If UPLO = 'L', the leading N-by-N lower 00111 *> triangular part of A contains the lower triangular part of 00112 *> the matrix A, and the strictly upper triangular part of A is 00113 *> not referenced. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDA 00117 *> \verbatim 00118 *> LDA is INTEGER 00119 *> The leading dimension of the array A. LDA >= max(1,N). 00120 *> \endverbatim 00121 *> 00122 *> \param[in] AF 00123 *> \verbatim 00124 *> AF is DOUBLE PRECISION array, dimension (LDAF,N) 00125 *> The triangular factor U or L from the Cholesky factorization 00126 *> A = U**T*U or A = L*L**T, as computed by DPOTRF. 00127 *> \endverbatim 00128 *> 00129 *> \param[in] LDAF 00130 *> \verbatim 00131 *> LDAF is INTEGER 00132 *> The leading dimension of the array AF. LDAF >= max(1,N). 00133 *> \endverbatim 00134 *> 00135 *> \param[in,out] S 00136 *> \verbatim 00137 *> S is DOUBLE PRECISION array, dimension (N) 00138 *> The row scale factors for A. If EQUED = 'Y', A is multiplied on 00139 *> the left and right by diag(S). S is an input argument if FACT = 00140 *> 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED 00141 *> = 'Y', each element of S must be positive. If S is output, each 00142 *> element of S is a power of the radix. If S is input, each element 00143 *> of S should be a power of the radix to ensure a reliable solution 00144 *> and error estimates. Scaling by powers of the radix does not cause 00145 *> rounding errors unless the result underflows or overflows. 00146 *> Rounding errors during scaling lead to refining with a matrix that 00147 *> is not equivalent to the input matrix, producing error estimates 00148 *> that may not be reliable. 00149 *> \endverbatim 00150 *> 00151 *> \param[in] B 00152 *> \verbatim 00153 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00154 *> The right hand side matrix B. 00155 *> \endverbatim 00156 *> 00157 *> \param[in] LDB 00158 *> \verbatim 00159 *> LDB is INTEGER 00160 *> The leading dimension of the array B. LDB >= max(1,N). 00161 *> \endverbatim 00162 *> 00163 *> \param[in,out] X 00164 *> \verbatim 00165 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00166 *> On entry, the solution matrix X, as computed by DGETRS. 00167 *> On exit, the improved solution matrix X. 00168 *> \endverbatim 00169 *> 00170 *> \param[in] LDX 00171 *> \verbatim 00172 *> LDX is INTEGER 00173 *> The leading dimension of the array X. LDX >= max(1,N). 00174 *> \endverbatim 00175 *> 00176 *> \param[out] RCOND 00177 *> \verbatim 00178 *> RCOND is DOUBLE PRECISION 00179 *> Reciprocal scaled condition number. This is an estimate of the 00180 *> reciprocal Skeel condition number of the matrix A after 00181 *> equilibration (if done). If this is less than the machine 00182 *> precision (in particular, if it is zero), the matrix is singular 00183 *> to working precision. Note that the error may still be small even 00184 *> if this number is very small and the matrix appears ill- 00185 *> conditioned. 00186 *> \endverbatim 00187 *> 00188 *> \param[out] BERR 00189 *> \verbatim 00190 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00191 *> Componentwise relative backward error. This is the 00192 *> componentwise relative backward error of each solution vector X(j) 00193 *> (i.e., the smallest relative change in any element of A or B that 00194 *> makes X(j) an exact solution). 00195 *> \endverbatim 00196 *> 00197 *> \param[in] N_ERR_BNDS 00198 *> \verbatim 00199 *> N_ERR_BNDS is INTEGER 00200 *> Number of error bounds to return for each right hand side 00201 *> and each type (normwise or componentwise). See ERR_BNDS_NORM and 00202 *> ERR_BNDS_COMP below. 00203 *> \endverbatim 00204 *> 00205 *> \param[out] ERR_BNDS_NORM 00206 *> \verbatim 00207 *> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) 00208 *> For each right-hand side, this array contains information about 00209 *> various error bounds and condition numbers corresponding to the 00210 *> normwise relative error, which is defined as follows: 00211 *> 00212 *> Normwise relative error in the ith solution vector: 00213 *> max_j (abs(XTRUE(j,i) - X(j,i))) 00214 *> ------------------------------ 00215 *> max_j abs(X(j,i)) 00216 *> 00217 *> The array is indexed by the type of error information as described 00218 *> below. There currently are up to three pieces of information 00219 *> returned. 00220 *> 00221 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 00222 *> right-hand side. 00223 *> 00224 *> The second index in ERR_BNDS_NORM(:,err) contains the following 00225 *> three fields: 00226 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the 00227 *> reciprocal condition number is less than the threshold 00228 *> sqrt(n) * dlamch('Epsilon'). 00229 *> 00230 *> err = 2 "Guaranteed" error bound: The estimated forward error, 00231 *> almost certainly within a factor of 10 of the true error 00232 *> so long as the next entry is greater than the threshold 00233 *> sqrt(n) * dlamch('Epsilon'). This error bound should only 00234 *> be trusted if the previous boolean is true. 00235 *> 00236 *> err = 3 Reciprocal condition number: Estimated normwise 00237 *> reciprocal condition number. Compared with the threshold 00238 *> sqrt(n) * dlamch('Epsilon') to determine if the error 00239 *> estimate is "guaranteed". These reciprocal condition 00240 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 00241 *> appropriately scaled matrix Z. 00242 *> Let Z = S*A, where S scales each row by a power of the 00243 *> radix so all absolute row sums of Z are approximately 1. 00244 *> 00245 *> See Lapack Working Note 165 for further details and extra 00246 *> cautions. 00247 *> \endverbatim 00248 *> 00249 *> \param[out] ERR_BNDS_COMP 00250 *> \verbatim 00251 *> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) 00252 *> For each right-hand side, this array contains information about 00253 *> various error bounds and condition numbers corresponding to the 00254 *> componentwise relative error, which is defined as follows: 00255 *> 00256 *> Componentwise relative error in the ith solution vector: 00257 *> abs(XTRUE(j,i) - X(j,i)) 00258 *> max_j ---------------------- 00259 *> abs(X(j,i)) 00260 *> 00261 *> The array is indexed by the right-hand side i (on which the 00262 *> componentwise relative error depends), and the type of error 00263 *> information as described below. There currently are up to three 00264 *> pieces of information returned for each right-hand side. If 00265 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then 00266 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most 00267 *> the first (:,N_ERR_BNDS) entries are returned. 00268 *> 00269 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith 00270 *> right-hand side. 00271 *> 00272 *> The second index in ERR_BNDS_COMP(:,err) contains the following 00273 *> three fields: 00274 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the 00275 *> reciprocal condition number is less than the threshold 00276 *> sqrt(n) * dlamch('Epsilon'). 00277 *> 00278 *> err = 2 "Guaranteed" error bound: The estimated forward error, 00279 *> almost certainly within a factor of 10 of the true error 00280 *> so long as the next entry is greater than the threshold 00281 *> sqrt(n) * dlamch('Epsilon'). This error bound should only 00282 *> be trusted if the previous boolean is true. 00283 *> 00284 *> err = 3 Reciprocal condition number: Estimated componentwise 00285 *> reciprocal condition number. Compared with the threshold 00286 *> sqrt(n) * dlamch('Epsilon') to determine if the error 00287 *> estimate is "guaranteed". These reciprocal condition 00288 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 00289 *> appropriately scaled matrix Z. 00290 *> Let Z = S*(A*diag(x)), where x is the solution for the 00291 *> current right-hand side and S scales each row of 00292 *> A*diag(x) by a power of the radix so all absolute row 00293 *> sums of Z are approximately 1. 00294 *> 00295 *> See Lapack Working Note 165 for further details and extra 00296 *> cautions. 00297 *> \endverbatim 00298 *> 00299 *> \param[in] NPARAMS 00300 *> \verbatim 00301 *> NPARAMS is INTEGER 00302 *> Specifies the number of parameters set in PARAMS. If .LE. 0, the 00303 *> PARAMS array is never referenced and default values are used. 00304 *> \endverbatim 00305 *> 00306 *> \param[in,out] PARAMS 00307 *> \verbatim 00308 *> PARAMS is / output) DOUBLE PRECISION array, dimension (NPARAMS) 00309 *> Specifies algorithm parameters. If an entry is .LT. 0.0, then 00310 *> that entry will be filled with default value used for that 00311 *> parameter. Only positions up to NPARAMS are accessed; defaults 00312 *> are used for higher-numbered parameters. 00313 *> 00314 *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative 00315 *> refinement or not. 00316 *> Default: 1.0D+0 00317 *> = 0.0 : No refinement is performed, and no error bounds are 00318 *> computed. 00319 *> = 1.0 : Use the double-precision refinement algorithm, 00320 *> possibly with doubled-single computations if the 00321 *> compilation environment does not support DOUBLE 00322 *> PRECISION. 00323 *> (other values are reserved for future use) 00324 *> 00325 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual 00326 *> computations allowed for refinement. 00327 *> Default: 10 00328 *> Aggressive: Set to 100 to permit convergence using approximate 00329 *> factorizations or factorizations other than LU. If 00330 *> the factorization uses a technique other than 00331 *> Gaussian elimination, the guarantees in 00332 *> err_bnds_norm and err_bnds_comp may no longer be 00333 *> trustworthy. 00334 *> 00335 *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code 00336 *> will attempt to find a solution with small componentwise 00337 *> relative error in the double-precision algorithm. Positive 00338 *> is true, 0.0 is false. 00339 *> Default: 1.0 (attempt componentwise convergence) 00340 *> \endverbatim 00341 *> 00342 *> \param[out] WORK 00343 *> \verbatim 00344 *> WORK is DOUBLE PRECISION array, dimension (4*N) 00345 *> \endverbatim 00346 *> 00347 *> \param[out] IWORK 00348 *> \verbatim 00349 *> IWORK is INTEGER array, dimension (N) 00350 *> \endverbatim 00351 *> 00352 *> \param[out] INFO 00353 *> \verbatim 00354 *> INFO is INTEGER 00355 *> = 0: Successful exit. The solution to every right-hand side is 00356 *> guaranteed. 00357 *> < 0: If INFO = -i, the i-th argument had an illegal value 00358 *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization 00359 *> has been completed, but the factor U is exactly singular, so 00360 *> the solution and error bounds could not be computed. RCOND = 0 00361 *> is returned. 00362 *> = N+J: The solution corresponding to the Jth right-hand side is 00363 *> not guaranteed. The solutions corresponding to other right- 00364 *> hand sides K with K > J may not be guaranteed as well, but 00365 *> only the first such right-hand side is reported. If a small 00366 *> componentwise error is not requested (PARAMS(3) = 0.0) then 00367 *> the Jth right-hand side is the first with a normwise error 00368 *> bound that is not guaranteed (the smallest J such 00369 *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) 00370 *> the Jth right-hand side is the first with either a normwise or 00371 *> componentwise error bound that is not guaranteed (the smallest 00372 *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or 00373 *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of 00374 *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information 00375 *> about all of the right-hand sides check ERR_BNDS_NORM or 00376 *> ERR_BNDS_COMP. 00377 *> \endverbatim 00378 * 00379 * Authors: 00380 * ======== 00381 * 00382 *> \author Univ. of Tennessee 00383 *> \author Univ. of California Berkeley 00384 *> \author Univ. of Colorado Denver 00385 *> \author NAG Ltd. 00386 * 00387 *> \date April 2012 00388 * 00389 *> \ingroup doublePOcomputational 00390 * 00391 * ===================================================================== 00392 SUBROUTINE DPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B, 00393 $ LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, 00394 $ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, 00395 $ WORK, IWORK, INFO ) 00396 * 00397 * -- LAPACK computational routine (version 3.4.1) -- 00398 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00399 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00400 * April 2012 00401 * 00402 * .. Scalar Arguments .. 00403 CHARACTER UPLO, EQUED 00404 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 00405 $ N_ERR_BNDS 00406 DOUBLE PRECISION RCOND 00407 * .. 00408 * .. Array Arguments .. 00409 INTEGER IWORK( * ) 00410 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00411 $ X( LDX, * ), WORK( * ) 00412 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), 00413 $ ERR_BNDS_NORM( NRHS, * ), 00414 $ ERR_BNDS_COMP( NRHS, * ) 00415 * .. 00416 * 00417 * ================================================================== 00418 * 00419 * .. Parameters .. 00420 DOUBLE PRECISION ZERO, ONE 00421 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00422 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT 00423 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT 00424 DOUBLE PRECISION DZTHRESH_DEFAULT 00425 PARAMETER ( ITREF_DEFAULT = 1.0D+0 ) 00426 PARAMETER ( ITHRESH_DEFAULT = 10.0D+0 ) 00427 PARAMETER ( COMPONENTWISE_DEFAULT = 1.0D+0 ) 00428 PARAMETER ( RTHRESH_DEFAULT = 0.5D+0 ) 00429 PARAMETER ( DZTHRESH_DEFAULT = 0.25D+0 ) 00430 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 00431 $ LA_LINRX_CWISE_I 00432 PARAMETER ( LA_LINRX_ITREF_I = 1, 00433 $ LA_LINRX_ITHRESH_I = 2 ) 00434 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 00435 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 00436 $ LA_LINRX_RCOND_I 00437 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 00438 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 00439 * .. 00440 * .. Local Scalars .. 00441 CHARACTER(1) NORM 00442 LOGICAL RCEQU 00443 INTEGER J, PREC_TYPE, REF_TYPE 00444 INTEGER N_NORMS 00445 DOUBLE PRECISION ANORM, RCOND_TMP 00446 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG 00447 LOGICAL IGNORE_CWISE 00448 INTEGER ITHRESH 00449 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH 00450 * .. 00451 * .. External Subroutines .. 00452 EXTERNAL XERBLA, DPOCON, DLA_PORFSX_EXTENDED 00453 * .. 00454 * .. Intrinsic Functions .. 00455 INTRINSIC MAX, SQRT 00456 * .. 00457 * .. External Functions .. 00458 EXTERNAL LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC 00459 EXTERNAL DLAMCH, DLANSY, DLA_PORCOND 00460 DOUBLE PRECISION DLAMCH, DLANSY, DLA_PORCOND 00461 LOGICAL LSAME 00462 INTEGER BLAS_FPINFO_X 00463 INTEGER ILATRANS, ILAPREC 00464 * .. 00465 * .. Executable Statements .. 00466 * 00467 * Check the input parameters. 00468 * 00469 INFO = 0 00470 REF_TYPE = INT( ITREF_DEFAULT ) 00471 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN 00472 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN 00473 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT 00474 ELSE 00475 REF_TYPE = PARAMS( LA_LINRX_ITREF_I ) 00476 END IF 00477 END IF 00478 * 00479 * Set default parameters. 00480 * 00481 ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' ) 00482 ITHRESH = INT( ITHRESH_DEFAULT ) 00483 RTHRESH = RTHRESH_DEFAULT 00484 UNSTABLE_THRESH = DZTHRESH_DEFAULT 00485 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0 00486 * 00487 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 00488 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN 00489 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH 00490 ELSE 00491 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) ) 00492 END IF 00493 END IF 00494 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN 00495 IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN 00496 IF ( IGNORE_CWISE ) THEN 00497 PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0 00498 ELSE 00499 PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0 00500 END IF 00501 ELSE 00502 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0 00503 END IF 00504 END IF 00505 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN 00506 N_NORMS = 0 00507 ELSE IF ( IGNORE_CWISE ) THEN 00508 N_NORMS = 1 00509 ELSE 00510 N_NORMS = 2 00511 END IF 00512 * 00513 RCEQU = LSAME( EQUED, 'Y' ) 00514 * 00515 * Test input parameters. 00516 * 00517 IF (.NOT.LSAME(UPLO, 'U') .AND. .NOT.LSAME(UPLO, 'L')) THEN 00518 INFO = -1 00519 ELSE IF( .NOT.RCEQU .AND. .NOT.LSAME( EQUED, 'N' ) ) THEN 00520 INFO = -2 00521 ELSE IF( N.LT.0 ) THEN 00522 INFO = -3 00523 ELSE IF( NRHS.LT.0 ) THEN 00524 INFO = -4 00525 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00526 INFO = -6 00527 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00528 INFO = -8 00529 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00530 INFO = -11 00531 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00532 INFO = -13 00533 END IF 00534 IF( INFO.NE.0 ) THEN 00535 CALL XERBLA( 'DPORFSX', -INFO ) 00536 RETURN 00537 END IF 00538 * 00539 * Quick return if possible. 00540 * 00541 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 00542 RCOND = 1.0D+0 00543 DO J = 1, NRHS 00544 BERR( J ) = 0.0D+0 00545 IF ( N_ERR_BNDS .GE. 1 ) THEN 00546 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00547 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00548 END IF 00549 IF ( N_ERR_BNDS .GE. 2 ) THEN 00550 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0 00551 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0 00552 END IF 00553 IF ( N_ERR_BNDS .GE. 3 ) THEN 00554 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0 00555 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0 00556 END IF 00557 END DO 00558 RETURN 00559 END IF 00560 * 00561 * Default to failure. 00562 * 00563 RCOND = 0.0D+0 00564 DO J = 1, NRHS 00565 BERR( J ) = 1.0D+0 00566 IF ( N_ERR_BNDS .GE. 1 ) THEN 00567 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00568 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00569 END IF 00570 IF ( N_ERR_BNDS .GE. 2 ) THEN 00571 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00572 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00573 END IF 00574 IF ( N_ERR_BNDS .GE. 3 ) THEN 00575 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0 00576 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0 00577 END IF 00578 END DO 00579 * 00580 * Compute the norm of A and the reciprocal of the condition 00581 * number of A. 00582 * 00583 NORM = 'I' 00584 ANORM = DLANSY( NORM, UPLO, N, A, LDA, WORK ) 00585 CALL DPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, 00586 $ IWORK, INFO ) 00587 * 00588 * Perform refinement on each right-hand side 00589 * 00590 IF ( REF_TYPE .NE. 0 ) THEN 00591 00592 PREC_TYPE = ILAPREC( 'E' ) 00593 00594 CALL DLA_PORFSX_EXTENDED( PREC_TYPE, UPLO, N, 00595 $ NRHS, A, LDA, AF, LDAF, RCEQU, S, B, 00596 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, 00597 $ WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ), WORK( 1 ), RCOND, 00598 $ ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 00599 $ INFO ) 00600 END IF 00601 00602 ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' ) 00603 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN 00604 * 00605 * Compute scaled normwise condition number cond(A*C). 00606 * 00607 IF ( RCEQU ) THEN 00608 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 00609 $ -1, S, INFO, WORK, IWORK ) 00610 ELSE 00611 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 00612 $ 0, S, INFO, WORK, IWORK ) 00613 END IF 00614 DO J = 1, NRHS 00615 * 00616 * Cap the error at 1.0. 00617 * 00618 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 00619 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 00620 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00621 * 00622 * Threshold the error (see LAWN). 00623 * 00624 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 00625 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00626 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0 00627 IF ( INFO .LE. N ) INFO = N + J 00628 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND ) 00629 $ THEN 00630 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 00631 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00632 END IF 00633 * 00634 * Save the condition number. 00635 * 00636 IF (N_ERR_BNDS .GE. LA_LINRX_RCOND_I) THEN 00637 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 00638 END IF 00639 END DO 00640 END IF 00641 00642 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2 ) THEN 00643 * 00644 * Compute componentwise condition number cond(A*diag(Y(:,J))) for 00645 * each right-hand side using the current solution as an estimate of 00646 * the true solution. If the componentwise error estimate is too 00647 * large, then the solution is a lousy estimate of truth and the 00648 * estimated RCOND may be too optimistic. To avoid misleading users, 00649 * the inverse condition number is set to 0.0 when the estimated 00650 * cwise error is at least CWISE_WRONG. 00651 * 00652 CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) ) 00653 DO J = 1, NRHS 00654 IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 00655 $ THEN 00656 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 1, 00657 $ X( 1, J ), INFO, WORK, IWORK ) 00658 ELSE 00659 RCOND_TMP = 0.0D+0 00660 END IF 00661 * 00662 * Cap the error at 1.0. 00663 * 00664 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 00665 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 00666 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00667 * 00668 * Threshold the error (see LAWN). 00669 * 00670 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 00671 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00672 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0 00673 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0 00674 $ .AND. INFO.LT.N + J ) INFO = N + J 00675 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) 00676 $ .LT. ERR_LBND ) THEN 00677 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND 00678 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00679 END IF 00680 * 00681 * Save the condition number. 00682 * 00683 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 00684 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP 00685 END IF 00686 00687 END DO 00688 END IF 00689 * 00690 RETURN 00691 * 00692 * End of DPORFSX 00693 * 00694 END