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