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