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