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