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