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