![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGERFSX 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGERFSX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgerfsx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgerfsx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgerfsx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGERFSX( 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, RWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * CHARACTER TRANS, 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 IPIV( * ) 00034 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00035 * $ X( LDX , * ), WORK( * ) 00036 * DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ), 00037 * $ ERR_BNDS_NORM( NRHS, * ), 00038 * $ ERR_BNDS_COMP( NRHS, * ), RWORK( * ) 00039 * .. 00040 * 00041 * 00042 *> \par Purpose: 00043 * ============= 00044 *> 00045 *> \verbatim 00046 *> 00047 *> ZGERFSX 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDAF,N) 00124 *> The factors L and U from the factorization A = P*L*U 00125 *> as computed by ZGETRF. 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 ZGETRF; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDX,NRHS) 00186 *> On entry, the solution matrix X, as computed by ZGETRS. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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) DOUBLE PRECISION 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.0D+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 COMPLEX*16 array, dimension (2*N) 00365 *> \endverbatim 00366 *> 00367 *> \param[out] RWORK 00368 *> \verbatim 00369 *> RWORK is DOUBLE PRECISION array, dimension (2*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 complex16GEcomputational 00410 * 00411 * ===================================================================== 00412 SUBROUTINE ZGERFSX( 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, RWORK, 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 DOUBLE PRECISION RCOND 00427 * .. 00428 * .. Array Arguments .. 00429 INTEGER IPIV( * ) 00430 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00431 $ X( LDX , * ), WORK( * ) 00432 DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ), 00433 $ ERR_BNDS_NORM( NRHS, * ), 00434 $ ERR_BNDS_COMP( NRHS, * ), RWORK( * ) 00435 * .. 00436 * 00437 * ================================================================== 00438 * 00439 * .. Parameters .. 00440 DOUBLE PRECISION ZERO, ONE 00441 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00442 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT 00443 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT 00444 DOUBLE PRECISION DZTHRESH_DEFAULT 00445 PARAMETER ( ITREF_DEFAULT = 1.0D+0 ) 00446 PARAMETER ( ITHRESH_DEFAULT = 10.0D+0 ) 00447 PARAMETER ( COMPONENTWISE_DEFAULT = 1.0D+0 ) 00448 PARAMETER ( RTHRESH_DEFAULT = 0.5D+0 ) 00449 PARAMETER ( DZTHRESH_DEFAULT = 0.25D+0 ) 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 DOUBLE PRECISION ANORM, RCOND_TMP 00466 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG 00467 LOGICAL IGNORE_CWISE 00468 INTEGER ITHRESH 00469 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH 00470 * .. 00471 * .. External Subroutines .. 00472 EXTERNAL XERBLA, ZGECON, ZLA_GERFSX_EXTENDED 00473 * .. 00474 * .. Intrinsic Functions .. 00475 INTRINSIC MAX, SQRT, TRANSFER 00476 * .. 00477 * .. External Functions .. 00478 EXTERNAL LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC 00479 EXTERNAL DLAMCH, ZLANGE, ZLA_GERCOND_X, ZLA_GERCOND_C 00480 DOUBLE PRECISION DLAMCH, ZLANGE, ZLA_GERCOND_X, ZLA_GERCOND_C 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.0D+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 = DBLE( N ) * DLAMCH( 'Epsilon' ) 00503 ITHRESH = INT( ITHRESH_DEFAULT ) 00504 RTHRESH = RTHRESH_DEFAULT 00505 UNSTABLE_THRESH = DZTHRESH_DEFAULT 00506 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0 00507 * 00508 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 00509 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+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.0D+0 ) THEN 00517 IF ( IGNORE_CWISE ) THEN 00518 PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0 00519 ELSE 00520 PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0 00521 END IF 00522 ELSE 00523 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+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( 'ZGERFSX', -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.0D+0 00567 DO J = 1, NRHS 00568 BERR( J ) = 0.0D+0 00569 IF ( N_ERR_BNDS .GE. 1 ) THEN 00570 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00571 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00572 END IF 00573 IF ( N_ERR_BNDS .GE. 2 ) THEN 00574 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0 00575 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0 00576 END IF 00577 IF ( N_ERR_BNDS .GE. 3 ) THEN 00578 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0 00579 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0 00580 END IF 00581 END DO 00582 RETURN 00583 END IF 00584 * 00585 * Default to failure. 00586 * 00587 RCOND = 0.0D+0 00588 DO J = 1, NRHS 00589 BERR( J ) = 1.0D+0 00590 IF ( N_ERR_BNDS .GE. 1 ) THEN 00591 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00592 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00593 END IF 00594 IF ( N_ERR_BNDS .GE. 2 ) THEN 00595 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00596 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00597 END IF 00598 IF ( N_ERR_BNDS .GE. 3 ) THEN 00599 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0 00600 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+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 = ZLANGE( NORM, N, N, A, LDA, RWORK ) 00613 CALL ZGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 00614 * 00615 * Perform refinement on each right-hand side 00616 * 00617 IF ( REF_TYPE .NE. 0 ) THEN 00618 00619 PREC_TYPE = ILAPREC( 'E' ) 00620 00621 IF ( NOTRAN ) THEN 00622 CALL ZLA_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, RWORK, WORK(N+1), 00626 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), 00627 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 00628 $ INFO ) 00629 ELSE 00630 CALL ZLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, 00631 $ NRHS, A, LDA, AF, LDAF, IPIV, ROWEQU, R, B, 00632 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, 00633 $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1), 00634 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), 00635 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 00636 $ INFO ) 00637 END IF 00638 END IF 00639 00640 ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' ) 00641 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN 00642 * 00643 * Compute scaled normwise condition number cond(A*C). 00644 * 00645 IF ( COLEQU .AND. NOTRAN ) THEN 00646 RCOND_TMP = ZLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV, 00647 $ C, .TRUE., INFO, WORK, RWORK ) 00648 ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN 00649 RCOND_TMP = ZLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV, 00650 $ R, .TRUE., INFO, WORK, RWORK ) 00651 ELSE 00652 RCOND_TMP = ZLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV, 00653 $ C, .FALSE., INFO, WORK, RWORK ) 00654 END IF 00655 DO J = 1, NRHS 00656 * 00657 * Cap the error at 1.0. 00658 * 00659 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 00660 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 00661 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00662 * 00663 * Threshold the error (see LAWN). 00664 * 00665 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 00666 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 00667 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0 00668 IF ( INFO .LE. N ) INFO = N + J 00669 ELSE IF (ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND) 00670 $ THEN 00671 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 00672 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 00673 END IF 00674 * 00675 * Save the condition number. 00676 * 00677 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 00678 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 00679 END IF 00680 END DO 00681 END IF 00682 00683 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2 ) THEN 00684 * 00685 * Compute componentwise condition number cond(A*diag(Y(:,J))) for 00686 * each right-hand side using the current solution as an estimate of 00687 * the true solution. If the componentwise error estimate is too 00688 * large, then the solution is a lousy estimate of truth and the 00689 * estimated RCOND may be too optimistic. To avoid misleading users, 00690 * the inverse condition number is set to 0.0 when the estimated 00691 * cwise error is at least CWISE_WRONG. 00692 * 00693 CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) ) 00694 DO J = 1, NRHS 00695 IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 00696 $ THEN 00697 RCOND_TMP = ZLA_GERCOND_X( TRANS, N, A, LDA, AF, LDAF, 00698 $ IPIV, X(1,J), INFO, WORK, RWORK ) 00699 ELSE 00700 RCOND_TMP = 0.0D+0 00701 END IF 00702 * 00703 * Cap the error at 1.0. 00704 * 00705 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 00706 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 00707 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00708 * 00709 * Threshold the error (see LAWN). 00710 * 00711 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 00712 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 00713 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0 00714 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0 00715 $ .AND. INFO.LT.N + J ) INFO = N + J 00716 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) 00717 $ .LT. ERR_LBND ) THEN 00718 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND 00719 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 00720 END IF 00721 * 00722 * Save the condition number. 00723 * 00724 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 00725 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP 00726 END IF 00727 00728 END DO 00729 END IF 00730 * 00731 RETURN 00732 * 00733 * End of ZGERFSX 00734 * 00735 END