LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
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 (LDY,NRHS)
181 *> On entry, the solution matrix X, as computed by DGBTRS.
182 *> On exit, the improved solution matrix Y.
183 *> \endverbatim
184 *>
185 *> \param[in] LDY
186 *> \verbatim
187 *> LDY is INTEGER
188 *> The leading dimension of the array Y. LDY >= max(1,N).
189 *> \endverbatim
190 *>
191 *> \param[out] BERR_OUT
192 *> \verbatim
193 *> BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
194 *> On exit, BERR_OUT(j) contains the componentwise relative backward
195 *> error for right-hand-side j from the formula
196 *> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
197 *> where abs(Z) is the componentwise absolute value of the matrix
198 *> or vector Z. This is computed by DLA_LIN_BERR.
199 *> \endverbatim
200 *>
201 *> \param[in] N_NORMS
202 *> \verbatim
203 *> N_NORMS is INTEGER
204 *> Determines which error bounds to return (see ERR_BNDS_NORM
205 *> and ERR_BNDS_COMP).
206 *> If N_NORMS >= 1 return normwise error bounds.
207 *> If N_NORMS >= 2 return componentwise error bounds.
208 *> \endverbatim
209 *>
210 *> \param[in,out] ERR_BNDS_NORM
211 *> \verbatim
212 *> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
213 *> For each right-hand side, this array contains information about
214 *> various error bounds and condition numbers corresponding to the
215 *> normwise relative error, which is defined as follows:
216 *>
217 *> Normwise relative error in the ith solution vector:
218 *> max_j (abs(XTRUE(j,i) - X(j,i)))
219 *> ------------------------------
220 *> max_j abs(X(j,i))
221 *>
222 *> The array is indexed by the type of error information as described
223 *> below. There currently are up to three pieces of information
224 *> returned.
225 *>
226 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
227 *> right-hand side.
228 *>
229 *> The second index in ERR_BNDS_NORM(:,err) contains the following
230 *> three fields:
231 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
232 *> reciprocal condition number is less than the threshold
233 *> sqrt(n) * slamch('Epsilon').
234 *>
235 *> err = 2 "Guaranteed" error bound: The estimated forward error,
236 *> almost certainly within a factor of 10 of the true error
237 *> so long as the next entry is greater than the threshold
238 *> sqrt(n) * slamch('Epsilon'). This error bound should only
239 *> be trusted if the previous boolean is true.
240 *>
241 *> err = 3 Reciprocal condition number: Estimated normwise
242 *> reciprocal condition number. Compared with the threshold
243 *> sqrt(n) * slamch('Epsilon') to determine if the error
244 *> estimate is "guaranteed". These reciprocal condition
245 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
246 *> appropriately scaled matrix Z.
247 *> Let Z = S*A, where S scales each row by a power of the
248 *> radix so all absolute row sums of Z are approximately 1.
249 *>
250 *> This subroutine is only responsible for setting the second field
251 *> above.
252 *> See Lapack Working Note 165 for further details and extra
253 *> cautions.
254 *> \endverbatim
255 *>
256 *> \param[in,out] ERR_BNDS_COMP
257 *> \verbatim
258 *> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
259 *> For each right-hand side, this array contains information about
260 *> various error bounds and condition numbers corresponding to the
261 *> componentwise relative error, which is defined as follows:
262 *>
263 *> Componentwise relative error in the ith solution vector:
264 *> abs(XTRUE(j,i) - X(j,i))
265 *> max_j ----------------------
266 *> abs(X(j,i))
267 *>
268 *> The array is indexed by the right-hand side i (on which the
269 *> componentwise relative error depends), and the type of error
270 *> information as described below. There currently are up to three
271 *> pieces of information returned for each right-hand side. If
272 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
273 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
274 *> the first (:,N_ERR_BNDS) entries are returned.
275 *>
276 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
277 *> right-hand side.
278 *>
279 *> The second index in ERR_BNDS_COMP(:,err) contains the following
280 *> three fields:
281 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
282 *> reciprocal condition number is less than the threshold
283 *> sqrt(n) * slamch('Epsilon').
284 *>
285 *> err = 2 "Guaranteed" error bound: The estimated forward error,
286 *> almost certainly within a factor of 10 of the true error
287 *> so long as the next entry is greater than the threshold
288 *> sqrt(n) * slamch('Epsilon'). This error bound should only
289 *> be trusted if the previous boolean is true.
290 *>
291 *> err = 3 Reciprocal condition number: Estimated componentwise
292 *> reciprocal condition number. Compared with the threshold
293 *> sqrt(n) * slamch('Epsilon') to determine if the error
294 *> estimate is "guaranteed". These reciprocal condition
295 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
296 *> appropriately scaled matrix Z.
297 *> Let Z = S*(A*diag(x)), where x is the solution for the
298 *> current right-hand side and S scales each row of
299 *> A*diag(x) by a power of the radix so all absolute row
300 *> sums of Z are approximately 1.
301 *>
302 *> This subroutine is only responsible for setting the second field
303 *> above.
304 *> See Lapack Working Note 165 for further details and extra
305 *> cautions.
306 *> \endverbatim
307 *>
308 *> \param[in] RES
309 *> \verbatim
310 *> RES is DOUBLE PRECISION array, dimension (N)
311 *> Workspace to hold the intermediate residual.
312 *> \endverbatim
313 *>
314 *> \param[in] AYB
315 *> \verbatim
316 *> AYB is DOUBLE PRECISION array, dimension (N)
317 *> Workspace. This can be the same workspace passed for Y_TAIL.
318 *> \endverbatim
319 *>
320 *> \param[in] DY
321 *> \verbatim
322 *> DY is DOUBLE PRECISION array, dimension (N)
323 *> Workspace to hold the intermediate solution.
324 *> \endverbatim
325 *>
326 *> \param[in] Y_TAIL
327 *> \verbatim
328 *> Y_TAIL is DOUBLE PRECISION array, dimension (N)
329 *> Workspace to hold the trailing bits of the intermediate solution.
330 *> \endverbatim
331 *>
332 *> \param[in] RCOND
333 *> \verbatim
334 *> RCOND is DOUBLE PRECISION
335 *> Reciprocal scaled condition number. This is an estimate of the
336 *> reciprocal Skeel condition number of the matrix A after
337 *> equilibration (if done). If this is less than the machine
338 *> precision (in particular, if it is zero), the matrix is singular
339 *> to working precision. Note that the error may still be small even
340 *> if this number is very small and the matrix appears ill-
341 *> conditioned.
342 *> \endverbatim
343 *>
344 *> \param[in] ITHRESH
345 *> \verbatim
346 *> ITHRESH is INTEGER
347 *> The maximum number of residual computations allowed for
348 *> refinement. The default is 10. For 'aggressive' set to 100 to
349 *> permit convergence using approximate factorizations or
350 *> factorizations other than LU. If the factorization uses a
351 *> technique other than Gaussian elimination, the guarantees in
352 *> ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
353 *> \endverbatim
354 *>
355 *> \param[in] RTHRESH
356 *> \verbatim
357 *> RTHRESH is DOUBLE PRECISION
358 *> Determines when to stop refinement if the error estimate stops
359 *> decreasing. Refinement will stop when the next solution no longer
360 *> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
361 *> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
362 *> default value is 0.5. For 'aggressive' set to 0.9 to permit
363 *> convergence on extremely ill-conditioned matrices. See LAWN 165
364 *> for more details.
365 *> \endverbatim
366 *>
367 *> \param[in] DZ_UB
368 *> \verbatim
369 *> DZ_UB is DOUBLE PRECISION
370 *> Determines when to start considering componentwise convergence.
371 *> Componentwise convergence is only considered after each component
372 *> of the solution Y is stable, which we definte as the relative
373 *> change in each component being less than DZ_UB. The default value
374 *> is 0.25, requiring the first bit to be stable. See LAWN 165 for
375 *> more details.
376 *> \endverbatim
377 *>
378 *> \param[in] IGNORE_CWISE
379 *> \verbatim
380 *> IGNORE_CWISE is LOGICAL
381 *> If .TRUE. then ignore componentwise convergence. Default value
382 *> is .FALSE..
383 *> \endverbatim
384 *>
385 *> \param[out] INFO
386 *> \verbatim
387 *> INFO is INTEGER
388 *> = 0: Successful exit.
389 *> < 0: if INFO = -i, the ith argument to DGBTRS had an illegal
390 *> value
391 *> \endverbatim
392 *
393 * Authors:
394 * ========
395 *
396 *> \author Univ. of Tennessee
397 *> \author Univ. of California Berkeley
398 *> \author Univ. of Colorado Denver
399 *> \author NAG Ltd.
400 *
401 *> \date June 2017
402 *
403 *> \ingroup doubleGBcomputational
404 *
405 * =====================================================================
406  SUBROUTINE dla_gbrfsx_extended( PREC_TYPE, TRANS_TYPE, N, KL, KU,
407  $ NRHS, AB, LDAB, AFB, LDAFB, IPIV,
408  $ COLEQU, C, B, LDB, Y, LDY,
409  $ BERR_OUT, N_NORMS, ERR_BNDS_NORM,
410  $ ERR_BNDS_COMP, RES, AYB, DY,
411  $ Y_TAIL, RCOND, ITHRESH, RTHRESH,
412  $ DZ_UB, IGNORE_CWISE, INFO )
413 *
414 * -- LAPACK computational routine (version 3.7.1) --
415 * -- LAPACK is a software package provided by Univ. of Tennessee, --
416 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
417 * June 2017
418 *
419 * .. Scalar Arguments ..
420  INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS,
421  $ prec_type, trans_type, n_norms, ithresh
422  LOGICAL COLEQU, IGNORE_CWISE
423  DOUBLE PRECISION RTHRESH, DZ_UB
424 * ..
425 * .. Array Arguments ..
426  INTEGER IPIV( * )
427  DOUBLE PRECISION AB( ldab, * ), AFB( ldafb, * ), B( ldb, * ),
428  $ y( ldy, * ), res(*), dy(*), y_tail(*)
429  DOUBLE PRECISION C( * ), AYB(*), RCOND, BERR_OUT(*),
430  $ err_bnds_norm( nrhs, * ),
431  $ err_bnds_comp( nrhs, * )
432 * ..
433 *
434 * =====================================================================
435 *
436 * .. Local Scalars ..
437  CHARACTER TRANS
438  INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE
439  DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
440  $ dzrat, prevnormdx, prev_dz_z, dxratmax,
441  $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
442  $ eps, hugeval, incr_thresh
443  LOGICAL INCR_PREC
444 * ..
445 * .. Parameters ..
446  INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
447  $ noprog_state, base_residual, extra_residual,
448  $ extra_y
449  parameter( unstable_state = 0, working_state = 1,
450  $ conv_state = 2, noprog_state = 3 )
451  parameter( base_residual = 0, extra_residual = 1,
452  $ extra_y = 2 )
453  INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
454  INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
455  INTEGER CMP_ERR_I, PIV_GROWTH_I
456  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
457  $ berr_i = 3 )
458  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
459  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
460  $ piv_growth_i = 9 )
461  INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
462  $ la_linrx_cwise_i
463  parameter( la_linrx_itref_i = 1,
464  $ la_linrx_ithresh_i = 2 )
465  parameter( la_linrx_cwise_i = 3 )
466  INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
467  $ la_linrx_rcond_i
468  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
469  parameter( la_linrx_rcond_i = 3 )
470 * ..
471 * .. External Subroutines ..
472  EXTERNAL daxpy, dcopy, dgbtrs, dgbmv, blas_dgbmv_x,
473  $ blas_dgbmv2_x, dla_gbamv, dla_wwaddw, dlamch,
475  DOUBLE PRECISION DLAMCH
476  CHARACTER CHLA_TRANSTYPE
477 * ..
478 * .. Intrinsic Functions ..
479  INTRINSIC abs, max, min
480 * ..
481 * .. Executable Statements ..
482 *
483  IF (info.NE.0) RETURN
484  trans = chla_transtype(trans_type)
485  eps = dlamch( 'Epsilon' )
486  hugeval = dlamch( 'Overflow' )
487 * Force HUGEVAL to Inf
488  hugeval = hugeval * hugeval
489 * Using HUGEVAL may lead to spurious underflows.
490  incr_thresh = dble( n ) * eps
491  m = kl+ku+1
492 
493  DO j = 1, nrhs
494  y_prec_state = extra_residual
495  IF ( y_prec_state .EQ. extra_y ) THEN
496  DO i = 1, n
497  y_tail( i ) = 0.0d+0
498  END DO
499  END IF
500 
501  dxrat = 0.0d+0
502  dxratmax = 0.0d+0
503  dzrat = 0.0d+0
504  dzratmax = 0.0d+0
505  final_dx_x = hugeval
506  final_dz_z = hugeval
507  prevnormdx = hugeval
508  prev_dz_z = hugeval
509  dz_z = hugeval
510  dx_x = hugeval
511 
512  x_state = working_state
513  z_state = unstable_state
514  incr_prec = .false.
515 
516  DO cnt = 1, ithresh
517 *
518 * Compute residual RES = B_s - op(A_s) * Y,
519 * op(A) = A, A**T, or A**H depending on TRANS (and type).
520 *
521  CALL dcopy( n, b( 1, j ), 1, res, 1 )
522  IF ( y_prec_state .EQ. base_residual ) THEN
523  CALL dgbmv( trans, m, n, kl, ku, -1.0d+0, ab, ldab,
524  $ y( 1, j ), 1, 1.0d+0, res, 1 )
525  ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
526  CALL blas_dgbmv_x( trans_type, n, n, kl, ku,
527  $ -1.0d+0, ab, ldab, y( 1, j ), 1, 1.0d+0, res, 1,
528  $ prec_type )
529  ELSE
530  CALL blas_dgbmv2_x( trans_type, n, n, kl, ku, -1.0d+0,
531  $ ab, ldab, y( 1, j ), y_tail, 1, 1.0d+0, res, 1,
532  $ prec_type )
533  END IF
534 
535 ! XXX: RES is no longer needed.
536  CALL dcopy( n, res, 1, dy, 1 )
537  CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
538  $ info )
539 *
540 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
541 *
542  normx = 0.0d+0
543  normy = 0.0d+0
544  normdx = 0.0d+0
545  dz_z = 0.0d+0
546  ymin = hugeval
547 
548  DO i = 1, n
549  yk = abs( y( i, j ) )
550  dyk = abs( dy( i ) )
551 
552  IF ( yk .NE. 0.0d+0 ) THEN
553  dz_z = max( dz_z, dyk / yk )
554  ELSE IF ( dyk .NE. 0.0d+0 ) THEN
555  dz_z = hugeval
556  END IF
557 
558  ymin = min( ymin, yk )
559 
560  normy = max( normy, yk )
561 
562  IF ( colequ ) THEN
563  normx = max( normx, yk * c( i ) )
564  normdx = max( normdx, dyk * c( i ) )
565  ELSE
566  normx = normy
567  normdx = max( normdx, dyk )
568  END IF
569  END DO
570 
571  IF ( normx .NE. 0.0d+0 ) THEN
572  dx_x = normdx / normx
573  ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
574  dx_x = 0.0d+0
575  ELSE
576  dx_x = hugeval
577  END IF
578 
579  dxrat = normdx / prevnormdx
580  dzrat = dz_z / prev_dz_z
581 *
582 * Check termination criteria.
583 *
584  IF ( .NOT.ignore_cwise
585  $ .AND. ymin*rcond .LT. incr_thresh*normy
586  $ .AND. y_prec_state .LT. extra_y )
587  $ incr_prec = .true.
588 
589  IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
590  $ x_state = working_state
591  IF ( x_state .EQ. working_state ) THEN
592  IF ( dx_x .LE. eps ) THEN
593  x_state = conv_state
594  ELSE IF ( dxrat .GT. rthresh ) THEN
595  IF ( y_prec_state .NE. extra_y ) THEN
596  incr_prec = .true.
597  ELSE
598  x_state = noprog_state
599  END IF
600  ELSE
601  IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
602  END IF
603  IF ( x_state .GT. working_state ) final_dx_x = dx_x
604  END IF
605 
606  IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
607  $ z_state = working_state
608  IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
609  $ z_state = working_state
610  IF ( z_state .EQ. working_state ) THEN
611  IF ( dz_z .LE. eps ) THEN
612  z_state = conv_state
613  ELSE IF ( dz_z .GT. dz_ub ) THEN
614  z_state = unstable_state
615  dzratmax = 0.0d+0
616  final_dz_z = hugeval
617  ELSE IF ( dzrat .GT. rthresh ) THEN
618  IF ( y_prec_state .NE. extra_y ) THEN
619  incr_prec = .true.
620  ELSE
621  z_state = noprog_state
622  END IF
623  ELSE
624  IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
625  END IF
626  IF ( z_state .GT. working_state ) final_dz_z = dz_z
627  END IF
628 *
629 * Exit if both normwise and componentwise stopped working,
630 * but if componentwise is unstable, let it go at least two
631 * iterations.
632 *
633  IF ( x_state.NE.working_state ) THEN
634  IF ( ignore_cwise ) GOTO 666
635  IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
636  $ GOTO 666
637  IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
638  END IF
639 
640  IF ( incr_prec ) THEN
641  incr_prec = .false.
642  y_prec_state = y_prec_state + 1
643  DO i = 1, n
644  y_tail( i ) = 0.0d+0
645  END DO
646  END IF
647 
648  prevnormdx = normdx
649  prev_dz_z = dz_z
650 *
651 * Update soluton.
652 *
653  IF (y_prec_state .LT. extra_y) THEN
654  CALL daxpy( n, 1.0d+0, dy, 1, y(1,j), 1 )
655  ELSE
656  CALL dla_wwaddw( n, y(1,j), y_tail, dy )
657  END IF
658 
659  END DO
660 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
661  666 CONTINUE
662 *
663 * Set final_* when cnt hits ithresh.
664 *
665  IF ( x_state .EQ. working_state ) final_dx_x = dx_x
666  IF ( z_state .EQ. working_state ) final_dz_z = dz_z
667 *
668 * Compute error bounds.
669 *
670  IF ( n_norms .GE. 1 ) THEN
671  err_bnds_norm( j, la_linrx_err_i ) =
672  $ final_dx_x / (1 - dxratmax)
673  END IF
674  IF (n_norms .GE. 2) THEN
675  err_bnds_comp( j, la_linrx_err_i ) =
676  $ final_dz_z / (1 - dzratmax)
677  END IF
678 *
679 * Compute componentwise relative backward error from formula
680 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
681 * where abs(Z) is the componentwise absolute value of the matrix
682 * or vector Z.
683 *
684 * Compute residual RES = B_s - op(A_s) * Y,
685 * op(A) = A, A**T, or A**H depending on TRANS (and type).
686 *
687  CALL dcopy( n, b( 1, j ), 1, res, 1 )
688  CALL dgbmv(trans, n, n, kl, ku, -1.0d+0, ab, ldab, y(1,j),
689  $ 1, 1.0d+0, res, 1 )
690 
691  DO i = 1, n
692  ayb( i ) = abs( b( i, j ) )
693  END DO
694 *
695 * Compute abs(op(A_s))*abs(Y) + abs(B_s).
696 *
697  CALL dla_gbamv( trans_type, n, n, kl, ku, 1.0d+0,
698  $ ab, ldab, y(1, j), 1, 1.0d+0, ayb, 1 )
699 
700  CALL dla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
701 *
702 * End of loop for each RHS
703 *
704  END DO
705 *
706  RETURN
707  END
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
subroutine dla_gbamv(TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, INCX, BETA, Y, INCY)
DLA_GBAMV performs a matrix-vector operation to calculate error bounds.
Definition: dla_gbamv.f:187
subroutine dgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGBMV
Definition: dgbmv.f:187
subroutine dla_wwaddw(N, X, Y, W)
DLA_WWADDW adds a vector into a doubled-single vector.
Definition: dla_wwaddw.f:83
subroutine dla_gbrfsx_extended(PREC_TYPE, TRANS_TYPE, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
DLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded...
subroutine dgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
DGBTRS
Definition: dgbtrs.f:140
subroutine dla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
DLA_LIN_BERR computes a component-wise relative backward error.
Definition: dla_lin_berr.f:103
character *1 function chla_transtype(TRANS)
CHLA_TRANSTYPE