00001 /* dla_gerfsx_extended.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 static doublereal c_b6 = -1.; 00020 static doublereal c_b8 = 1.; 00021 00022 /* Subroutine */ int dla_gerfsx_extended__(integer *prec_type__, integer * 00023 trans_type__, integer *n, integer *nrhs, doublereal *a, integer *lda, 00024 doublereal *af, integer *ldaf, integer *ipiv, logical *colequ, 00025 doublereal *c__, doublereal *b, integer *ldb, doublereal *y, integer * 00026 ldy, doublereal *berr_out__, integer *n_norms__, doublereal *errs_n__, 00027 doublereal *errs_c__, doublereal *res, doublereal *ayb, doublereal * 00028 dy, doublereal *y_tail__, doublereal *rcond, integer *ithresh, 00029 doublereal *rthresh, doublereal *dz_ub__, logical *ignore_cwise__, 00030 integer *info) 00031 { 00032 /* System generated locals */ 00033 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, y_dim1, 00034 y_offset, errs_n_dim1, errs_n_offset, errs_c_dim1, errs_c_offset, 00035 i__1, i__2, i__3; 00036 doublereal d__1, d__2; 00037 char ch__1[1]; 00038 00039 /* Local variables */ 00040 doublereal dxratmax, dzratmax; 00041 integer i__, j; 00042 extern /* Subroutine */ int dla_geamv__(integer *, integer *, integer *, 00043 doublereal *, doublereal *, integer *, doublereal *, integer *, 00044 doublereal *, doublereal *, integer *); 00045 logical incr_prec__; 00046 doublereal prev_dz_z__, yk, final_dx_x__; 00047 extern /* Subroutine */ int dla_wwaddw__(integer *, doublereal *, 00048 doublereal *, doublereal *); 00049 doublereal final_dz_z__, prevnormdx; 00050 integer cnt; 00051 doublereal dyk, eps, incr_thresh__, dx_x__, dz_z__; 00052 extern /* Subroutine */ int dla_lin_berr__(integer *, integer *, integer * 00053 , doublereal *, doublereal *, doublereal *); 00054 doublereal ymin; 00055 extern /* Subroutine */ int blas_dgemv_x__(integer *, integer *, integer * 00056 , doublereal *, doublereal *, integer *, doublereal *, integer *, 00057 doublereal *, doublereal *, integer *, integer *); 00058 integer y_prec_state__; 00059 extern /* Subroutine */ int blas_dgemv2_x__(integer *, integer *, integer 00060 *, doublereal *, doublereal *, integer *, doublereal *, 00061 doublereal *, integer *, doublereal *, doublereal *, integer *, 00062 integer *), dgemv_(char *, integer *, integer *, doublereal *, 00063 doublereal *, integer *, doublereal *, integer *, doublereal *, 00064 doublereal *, integer *), dcopy_(integer *, doublereal *, 00065 integer *, doublereal *, integer *); 00066 doublereal dxrat, dzrat; 00067 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 00068 integer *, doublereal *, integer *); 00069 char trans[1]; 00070 doublereal normx, normy; 00071 extern doublereal dlamch_(char *); 00072 extern /* Subroutine */ int dgetrs_(char *, integer *, integer *, 00073 doublereal *, integer *, integer *, doublereal *, integer *, 00074 integer *); 00075 doublereal normdx; 00076 extern /* Character */ VOID chla_transtype__(char *, ftnlen, integer *); 00077 doublereal hugeval; 00078 integer x_state__, z_state__; 00079 00080 00081 /* -- LAPACK routine (version 3.2.1) -- */ 00082 /* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */ 00083 /* -- Jason Riedy of Univ. of California Berkeley. -- */ 00084 /* -- April 2009 -- */ 00085 00086 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00087 /* -- Univ. of California Berkeley and NAG Ltd. -- */ 00088 00089 /* .. */ 00090 /* .. Scalar Arguments .. */ 00091 /* .. */ 00092 /* .. Array Arguments .. */ 00093 /* .. */ 00094 00095 /* Purpose */ 00096 /* ======= */ 00097 00098 /* DLA_GERFSX_EXTENDED improves the computed solution to a system of */ 00099 /* linear equations by performing extra-precise iterative refinement */ 00100 /* and provides error bounds and backward error estimates for the solution. */ 00101 /* This subroutine is called by DGERFSX to perform iterative refinement. */ 00102 /* In addition to normwise error bound, the code provides maximum */ 00103 /* componentwise error bound if possible. See comments for ERR_BNDS_NORM */ 00104 /* and ERR_BNDS_COMP for details of the error bounds. Note that this */ 00105 /* subroutine is only resonsible for setting the second fields of */ 00106 /* ERR_BNDS_NORM and ERR_BNDS_COMP. */ 00107 00108 /* Arguments */ 00109 /* ========= */ 00110 00111 /* PREC_TYPE (input) INTEGER */ 00112 /* Specifies the intermediate precision to be used in refinement. */ 00113 /* The value is defined by ILAPREC(P) where P is a CHARACTER and */ 00114 /* P = 'S': Single */ 00115 /* = 'D': Double */ 00116 /* = 'I': Indigenous */ 00117 /* = 'X', 'E': Extra */ 00118 00119 /* TRANS_TYPE (input) INTEGER */ 00120 /* Specifies the transposition operation on A. */ 00121 /* The value is defined by ILATRANS(T) where T is a CHARACTER and */ 00122 /* T = 'N': No transpose */ 00123 /* = 'T': Transpose */ 00124 /* = 'C': Conjugate transpose */ 00125 00126 /* N (input) INTEGER */ 00127 /* The number of linear equations, i.e., the order of the */ 00128 /* matrix A. N >= 0. */ 00129 00130 /* NRHS (input) INTEGER */ 00131 /* The number of right-hand-sides, i.e., the number of columns of the */ 00132 /* matrix B. */ 00133 00134 /* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ 00135 /* On entry, the N-by-N matrix A. */ 00136 00137 /* LDA (input) INTEGER */ 00138 /* The leading dimension of the array A. LDA >= max(1,N). */ 00139 00140 /* AF (input) DOUBLE PRECISION array, dimension (LDAF,N) */ 00141 /* The factors L and U from the factorization */ 00142 /* A = P*L*U as computed by DGETRF. */ 00143 00144 /* LDAF (input) INTEGER */ 00145 /* The leading dimension of the array AF. LDAF >= max(1,N). */ 00146 00147 /* IPIV (input) INTEGER array, dimension (N) */ 00148 /* The pivot indices from the factorization A = P*L*U */ 00149 /* as computed by DGETRF; row i of the matrix was interchanged */ 00150 /* with row IPIV(i). */ 00151 00152 /* COLEQU (input) LOGICAL */ 00153 /* If .TRUE. then column equilibration was done to A before calling */ 00154 /* this routine. This is needed to compute the solution and error */ 00155 /* bounds correctly. */ 00156 00157 /* C (input) DOUBLE PRECISION array, dimension (N) */ 00158 /* The column scale factors for A. If COLEQU = .FALSE., C */ 00159 /* is not accessed. If C is input, each element of C should be a power */ 00160 /* of the radix to ensure a reliable solution and error estimates. */ 00161 /* Scaling by powers of the radix does not cause rounding errors unless */ 00162 /* the result underflows or overflows. Rounding errors during scaling */ 00163 /* lead to refining with a matrix that is not equivalent to the */ 00164 /* input matrix, producing error estimates that may not be */ 00165 /* reliable. */ 00166 00167 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00168 /* The right-hand-side matrix B. */ 00169 00170 /* LDB (input) INTEGER */ 00171 /* The leading dimension of the array B. LDB >= max(1,N). */ 00172 00173 /* Y (input/output) DOUBLE PRECISION array, dimension */ 00174 /* (LDY,NRHS) */ 00175 /* On entry, the solution matrix X, as computed by DGETRS. */ 00176 /* On exit, the improved solution matrix Y. */ 00177 00178 /* LDY (input) INTEGER */ 00179 /* The leading dimension of the array Y. LDY >= max(1,N). */ 00180 00181 /* BERR_OUT (output) DOUBLE PRECISION array, dimension (NRHS) */ 00182 /* On exit, BERR_OUT(j) contains the componentwise relative backward */ 00183 /* error for right-hand-side j from the formula */ 00184 /* max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) */ 00185 /* where abs(Z) is the componentwise absolute value of the matrix */ 00186 /* or vector Z. This is computed by DLA_LIN_BERR. */ 00187 00188 /* N_NORMS (input) INTEGER */ 00189 /* Determines which error bounds to return (see ERR_BNDS_NORM */ 00190 /* and ERR_BNDS_COMP). */ 00191 /* If N_NORMS >= 1 return normwise error bounds. */ 00192 /* If N_NORMS >= 2 return componentwise error bounds. */ 00193 00194 /* ERR_BNDS_NORM (input/output) DOUBLE PRECISION array, dimension */ 00195 /* (NRHS, N_ERR_BNDS) */ 00196 /* For each right-hand side, this array contains information about */ 00197 /* various error bounds and condition numbers corresponding to the */ 00198 /* normwise relative error, which is defined as follows: */ 00199 00200 /* Normwise relative error in the ith solution vector: */ 00201 /* max_j (abs(XTRUE(j,i) - X(j,i))) */ 00202 /* ------------------------------ */ 00203 /* max_j abs(X(j,i)) */ 00204 00205 /* The array is indexed by the type of error information as described */ 00206 /* below. There currently are up to three pieces of information */ 00207 /* returned. */ 00208 00209 /* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith */ 00210 /* right-hand side. */ 00211 00212 /* The second index in ERR_BNDS_NORM(:,err) contains the following */ 00213 /* three fields: */ 00214 /* err = 1 "Trust/don't trust" boolean. Trust the answer if the */ 00215 /* reciprocal condition number is less than the threshold */ 00216 /* sqrt(n) * slamch('Epsilon'). */ 00217 00218 /* err = 2 "Guaranteed" error bound: The estimated forward error, */ 00219 /* almost certainly within a factor of 10 of the true error */ 00220 /* so long as the next entry is greater than the threshold */ 00221 /* sqrt(n) * slamch('Epsilon'). This error bound should only */ 00222 /* be trusted if the previous boolean is true. */ 00223 00224 /* err = 3 Reciprocal condition number: Estimated normwise */ 00225 /* reciprocal condition number. Compared with the threshold */ 00226 /* sqrt(n) * slamch('Epsilon') to determine if the error */ 00227 /* estimate is "guaranteed". These reciprocal condition */ 00228 /* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some */ 00229 /* appropriately scaled matrix Z. */ 00230 /* Let Z = S*A, where S scales each row by a power of the */ 00231 /* radix so all absolute row sums of Z are approximately 1. */ 00232 00233 /* This subroutine is only responsible for setting the second field */ 00234 /* above. */ 00235 /* See Lapack Working Note 165 for further details and extra */ 00236 /* cautions. */ 00237 00238 /* ERR_BNDS_COMP (input/output) DOUBLE PRECISION array, dimension */ 00239 /* (NRHS, N_ERR_BNDS) */ 00240 /* For each right-hand side, this array contains information about */ 00241 /* various error bounds and condition numbers corresponding to the */ 00242 /* componentwise relative error, which is defined as follows: */ 00243 00244 /* Componentwise relative error in the ith solution vector: */ 00245 /* abs(XTRUE(j,i) - X(j,i)) */ 00246 /* max_j ---------------------- */ 00247 /* abs(X(j,i)) */ 00248 00249 /* The array is indexed by the right-hand side i (on which the */ 00250 /* componentwise relative error depends), and the type of error */ 00251 /* information as described below. There currently are up to three */ 00252 /* pieces of information returned for each right-hand side. If */ 00253 /* componentwise accuracy is not requested (PARAMS(3) = 0.0), then */ 00254 /* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most */ 00255 /* the first (:,N_ERR_BNDS) entries are returned. */ 00256 00257 /* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith */ 00258 /* right-hand side. */ 00259 00260 /* The second index in ERR_BNDS_COMP(:,err) contains the following */ 00261 /* three fields: */ 00262 /* err = 1 "Trust/don't trust" boolean. Trust the answer if the */ 00263 /* reciprocal condition number is less than the threshold */ 00264 /* sqrt(n) * slamch('Epsilon'). */ 00265 00266 /* err = 2 "Guaranteed" error bound: The estimated forward error, */ 00267 /* almost certainly within a factor of 10 of the true error */ 00268 /* so long as the next entry is greater than the threshold */ 00269 /* sqrt(n) * slamch('Epsilon'). This error bound should only */ 00270 /* be trusted if the previous boolean is true. */ 00271 00272 /* err = 3 Reciprocal condition number: Estimated componentwise */ 00273 /* reciprocal condition number. Compared with the threshold */ 00274 /* sqrt(n) * slamch('Epsilon') to determine if the error */ 00275 /* estimate is "guaranteed". These reciprocal condition */ 00276 /* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some */ 00277 /* appropriately scaled matrix Z. */ 00278 /* Let Z = S*(A*diag(x)), where x is the solution for the */ 00279 /* current right-hand side and S scales each row of */ 00280 /* A*diag(x) by a power of the radix so all absolute row */ 00281 /* sums of Z are approximately 1. */ 00282 00283 /* This subroutine is only responsible for setting the second field */ 00284 /* above. */ 00285 /* See Lapack Working Note 165 for further details and extra */ 00286 /* cautions. */ 00287 00288 /* RES (input) DOUBLE PRECISION array, dimension (N) */ 00289 /* Workspace to hold the intermediate residual. */ 00290 00291 /* AYB (input) DOUBLE PRECISION array, dimension (N) */ 00292 /* Workspace. This can be the same workspace passed for Y_TAIL. */ 00293 00294 /* DY (input) DOUBLE PRECISION array, dimension (N) */ 00295 /* Workspace to hold the intermediate solution. */ 00296 00297 /* Y_TAIL (input) DOUBLE PRECISION array, dimension (N) */ 00298 /* Workspace to hold the trailing bits of the intermediate solution. */ 00299 00300 /* RCOND (input) DOUBLE PRECISION */ 00301 /* Reciprocal scaled condition number. This is an estimate of the */ 00302 /* reciprocal Skeel condition number of the matrix A after */ 00303 /* equilibration (if done). If this is less than the machine */ 00304 /* precision (in particular, if it is zero), the matrix is singular */ 00305 /* to working precision. Note that the error may still be small even */ 00306 /* if this number is very small and the matrix appears ill- */ 00307 /* conditioned. */ 00308 00309 /* ITHRESH (input) INTEGER */ 00310 /* The maximum number of residual computations allowed for */ 00311 /* refinement. The default is 10. For 'aggressive' set to 100 to */ 00312 /* permit convergence using approximate factorizations or */ 00313 /* factorizations other than LU. If the factorization uses a */ 00314 /* technique other than Gaussian elimination, the guarantees in */ 00315 /* ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy. */ 00316 00317 /* RTHRESH (input) DOUBLE PRECISION */ 00318 /* Determines when to stop refinement if the error estimate stops */ 00319 /* decreasing. Refinement will stop when the next solution no longer */ 00320 /* satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is */ 00321 /* the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The */ 00322 /* default value is 0.5. For 'aggressive' set to 0.9 to permit */ 00323 /* convergence on extremely ill-conditioned matrices. See LAWN 165 */ 00324 /* for more details. */ 00325 00326 /* DZ_UB (input) DOUBLE PRECISION */ 00327 /* Determines when to start considering componentwise convergence. */ 00328 /* Componentwise convergence is only considered after each component */ 00329 /* of the solution Y is stable, which we definte as the relative */ 00330 /* change in each component being less than DZ_UB. The default value */ 00331 /* is 0.25, requiring the first bit to be stable. See LAWN 165 for */ 00332 /* more details. */ 00333 00334 /* IGNORE_CWISE (input) LOGICAL */ 00335 /* If .TRUE. then ignore componentwise convergence. Default value */ 00336 /* is .FALSE.. */ 00337 00338 /* INFO (output) INTEGER */ 00339 /* = 0: Successful exit. */ 00340 /* < 0: if INFO = -i, the ith argument to DGETRS had an illegal */ 00341 /* value */ 00342 00343 /* ===================================================================== */ 00344 00345 /* .. Local Scalars .. */ 00346 /* .. */ 00347 /* .. Parameters .. */ 00348 /* .. */ 00349 /* .. External Subroutines .. */ 00350 /* .. */ 00351 /* .. Intrinsic Functions .. */ 00352 /* .. */ 00353 /* .. Executable Statements .. */ 00354 00355 /* Parameter adjustments */ 00356 errs_c_dim1 = *nrhs; 00357 errs_c_offset = 1 + errs_c_dim1; 00358 errs_c__ -= errs_c_offset; 00359 errs_n_dim1 = *nrhs; 00360 errs_n_offset = 1 + errs_n_dim1; 00361 errs_n__ -= errs_n_offset; 00362 a_dim1 = *lda; 00363 a_offset = 1 + a_dim1; 00364 a -= a_offset; 00365 af_dim1 = *ldaf; 00366 af_offset = 1 + af_dim1; 00367 af -= af_offset; 00368 --ipiv; 00369 --c__; 00370 b_dim1 = *ldb; 00371 b_offset = 1 + b_dim1; 00372 b -= b_offset; 00373 y_dim1 = *ldy; 00374 y_offset = 1 + y_dim1; 00375 y -= y_offset; 00376 --berr_out__; 00377 --res; 00378 --ayb; 00379 --dy; 00380 --y_tail__; 00381 00382 /* Function Body */ 00383 if (*info != 0) { 00384 return 0; 00385 } 00386 chla_transtype__(ch__1, (ftnlen)1, trans_type__); 00387 *(unsigned char *)trans = *(unsigned char *)&ch__1[0]; 00388 eps = dlamch_("Epsilon"); 00389 hugeval = dlamch_("Overflow"); 00390 /* Force HUGEVAL to Inf */ 00391 hugeval *= hugeval; 00392 /* Using HUGEVAL may lead to spurious underflows. */ 00393 incr_thresh__ = (doublereal) (*n) * eps; 00394 00395 i__1 = *nrhs; 00396 for (j = 1; j <= i__1; ++j) { 00397 y_prec_state__ = 1; 00398 if (y_prec_state__ == 2) { 00399 i__2 = *n; 00400 for (i__ = 1; i__ <= i__2; ++i__) { 00401 y_tail__[i__] = 0.; 00402 } 00403 } 00404 dxrat = 0.; 00405 dxratmax = 0.; 00406 dzrat = 0.; 00407 dzratmax = 0.; 00408 final_dx_x__ = hugeval; 00409 final_dz_z__ = hugeval; 00410 prevnormdx = hugeval; 00411 prev_dz_z__ = hugeval; 00412 dz_z__ = hugeval; 00413 dx_x__ = hugeval; 00414 x_state__ = 1; 00415 z_state__ = 0; 00416 incr_prec__ = FALSE_; 00417 i__2 = *ithresh; 00418 for (cnt = 1; cnt <= i__2; ++cnt) { 00419 00420 /* Compute residual RES = B_s - op(A_s) * Y, */ 00421 /* op(A) = A, A**T, or A**H depending on TRANS (and type). */ 00422 00423 dcopy_(n, &b[j * b_dim1 + 1], &c__1, &res[1], &c__1); 00424 if (y_prec_state__ == 0) { 00425 dgemv_(trans, n, n, &c_b6, &a[a_offset], lda, &y[j * y_dim1 + 00426 1], &c__1, &c_b8, &res[1], &c__1); 00427 } else if (y_prec_state__ == 1) { 00428 blas_dgemv_x__(trans_type__, n, n, &c_b6, &a[a_offset], lda, & 00429 y[j * y_dim1 + 1], &c__1, &c_b8, &res[1], &c__1, 00430 prec_type__); 00431 } else { 00432 blas_dgemv2_x__(trans_type__, n, n, &c_b6, &a[a_offset], lda, 00433 &y[j * y_dim1 + 1], &y_tail__[1], &c__1, &c_b8, &res[ 00434 1], &c__1, prec_type__); 00435 } 00436 /* XXX: RES is no longer needed. */ 00437 dcopy_(n, &res[1], &c__1, &dy[1], &c__1); 00438 dgetrs_(trans, n, &c__1, &af[af_offset], ldaf, &ipiv[1], &dy[1], 00439 n, info); 00440 00441 /* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. */ 00442 00443 normx = 0.; 00444 normy = 0.; 00445 normdx = 0.; 00446 dz_z__ = 0.; 00447 ymin = hugeval; 00448 00449 i__3 = *n; 00450 for (i__ = 1; i__ <= i__3; ++i__) { 00451 yk = (d__1 = y[i__ + j * y_dim1], abs(d__1)); 00452 dyk = (d__1 = dy[i__], abs(d__1)); 00453 if (yk != 0.) { 00454 /* Computing MAX */ 00455 d__1 = dz_z__, d__2 = dyk / yk; 00456 dz_z__ = max(d__1,d__2); 00457 } else if (dyk != 0.) { 00458 dz_z__ = hugeval; 00459 } 00460 ymin = min(ymin,yk); 00461 normy = max(normy,yk); 00462 if (*colequ) { 00463 /* Computing MAX */ 00464 d__1 = normx, d__2 = yk * c__[i__]; 00465 normx = max(d__1,d__2); 00466 /* Computing MAX */ 00467 d__1 = normdx, d__2 = dyk * c__[i__]; 00468 normdx = max(d__1,d__2); 00469 } else { 00470 normx = normy; 00471 normdx = max(normdx,dyk); 00472 } 00473 } 00474 if (normx != 0.) { 00475 dx_x__ = normdx / normx; 00476 } else if (normdx == 0.) { 00477 dx_x__ = 0.; 00478 } else { 00479 dx_x__ = hugeval; 00480 } 00481 dxrat = normdx / prevnormdx; 00482 dzrat = dz_z__ / prev_dz_z__; 00483 00484 /* Check termination criteria */ 00485 00486 if (! (*ignore_cwise__) && ymin * *rcond < incr_thresh__ * normy 00487 && y_prec_state__ < 2) { 00488 incr_prec__ = TRUE_; 00489 } 00490 if (x_state__ == 3 && dxrat <= *rthresh) { 00491 x_state__ = 1; 00492 } 00493 if (x_state__ == 1) { 00494 if (dx_x__ <= eps) { 00495 x_state__ = 2; 00496 } else if (dxrat > *rthresh) { 00497 if (y_prec_state__ != 2) { 00498 incr_prec__ = TRUE_; 00499 } else { 00500 x_state__ = 3; 00501 } 00502 } else { 00503 if (dxrat > dxratmax) { 00504 dxratmax = dxrat; 00505 } 00506 } 00507 if (x_state__ > 1) { 00508 final_dx_x__ = dx_x__; 00509 } 00510 } 00511 if (z_state__ == 0 && dz_z__ <= *dz_ub__) { 00512 z_state__ = 1; 00513 } 00514 if (z_state__ == 3 && dzrat <= *rthresh) { 00515 z_state__ = 1; 00516 } 00517 if (z_state__ == 1) { 00518 if (dz_z__ <= eps) { 00519 z_state__ = 2; 00520 } else if (dz_z__ > *dz_ub__) { 00521 z_state__ = 0; 00522 dzratmax = 0.; 00523 final_dz_z__ = hugeval; 00524 } else if (dzrat > *rthresh) { 00525 if (y_prec_state__ != 2) { 00526 incr_prec__ = TRUE_; 00527 } else { 00528 z_state__ = 3; 00529 } 00530 } else { 00531 if (dzrat > dzratmax) { 00532 dzratmax = dzrat; 00533 } 00534 } 00535 if (z_state__ > 1) { 00536 final_dz_z__ = dz_z__; 00537 } 00538 } 00539 00540 /* Exit if both normwise and componentwise stopped working, */ 00541 /* but if componentwise is unstable, let it go at least two */ 00542 /* iterations. */ 00543 00544 if (x_state__ != 1) { 00545 if (*ignore_cwise__) { 00546 goto L666; 00547 } 00548 if (z_state__ == 3 || z_state__ == 2) { 00549 goto L666; 00550 } 00551 if (z_state__ == 0 && cnt > 1) { 00552 goto L666; 00553 } 00554 } 00555 if (incr_prec__) { 00556 incr_prec__ = FALSE_; 00557 ++y_prec_state__; 00558 i__3 = *n; 00559 for (i__ = 1; i__ <= i__3; ++i__) { 00560 y_tail__[i__] = 0.; 00561 } 00562 } 00563 prevnormdx = normdx; 00564 prev_dz_z__ = dz_z__; 00565 00566 /* Update soluton. */ 00567 00568 if (y_prec_state__ < 2) { 00569 daxpy_(n, &c_b8, &dy[1], &c__1, &y[j * y_dim1 + 1], &c__1); 00570 } else { 00571 dla_wwaddw__(n, &y[j * y_dim1 + 1], &y_tail__[1], &dy[1]); 00572 } 00573 } 00574 /* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. */ 00575 L666: 00576 00577 /* Set final_* when cnt hits ithresh. */ 00578 00579 if (x_state__ == 1) { 00580 final_dx_x__ = dx_x__; 00581 } 00582 if (z_state__ == 1) { 00583 final_dz_z__ = dz_z__; 00584 } 00585 00586 /* Compute error bounds */ 00587 00588 if (*n_norms__ >= 1) { 00589 errs_n__[j + (errs_n_dim1 << 1)] = final_dx_x__ / (1 - dxratmax); 00590 } 00591 if (*n_norms__ >= 2) { 00592 errs_c__[j + (errs_c_dim1 << 1)] = final_dz_z__ / (1 - dzratmax); 00593 } 00594 00595 /* Compute componentwise relative backward error from formula */ 00596 /* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) */ 00597 /* where abs(Z) is the componentwise absolute value of the matrix */ 00598 /* or vector Z. */ 00599 00600 /* Compute residual RES = B_s - op(A_s) * Y, */ 00601 /* op(A) = A, A**T, or A**H depending on TRANS (and type). */ 00602 00603 dcopy_(n, &b[j * b_dim1 + 1], &c__1, &res[1], &c__1); 00604 dgemv_(trans, n, n, &c_b6, &a[a_offset], lda, &y[j * y_dim1 + 1], & 00605 c__1, &c_b8, &res[1], &c__1); 00606 i__2 = *n; 00607 for (i__ = 1; i__ <= i__2; ++i__) { 00608 ayb[i__] = (d__1 = b[i__ + j * b_dim1], abs(d__1)); 00609 } 00610 00611 /* Compute abs(op(A_s))*abs(Y) + abs(B_s). */ 00612 00613 dla_geamv__(trans_type__, n, n, &c_b8, &a[a_offset], lda, &y[j * 00614 y_dim1 + 1], &c__1, &c_b8, &ayb[1], &c__1); 00615 dla_lin_berr__(n, n, &c__1, &res[1], &ayb[1], &berr_out__[j]); 00616 00617 /* End of loop for each RHS. */ 00618 00619 } 00620 00621 return 0; 00622 } /* dla_gerfsx_extended__ */