1 #ifdef FLA_ENABLE_XBLAS
2 /* ../netlib/sla_gerfsx_extended.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
3 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
4 #include "FLA_f2c.h" /* Table of constant values */
5 static integer c__1 = 1;
6 static real c_b6 = -1.f;
7 static real c_b8 = 1.f;
8 /* > \brief \b SLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. */
9 /* =========== DOCUMENTATION =========== */
10 /* Online html documentation available at */
11 /* http://www.netlib.org/lapack/explore-html/ */
12 /* > \htmlonly */
13 /* > Download SLA_GERFSX_EXTENDED + dependencies */
14 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_ger fsx_extended.f"> */
15 /* > [TGZ]</a> */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_ger fsx_extended.f"> */
17 /* > [ZIP]</a> */
18 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_ger fsx_extended.f"> */
19 /* > [TXT]</a> */
20 /* > \endhtmlonly */
21 /* Definition: */
22 /* =========== */
23 /* SUBROUTINE SLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS, A, */
24 /* LDA, AF, LDAF, IPIV, COLEQU, C, B, */
25 /* LDB, Y, LDY, BERR_OUT, N_NORMS, */
26 /* ERRS_N, ERRS_C, RES, */
27 /* AYB, DY, Y_TAIL, RCOND, ITHRESH, */
28 /* RTHRESH, DZ_UB, IGNORE_CWISE, */
29 /* INFO ) */
30 /* .. Scalar Arguments .. */
31 /* INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE, */
32 /* $ TRANS_TYPE, N_NORMS, ITHRESH */
33 /* LOGICAL COLEQU, IGNORE_CWISE */
34 /* REAL RTHRESH, DZ_UB */
35 /* .. */
36 /* .. Array Arguments .. */
37 /* INTEGER IPIV( * ) */
38 /* REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ), */
39 /* $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * ) */
40 /* REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ), */
41 /* $ ERRS_N( NRHS, * ), */
42 /* $ ERRS_C( NRHS, * ) */
43 /* .. */
44 /* > \par Purpose: */
45 /* ============= */
46 /* > */
47 /* > \verbatim */
48 /* > */
49 /* > SLA_GERFSX_EXTENDED improves the computed solution to a system of */
50 /* > linear equations by performing extra-precise iterative refinement */
51 /* > and provides error bounds and backward error estimates for the solution. */
52 /* > This subroutine is called by SGERFSX to perform iterative refinement. */
53 /* > In addition to normwise error bound, the code provides maximum */
54 /* > componentwise error bound if possible. See comments for ERRS_N */
55 /* > and ERRS_C for details of the error bounds. Note that this */
56 /* > subroutine is only resonsible for setting the second fields of */
57 /* > ERRS_N and ERRS_C. */
58 /* > \endverbatim */
59 /* Arguments: */
60 /* ========== */
61 /* > \param[in] PREC_TYPE */
62 /* > \verbatim */
63 /* > PREC_TYPE is INTEGER */
64 /* > Specifies the intermediate precision to be used in refinement. */
65 /* > The value is defined by ILAPREC(P) where P is a CHARACTER and */
66 /* > P = 'S': Single */
67 /* > = 'D': Double */
68 /* > = 'I': Indigenous */
69 /* > = 'X', 'E': Extra */
70 /* > \endverbatim */
71 /* > */
72 /* > \param[in] TRANS_TYPE */
73 /* > \verbatim */
74 /* > TRANS_TYPE is INTEGER */
75 /* > Specifies the transposition operation on A. */
76 /* > The value is defined by ILATRANS(T) where T is a CHARACTER and */
77 /* > T = 'N': No transpose */
78 /* > = 'T': Transpose */
79 /* > = 'C': Conjugate transpose */
80 /* > \endverbatim */
81 /* > */
82 /* > \param[in] N */
83 /* > \verbatim */
84 /* > N is INTEGER */
85 /* > The number of linear equations, i.e., the order of the */
86 /* > matrix A. N >= 0. */
87 /* > \endverbatim */
88 /* > */
89 /* > \param[in] NRHS */
90 /* > \verbatim */
91 /* > NRHS is INTEGER */
92 /* > The number of right-hand-sides, i.e., the number of columns of the */
93 /* > matrix B. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[in] A */
97 /* > \verbatim */
98 /* > A is REAL array, dimension (LDA,N) */
99 /* > On entry, the N-by-N matrix A. */
100 /* > \endverbatim */
101 /* > */
102 /* > \param[in] LDA */
103 /* > \verbatim */
104 /* > LDA is INTEGER */
105 /* > The leading dimension of the array A. LDA >= max(1,N). */
106 /* > \endverbatim */
107 /* > */
108 /* > \param[in] AF */
109 /* > \verbatim */
110 /* > AF is REAL array, dimension (LDAF,N) */
111 /* > The factors L and U from the factorization */
112 /* > A = P*L*U as computed by SGETRF. */
113 /* > \endverbatim */
114 /* > */
115 /* > \param[in] LDAF */
116 /* > \verbatim */
117 /* > LDAF is INTEGER */
118 /* > The leading dimension of the array AF. LDAF >= max(1,N). */
119 /* > \endverbatim */
120 /* > */
121 /* > \param[in] IPIV */
122 /* > \verbatim */
123 /* > IPIV is INTEGER array, dimension (N) */
124 /* > The pivot indices from the factorization A = P*L*U */
125 /* > as computed by SGETRF;
126 row i of the matrix was interchanged */
127 /* > with row IPIV(i). */
128 /* > \endverbatim */
129 /* > */
130 /* > \param[in] COLEQU */
131 /* > \verbatim */
132 /* > COLEQU is LOGICAL */
133 /* > If .TRUE. then column equilibration was done to A before calling */
134 /* > this routine. This is needed to compute the solution and error */
135 /* > bounds correctly. */
136 /* > \endverbatim */
137 /* > */
138 /* > \param[in] C */
139 /* > \verbatim */
140 /* > C is REAL array, dimension (N) */
141 /* > The column scale factors for A. If COLEQU = .FALSE., C */
142 /* > is not accessed. If C is input, each element of C should be a power */
143 /* > of the radix to ensure a reliable solution and error estimates. */
144 /* > Scaling by powers of the radix does not cause rounding errors unless */
145 /* > the result underflows or overflows. Rounding errors during scaling */
146 /* > lead to refining with a matrix that is not equivalent to the */
147 /* > input matrix, producing error estimates that may not be */
148 /* > reliable. */
149 /* > \endverbatim */
150 /* > */
151 /* > \param[in] B */
152 /* > \verbatim */
153 /* > B is REAL array, dimension (LDB,NRHS) */
154 /* > The right-hand-side matrix B. */
155 /* > \endverbatim */
156 /* > */
157 /* > \param[in] LDB */
158 /* > \verbatim */
159 /* > LDB is INTEGER */
160 /* > The leading dimension of the array B. LDB >= max(1,N). */
161 /* > \endverbatim */
162 /* > */
163 /* > \param[in,out] Y */
164 /* > \verbatim */
165 /* > Y is REAL array, dimension (LDY,NRHS) */
166 /* > On entry, the solution matrix X, as computed by SGETRS. */
167 /* > On exit, the improved solution matrix Y. */
168 /* > \endverbatim */
169 /* > */
170 /* > \param[in] LDY */
171 /* > \verbatim */
172 /* > LDY is INTEGER */
173 /* > The leading dimension of the array Y. LDY >= max(1,N). */
174 /* > \endverbatim */
175 /* > */
176 /* > \param[out] BERR_OUT */
177 /* > \verbatim */
178 /* > BERR_OUT is REAL array, dimension (NRHS) */
179 /* > On exit, BERR_OUT(j) contains the componentwise relative backward */
180 /* > error for right-hand-side j from the formula */
181 /* > max(i) ( f2c_abs(RES(i)) / ( f2c_abs(op(A_s))*f2c_abs(Y) + f2c_abs(B_s) )(i) ) */
182 /* > where f2c_abs(Z) is the componentwise absolute value of the matrix */
183 /* > or vector Z. This is computed by SLA_LIN_BERR. */
184 /* > \endverbatim */
185 /* > */
186 /* > \param[in] N_NORMS */
187 /* > \verbatim */
188 /* > N_NORMS is INTEGER */
189 /* > Determines which error bounds to return (see ERRS_N */
190 /* > and ERRS_C). */
191 /* > If N_NORMS >= 1 return normwise error bounds. */
192 /* > If N_NORMS >= 2 return componentwise error bounds. */
193 /* > \endverbatim */
194 /* > */
195 /* > \param[in,out] ERRS_N */
196 /* > \verbatim */
197 /* > ERRS_N is REAL array, dimension (NRHS, N_ERR_BNDS) */
198 /* > For each right-hand side, this array contains information about */
199 /* > various error bounds and condition numbers corresponding to the */
200 /* > normwise relative error, which is defined as follows: */
201 /* > */
202 /* > Normwise relative error in the ith solution vector: */
203 /* > max_j (f2c_abs(XTRUE(j,i) - X(j,i))) */
204 /* > ------------------------------ */
205 /* > max_j f2c_abs(X(j,i)) */
206 /* > */
207 /* > The array is indexed by the type of error information as described */
208 /* > below. There currently are up to three pieces of information */
209 /* > returned. */
210 /* > */
211 /* > The first index in ERRS_N(i,:) corresponds to the ith */
212 /* > right-hand side. */
213 /* > */
214 /* > The second index in ERRS_N(:,err) contains the following */
215 /* > three fields: */
216 /* > err = 1 "Trust/don't trust" boolean. Trust the answer if the */
217 /* > reciprocal condition number is less than the threshold */
218 /* > sqrt(n) * slamch('Epsilon'). */
219 /* > */
220 /* > err = 2 "Guaranteed" error bound: The estimated forward error, */
221 /* > almost certainly within a factor of 10 of the true error */
222 /* > so long as the next entry is greater than the threshold */
223 /* > sqrt(n) * slamch('Epsilon'). This error bound should only */
224 /* > be trusted if the previous boolean is true. */
225 /* > */
226 /* > err = 3 Reciprocal condition number: Estimated normwise */
227 /* > reciprocal condition number. Compared with the threshold */
228 /* > sqrt(n) * slamch('Epsilon') to determine if the error */
229 /* > estimate is "guaranteed". These reciprocal condition */
230 /* > numbers are 1 / (norm(Z^{
231 -1}
232 ,inf) * norm(Z,inf)) for some */
233 /* > appropriately scaled matrix Z. */
234 /* > Let Z = S*A, where S scales each row by a power of the */
235 /* > radix so all absolute row sums of Z are approximately 1. */
236 /* > */
237 /* > This subroutine is only responsible for setting the second field */
238 /* > above. */
239 /* > See Lapack Working Note 165 for further details and extra */
240 /* > cautions. */
241 /* > \endverbatim */
242 /* > */
243 /* > \param[in,out] ERRS_C */
244 /* > \verbatim */
245 /* > ERRS_C is REAL array, dimension (NRHS, N_ERR_BNDS) */
246 /* > For each right-hand side, this array contains information about */
247 /* > various error bounds and condition numbers corresponding to the */
248 /* > componentwise relative error, which is defined as follows: */
249 /* > */
250 /* > Componentwise relative error in the ith solution vector: */
251 /* > f2c_abs(XTRUE(j,i) - X(j,i)) */
252 /* > max_j ---------------------- */
253 /* > f2c_abs(X(j,i)) */
254 /* > */
255 /* > The array is indexed by the right-hand side i (on which the */
256 /* > componentwise relative error depends), and the type of error */
257 /* > information as described below. There currently are up to three */
258 /* > pieces of information returned for each right-hand side. If */
259 /* > componentwise accuracy is not requested (PARAMS(3) = 0.0), then */
260 /* > ERRS_C is not accessed. If N_ERR_BNDS .LT. 3, then at most */
261 /* > the first (:,N_ERR_BNDS) entries are returned. */
262 /* > */
263 /* > The first index in ERRS_C(i,:) corresponds to the ith */
264 /* > right-hand side. */
265 /* > */
266 /* > The second index in ERRS_C(:,err) contains the following */
267 /* > three fields: */
268 /* > err = 1 "Trust/don't trust" boolean. Trust the answer if the */
269 /* > reciprocal condition number is less than the threshold */
270 /* > sqrt(n) * slamch('Epsilon'). */
271 /* > */
272 /* > err = 2 "Guaranteed" error bound: The estimated forward error, */
273 /* > almost certainly within a factor of 10 of the true error */
274 /* > so long as the next entry is greater than the threshold */
275 /* > sqrt(n) * slamch('Epsilon'). This error bound should only */
276 /* > be trusted if the previous boolean is true. */
277 /* > */
278 /* > err = 3 Reciprocal condition number: Estimated componentwise */
279 /* > reciprocal condition number. Compared with the threshold */
280 /* > sqrt(n) * slamch('Epsilon') to determine if the error */
281 /* > estimate is "guaranteed". These reciprocal condition */
282 /* > numbers are 1 / (norm(Z^{
283 -1}
284 ,inf) * norm(Z,inf)) for some */
285 /* > appropriately scaled matrix Z. */
286 /* > Let Z = S*(A*diag(x)), where x is the solution for the */
287 /* > current right-hand side and S scales each row of */
288 /* > A*diag(x) by a power of the radix so all absolute row */
289 /* > sums of Z are approximately 1. */
290 /* > */
291 /* > This subroutine is only responsible for setting the second field */
292 /* > above. */
293 /* > See Lapack Working Note 165 for further details and extra */
294 /* > cautions. */
295 /* > \endverbatim */
296 /* > */
297 /* > \param[in] RES */
298 /* > \verbatim */
299 /* > RES is REAL array, dimension (N) */
300 /* > Workspace to hold the intermediate residual. */
301 /* > \endverbatim */
302 /* > */
303 /* > \param[in] AYB */
304 /* > \verbatim */
305 /* > AYB is REAL array, dimension (N) */
306 /* > Workspace. This can be the same workspace passed for Y_TAIL. */
307 /* > \endverbatim */
308 /* > */
309 /* > \param[in] DY */
310 /* > \verbatim */
311 /* > DY is REAL array, dimension (N) */
312 /* > Workspace to hold the intermediate solution. */
313 /* > \endverbatim */
314 /* > */
315 /* > \param[in] Y_TAIL */
316 /* > \verbatim */
317 /* > Y_TAIL is REAL array, dimension (N) */
318 /* > Workspace to hold the trailing bits of the intermediate solution. */
319 /* > \endverbatim */
320 /* > */
321 /* > \param[in] RCOND */
322 /* > \verbatim */
323 /* > RCOND is REAL */
324 /* > Reciprocal scaled condition number. This is an estimate of the */
325 /* > reciprocal Skeel condition number of the matrix A after */
326 /* > equilibration (if done). If this is less than the machine */
327 /* > precision (in particular, if it is zero), the matrix is singular */
328 /* > to working precision. Note that the error may still be small even */
329 /* > if this number is very small and the matrix appears ill- */
330 /* > conditioned. */
331 /* > \endverbatim */
332 /* > */
333 /* > \param[in] ITHRESH */
334 /* > \verbatim */
335 /* > ITHRESH is INTEGER */
336 /* > The maximum number of residual computations allowed for */
337 /* > refinement. The default is 10. For 'aggressive' set to 100 to */
338 /* > permit convergence using approximate factorizations or */
339 /* > factorizations other than LU. If the factorization uses a */
340 /* > technique other than Gaussian elimination, the guarantees in */
341 /* > ERRS_N and ERRS_C may no longer be trustworthy. */
342 /* > \endverbatim */
343 /* > */
344 /* > \param[in] RTHRESH */
345 /* > \verbatim */
346 /* > RTHRESH is REAL */
347 /* > Determines when to stop refinement if the error estimate stops */
348 /* > decreasing. Refinement will stop when the next solution no longer */
349 /* > satisfies norm(dx_{
350 i+1}
351 ) < RTHRESH * norm(dx_i) where norm(Z) is */
352 /* > the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The */
353 /* > default value is 0.5. For 'aggressive' set to 0.9 to permit */
354 /* > convergence on extremely ill-conditioned matrices. See LAWN 165 */
355 /* > for more details. */
356 /* > \endverbatim */
357 /* > */
358 /* > \param[in] DZ_UB */
359 /* > \verbatim */
360 /* > DZ_UB is REAL */
361 /* > Determines when to start considering componentwise convergence. */
362 /* > Componentwise convergence is only considered after each component */
363 /* > of the solution Y is stable, which we definte as the relative */
364 /* > change in each component being less than DZ_UB. The default value */
365 /* > is 0.25, requiring the first bit to be stable. See LAWN 165 for */
366 /* > more details. */
367 /* > \endverbatim */
368 /* > */
369 /* > \param[in] IGNORE_CWISE */
370 /* > \verbatim */
371 /* > IGNORE_CWISE is LOGICAL */
372 /* > If .TRUE. then ignore componentwise convergence. Default value */
373 /* > is .FALSE.. */
374 /* > \endverbatim */
375 /* > */
376 /* > \param[out] INFO */
377 /* > \verbatim */
378 /* > INFO is INTEGER */
379 /* > = 0: Successful exit. */
380 /* > < 0: if INFO = -i, the ith argument to SGETRS had an illegal */
381 /* > value */
382 /* > \endverbatim */
383 /* Authors: */
384 /* ======== */
385 /* > \author Univ. of Tennessee */
386 /* > \author Univ. of California Berkeley */
387 /* > \author Univ. of Colorado Denver */
388 /* > \author NAG Ltd. */
389 /* > \date September 2012 */
390 /* > \ingroup realGEcomputational */
391 /* ===================================================================== */
392 /* Subroutine */
sla_gerfsx_extended_(integer * prec_type__,integer * trans_type__,integer * n,integer * nrhs,real * a,integer * lda,real * af,integer * ldaf,integer * ipiv,logical * colequ,real * c__,real * b,integer * ldb,real * y,integer * ldy,real * berr_out__,integer * n_norms__,real * errs_n__,real * errs_c__,real * res,real * ayb,real * dy,real * y_tail__,real * rcond,integer * ithresh,real * rthresh,real * dz_ub__,logical * ignore_cwise__,integer * info)393 int sla_gerfsx_extended_(integer *prec_type__, integer * trans_type__, integer *n, integer *nrhs, real *a, integer *lda, real * af, integer *ldaf, integer *ipiv, logical *colequ, real *c__, real *b, integer *ldb, real *y, integer *ldy, real *berr_out__, integer * n_norms__, real *errs_n__, real *errs_c__, real *res, real *ayb, real *dy, real *y_tail__, real *rcond, integer *ithresh, real *rthresh, real *dz_ub__, logical *ignore_cwise__, integer *info)
394 {
395 /* System generated locals */
396 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, y_dim1, y_offset, errs_n_dim1, errs_n_offset, errs_c_dim1, errs_c_offset, i__1, i__2, i__3;
397 real r__1, r__2;
398 char ch__1[1];
399 /* Local variables */
400 real dxratmax, dzratmax;
401 integer i__, j;
402 extern /* Subroutine */
403 int sla_geamv_(integer *, integer *, integer *, real *, real *, integer *, real *, integer *, real *, real *, integer *);
404 logical incr_prec__;
405 real prev_dz_z__, yk, final_dx_x__, final_dz_z__;
406 extern /* Subroutine */
407 int sla_wwaddw_(integer *, real *, real *, real * );
408 real prevnormdx;
409 integer cnt;
410 real dyk, eps, incr_thresh__, dx_x__, dz_z__, ymin;
411 extern /* Subroutine */
412 int sla_lin_berr_(integer *, integer *, integer * , real *, real *, real *), blas_sgemv_x_(integer *, integer *, integer *, real *, real *, integer *, real *, integer *, real *, real *, integer *, integer *);
413 integer y_prec_state__;
414 extern /* Subroutine */
415 int blas_sgemv2_x_(integer *, integer *, integer *, real *, real *, integer *, real *, real *, integer *, real *, real *, integer *, integer *), sgemv_(char *, integer *, integer * , real *, real *, integer *, real *, integer *, real *, real *, integer *);
416 real dxrat, dzrat;
417 char trans[1];
418 extern /* Subroutine */
419 int scopy_(integer *, real *, integer *, real *, integer *);
420 real normx, normy;
421 extern /* Subroutine */
422 int saxpy_(integer *, real *, real *, integer *, real *, integer *);
423 extern real slamch_(char *);
424 real normdx;
425 extern /* Subroutine */
426 int sgetrs_(char *, integer *, integer *, real *, integer *, integer *, real *, integer *, integer *);
427 extern /* Character */
428 VOID chla_transtype_(char *, integer *);
429 real hugeval;
430 integer x_state__, z_state__;
431 /* -- LAPACK computational routine (version 3.4.2) -- */
432 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
433 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
434 /* September 2012 */
435 /* .. Scalar Arguments .. */
436 /* .. */
437 /* .. Array Arguments .. */
438 /* .. */
439 /* ===================================================================== */
440 /* .. Local Scalars .. */
441 /* .. */
442 /* .. Parameters .. */
443 /* .. */
444 /* .. External Subroutines .. */
445 /* .. */
446 /* .. Intrinsic Functions .. */
447 /* .. */
448 /* .. Executable Statements .. */
449 /* Parameter adjustments */
450 errs_c_dim1 = *nrhs;
451 errs_c_offset = 1 + errs_c_dim1;
452 errs_c__ -= errs_c_offset;
453 errs_n_dim1 = *nrhs;
454 errs_n_offset = 1 + errs_n_dim1;
455 errs_n__ -= errs_n_offset;
456 a_dim1 = *lda;
457 a_offset = 1 + a_dim1;
458 a -= a_offset;
459 af_dim1 = *ldaf;
460 af_offset = 1 + af_dim1;
461 af -= af_offset;
462 --ipiv;
463 --c__;
464 b_dim1 = *ldb;
465 b_offset = 1 + b_dim1;
466 b -= b_offset;
467 y_dim1 = *ldy;
468 y_offset = 1 + y_dim1;
469 y -= y_offset;
470 --berr_out__;
471 --res;
472 --ayb;
473 --dy;
474 --y_tail__;
475 /* Function Body */
476 if (*info != 0)
477 {
478 return 0;
479 }
480 chla_transtype_(ch__1, trans_type__);
481 *(unsigned char *)trans = *(unsigned char *)&ch__1[0];
482 eps = slamch_("Epsilon");
483 hugeval = slamch_("Overflow");
484 /* Force HUGEVAL to Inf */
485 hugeval *= hugeval;
486 /* Using HUGEVAL may lead to spurious underflows. */
487 incr_thresh__ = (real) (*n) * eps;
488 i__1 = *nrhs;
489 for (j = 1;
490 j <= i__1;
491 ++j)
492 {
493 y_prec_state__ = 1;
494 if (y_prec_state__ == 2)
495 {
496 i__2 = *n;
497 for (i__ = 1;
498 i__ <= i__2;
499 ++i__)
500 {
501 y_tail__[i__] = 0.f;
502 }
503 }
504 dxrat = 0.f;
505 dxratmax = 0.f;
506 dzrat = 0.f;
507 dzratmax = 0.f;
508 final_dx_x__ = hugeval;
509 final_dz_z__ = hugeval;
510 prevnormdx = hugeval;
511 prev_dz_z__ = hugeval;
512 dz_z__ = hugeval;
513 dx_x__ = hugeval;
514 x_state__ = 1;
515 z_state__ = 0;
516 incr_prec__ = FALSE_;
517 i__2 = *ithresh;
518 for (cnt = 1;
519 cnt <= i__2;
520 ++cnt)
521 {
522 /* Compute residual RES = B_s - op(A_s) * Y, */
523 /* op(A) = A, A**T, or A**H depending on TRANS (and type). */
524 scopy_(n, &b[j * b_dim1 + 1], &c__1, &res[1], &c__1);
525 if (y_prec_state__ == 0)
526 {
527 sgemv_(trans, n, n, &c_b6, &a[a_offset], lda, &y[j * y_dim1 + 1], &c__1, &c_b8, &res[1], &c__1);
528 }
529 else if (y_prec_state__ == 1)
530 {
531 blas_sgemv_x_(trans_type__, n, n, &c_b6, &a[a_offset], lda, & y[j * y_dim1 + 1], &c__1, &c_b8, &res[1], &c__1, prec_type__);
532 }
533 else
534 {
535 blas_sgemv2_x_(trans_type__, n, n, &c_b6, &a[a_offset], lda, &y[j * y_dim1 + 1], &y_tail__[1], &c__1, &c_b8, &res[ 1], &c__1, prec_type__);
536 }
537 /* XXX: RES is no longer needed. */
538 scopy_(n, &res[1], &c__1, &dy[1], &c__1);
539 sgetrs_(trans, n, &c__1, &af[af_offset], ldaf, &ipiv[1], &dy[1], n, info);
540 /* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. */
541 normx = 0.f;
542 normy = 0.f;
543 normdx = 0.f;
544 dz_z__ = 0.f;
545 ymin = hugeval;
546 i__3 = *n;
547 for (i__ = 1;
548 i__ <= i__3;
549 ++i__)
550 {
551 yk = (r__1 = y[i__ + j * y_dim1], f2c_abs(r__1));
552 dyk = (r__1 = dy[i__], f2c_abs(r__1));
553 if (yk != 0.f)
554 {
555 /* Computing MAX */
556 r__1 = dz_z__;
557 r__2 = dyk / yk; // , expr subst
558 dz_z__ = max(r__1,r__2);
559 }
560 else if (dyk != 0.f)
561 {
562 dz_z__ = hugeval;
563 }
564 ymin = min(ymin,yk);
565 normy = max(normy,yk);
566 if (*colequ)
567 {
568 /* Computing MAX */
569 r__1 = normx;
570 r__2 = yk * c__[i__]; // , expr subst
571 normx = max(r__1,r__2);
572 /* Computing MAX */
573 r__1 = normdx;
574 r__2 = dyk * c__[i__]; // , expr subst
575 normdx = max(r__1,r__2);
576 }
577 else
578 {
579 normx = normy;
580 normdx = max(normdx,dyk);
581 }
582 }
583 if (normx != 0.f)
584 {
585 dx_x__ = normdx / normx;
586 }
587 else if (normdx == 0.f)
588 {
589 dx_x__ = 0.f;
590 }
591 else
592 {
593 dx_x__ = hugeval;
594 }
595 dxrat = normdx / prevnormdx;
596 dzrat = dz_z__ / prev_dz_z__;
597 /* Check termination criteria */
598 if (! (*ignore_cwise__) && ymin * *rcond < incr_thresh__ * normy && y_prec_state__ < 2)
599 {
600 incr_prec__ = TRUE_;
601 }
602 if (x_state__ == 3 && dxrat <= *rthresh)
603 {
604 x_state__ = 1;
605 }
606 if (x_state__ == 1)
607 {
608 if (dx_x__ <= eps)
609 {
610 x_state__ = 2;
611 }
612 else if (dxrat > *rthresh)
613 {
614 if (y_prec_state__ != 2)
615 {
616 incr_prec__ = TRUE_;
617 }
618 else
619 {
620 x_state__ = 3;
621 }
622 }
623 else
624 {
625 if (dxrat > dxratmax)
626 {
627 dxratmax = dxrat;
628 }
629 }
630 if (x_state__ > 1)
631 {
632 final_dx_x__ = dx_x__;
633 }
634 }
635 if (z_state__ == 0 && dz_z__ <= *dz_ub__)
636 {
637 z_state__ = 1;
638 }
639 if (z_state__ == 3 && dzrat <= *rthresh)
640 {
641 z_state__ = 1;
642 }
643 if (z_state__ == 1)
644 {
645 if (dz_z__ <= eps)
646 {
647 z_state__ = 2;
648 }
649 else if (dz_z__ > *dz_ub__)
650 {
651 z_state__ = 0;
652 dzratmax = 0.f;
653 final_dz_z__ = hugeval;
654 }
655 else if (dzrat > *rthresh)
656 {
657 if (y_prec_state__ != 2)
658 {
659 incr_prec__ = TRUE_;
660 }
661 else
662 {
663 z_state__ = 3;
664 }
665 }
666 else
667 {
668 if (dzrat > dzratmax)
669 {
670 dzratmax = dzrat;
671 }
672 }
673 if (z_state__ > 1)
674 {
675 final_dz_z__ = dz_z__;
676 }
677 }
678 /* Exit if both normwise and componentwise stopped working, */
679 /* but if componentwise is unstable, let it go at least two */
680 /* iterations. */
681 if (x_state__ != 1)
682 {
683 if (*ignore_cwise__)
684 {
685 goto L666;
686 }
687 if (z_state__ == 3 || z_state__ == 2)
688 {
689 goto L666;
690 }
691 if (z_state__ == 0 && cnt > 1)
692 {
693 goto L666;
694 }
695 }
696 if (incr_prec__)
697 {
698 incr_prec__ = FALSE_;
699 ++y_prec_state__;
700 i__3 = *n;
701 for (i__ = 1;
702 i__ <= i__3;
703 ++i__)
704 {
705 y_tail__[i__] = 0.f;
706 }
707 }
708 prevnormdx = normdx;
709 prev_dz_z__ = dz_z__;
710 /* Update soluton. */
711 if (y_prec_state__ < 2)
712 {
713 saxpy_(n, &c_b8, &dy[1], &c__1, &y[j * y_dim1 + 1], &c__1);
714 }
715 else
716 {
717 sla_wwaddw_(n, &y[j * y_dim1 + 1], &y_tail__[1], &dy[1]);
718 }
719 }
720 /* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't CALL F90_EXIT. */
721 L666: /* Set final_* when cnt hits ithresh. */
722 if (x_state__ == 1)
723 {
724 final_dx_x__ = dx_x__;
725 }
726 if (z_state__ == 1)
727 {
728 final_dz_z__ = dz_z__;
729 }
730 /* Compute error bounds */
731 if (*n_norms__ >= 1)
732 {
733 errs_n__[j + (errs_n_dim1 << 1)] = final_dx_x__ / (1 - dxratmax);
734 }
735 if (*n_norms__ >= 2)
736 {
737 errs_c__[j + (errs_c_dim1 << 1)] = final_dz_z__ / (1 - dzratmax);
738 }
739 /* Compute componentwise relative backward error from formula */
740 /* max(i) ( f2c_abs(R(i)) / ( f2c_abs(op(A_s))*f2c_abs(Y) + f2c_abs(B_s) )(i) ) */
741 /* where f2c_abs(Z) is the componentwise absolute value of the matrix */
742 /* or vector Z. */
743 /* Compute residual RES = B_s - op(A_s) * Y, */
744 /* op(A) = A, A**T, or A**H depending on TRANS (and type). */
745 scopy_(n, &b[j * b_dim1 + 1], &c__1, &res[1], &c__1);
746 sgemv_(trans, n, n, &c_b6, &a[a_offset], lda, &y[j * y_dim1 + 1], & c__1, &c_b8, &res[1], &c__1);
747 i__2 = *n;
748 for (i__ = 1;
749 i__ <= i__2;
750 ++i__)
751 {
752 ayb[i__] = (r__1 = b[i__ + j * b_dim1], f2c_abs(r__1));
753 }
754 /* Compute f2c_abs(op(A_s))*f2c_abs(Y) + f2c_abs(B_s). */
755 sla_geamv_(trans_type__, n, n, &c_b8, &a[a_offset], lda, &y[j * y_dim1 + 1], &c__1, &c_b8, &ayb[1], &c__1);
756 sla_lin_berr_(n, n, &c__1, &res[1], &ayb[1], &berr_out__[j]);
757 /* End of loop for each RHS. */
758 }
759 return 0;
760 }
761 /* sla_gerfsx_extended__ */
762 #endif
763