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