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