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