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