LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 responsible 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 P
70*> = 'S': Single
71*> = 'D': Double
72*> = 'I': Indigenous
73*> = 'X' or '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 T
81*> = '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 < 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 define 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*> \ingroup doubleGBcomputational
402*
403* =====================================================================
404 SUBROUTINE dla_gbrfsx_extended( PREC_TYPE, TRANS_TYPE, N, KL, KU,
405 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV,
406 $ COLEQU, C, B, LDB, Y, LDY,
407 $ BERR_OUT, N_NORMS, ERR_BNDS_NORM,
408 $ ERR_BNDS_COMP, RES, AYB, DY,
409 $ Y_TAIL, RCOND, ITHRESH, RTHRESH,
410 $ DZ_UB, IGNORE_CWISE, INFO )
411*
412* -- LAPACK computational routine --
413* -- LAPACK is a software package provided by Univ. of Tennessee, --
414* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
415*
416* .. Scalar Arguments ..
417 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS,
418 $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH
419 LOGICAL COLEQU, IGNORE_CWISE
420 DOUBLE PRECISION RTHRESH, DZ_UB
421* ..
422* .. Array Arguments ..
423 INTEGER IPIV( * )
424 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
425 $ y( ldy, * ), res(*), dy(*), y_tail(*)
426 DOUBLE PRECISION C( * ), AYB(*), RCOND, BERR_OUT(*),
427 $ ERR_BNDS_NORM( NRHS, * ),
428 $ ERR_BNDS_COMP( NRHS, * )
429* ..
430*
431* =====================================================================
432*
433* .. Local Scalars ..
434 CHARACTER TRANS
435 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE
436 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
437 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
438 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
439 $ EPS, HUGEVAL, INCR_THRESH
440 LOGICAL INCR_PREC
441* ..
442* .. Parameters ..
443 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
444 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
445 $ EXTRA_Y
446 parameter( unstable_state = 0, working_state = 1,
447 $ conv_state = 2, noprog_state = 3 )
448 parameter( base_residual = 0, extra_residual = 1,
449 $ extra_y = 2 )
450 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
451 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
452 INTEGER CMP_ERR_I, PIV_GROWTH_I
453 PARAMETER ( FINAL_NRM_ERR_I = 1, final_cmp_err_i = 2,
454 $ berr_i = 3 )
455 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
456 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
457 $ piv_growth_i = 9 )
458 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
459 $ la_linrx_cwise_i
460 parameter( la_linrx_itref_i = 1,
461 $ la_linrx_ithresh_i = 2 )
462 parameter( la_linrx_cwise_i = 3 )
463 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
464 $ la_linrx_rcond_i
465 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
466 parameter( la_linrx_rcond_i = 3 )
467* ..
468* .. External Subroutines ..
469 EXTERNAL daxpy, dcopy, dgbtrs, dgbmv, blas_dgbmv_x,
470 $ blas_dgbmv2_x, dla_gbamv, dla_wwaddw, dlamch,
472 DOUBLE PRECISION DLAMCH
473 CHARACTER CHLA_TRANSTYPE
474* ..
475* .. Intrinsic Functions ..
476 INTRINSIC abs, max, min
477* ..
478* .. Executable Statements ..
479*
480 IF (info.NE.0) RETURN
481 trans = chla_transtype(trans_type)
482 eps = dlamch( 'Epsilon' )
483 hugeval = dlamch( 'Overflow' )
484* Force HUGEVAL to Inf
485 hugeval = hugeval * hugeval
486* Using HUGEVAL may lead to spurious underflows.
487 incr_thresh = dble( n ) * eps
488 m = kl+ku+1
489
490 DO j = 1, nrhs
491 y_prec_state = extra_residual
492 IF ( y_prec_state .EQ. extra_y ) THEN
493 DO i = 1, n
494 y_tail( i ) = 0.0d+0
495 END DO
496 END IF
497
498 dxrat = 0.0d+0
499 dxratmax = 0.0d+0
500 dzrat = 0.0d+0
501 dzratmax = 0.0d+0
502 final_dx_x = hugeval
503 final_dz_z = hugeval
504 prevnormdx = hugeval
505 prev_dz_z = hugeval
506 dz_z = hugeval
507 dx_x = hugeval
508
509 x_state = working_state
510 z_state = unstable_state
511 incr_prec = .false.
512
513 DO cnt = 1, ithresh
514*
515* Compute residual RES = B_s - op(A_s) * Y,
516* op(A) = A, A**T, or A**H depending on TRANS (and type).
517*
518 CALL dcopy( n, b( 1, j ), 1, res, 1 )
519 IF ( y_prec_state .EQ. base_residual ) THEN
520 CALL dgbmv( trans, m, n, kl, ku, -1.0d+0, ab, ldab,
521 $ y( 1, j ), 1, 1.0d+0, res, 1 )
522 ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
523 CALL blas_dgbmv_x( trans_type, n, n, kl, ku,
524 $ -1.0d+0, ab, ldab, y( 1, j ), 1, 1.0d+0, res, 1,
525 $ prec_type )
526 ELSE
527 CALL blas_dgbmv2_x( trans_type, n, n, kl, ku, -1.0d+0,
528 $ ab, ldab, y( 1, j ), y_tail, 1, 1.0d+0, res, 1,
529 $ prec_type )
530 END IF
531
532! XXX: RES is no longer needed.
533 CALL dcopy( n, res, 1, dy, 1 )
534 CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
535 $ info )
536*
537* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
538*
539 normx = 0.0d+0
540 normy = 0.0d+0
541 normdx = 0.0d+0
542 dz_z = 0.0d+0
543 ymin = hugeval
544
545 DO i = 1, n
546 yk = abs( y( i, j ) )
547 dyk = abs( dy( i ) )
548
549 IF ( yk .NE. 0.0d+0 ) THEN
550 dz_z = max( dz_z, dyk / yk )
551 ELSE IF ( dyk .NE. 0.0d+0 ) THEN
552 dz_z = hugeval
553 END IF
554
555 ymin = min( ymin, yk )
556
557 normy = max( normy, yk )
558
559 IF ( colequ ) THEN
560 normx = max( normx, yk * c( i ) )
561 normdx = max( normdx, dyk * c( i ) )
562 ELSE
563 normx = normy
564 normdx = max( normdx, dyk )
565 END IF
566 END DO
567
568 IF ( normx .NE. 0.0d+0 ) THEN
569 dx_x = normdx / normx
570 ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
571 dx_x = 0.0d+0
572 ELSE
573 dx_x = hugeval
574 END IF
575
576 dxrat = normdx / prevnormdx
577 dzrat = dz_z / prev_dz_z
578*
579* Check termination criteria.
580*
581 IF ( .NOT.ignore_cwise
582 $ .AND. ymin*rcond .LT. incr_thresh*normy
583 $ .AND. y_prec_state .LT. extra_y )
584 $ incr_prec = .true.
585
586 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
587 $ x_state = working_state
588 IF ( x_state .EQ. working_state ) THEN
589 IF ( dx_x .LE. eps ) THEN
590 x_state = conv_state
591 ELSE IF ( dxrat .GT. rthresh ) THEN
592 IF ( y_prec_state .NE. extra_y ) THEN
593 incr_prec = .true.
594 ELSE
595 x_state = noprog_state
596 END IF
597 ELSE
598 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
599 END IF
600 IF ( x_state .GT. working_state ) final_dx_x = dx_x
601 END IF
602
603 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
604 $ z_state = working_state
605 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
606 $ z_state = working_state
607 IF ( z_state .EQ. working_state ) THEN
608 IF ( dz_z .LE. eps ) THEN
609 z_state = conv_state
610 ELSE IF ( dz_z .GT. dz_ub ) THEN
611 z_state = unstable_state
612 dzratmax = 0.0d+0
613 final_dz_z = hugeval
614 ELSE IF ( dzrat .GT. rthresh ) THEN
615 IF ( y_prec_state .NE. extra_y ) THEN
616 incr_prec = .true.
617 ELSE
618 z_state = noprog_state
619 END IF
620 ELSE
621 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
622 END IF
623 IF ( z_state .GT. working_state ) final_dz_z = dz_z
624 END IF
625*
626* Exit if both normwise and componentwise stopped working,
627* but if componentwise is unstable, let it go at least two
628* iterations.
629*
630 IF ( x_state.NE.working_state ) THEN
631 IF ( ignore_cwise ) GOTO 666
632 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
633 $ GOTO 666
634 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
635 END IF
636
637 IF ( incr_prec ) THEN
638 incr_prec = .false.
639 y_prec_state = y_prec_state + 1
640 DO i = 1, n
641 y_tail( i ) = 0.0d+0
642 END DO
643 END IF
644
645 prevnormdx = normdx
646 prev_dz_z = dz_z
647*
648* Update soluton.
649*
650 IF (y_prec_state .LT. extra_y) THEN
651 CALL daxpy( n, 1.0d+0, dy, 1, y(1,j), 1 )
652 ELSE
653 CALL dla_wwaddw( n, y(1,j), y_tail, dy )
654 END IF
655
656 END DO
657* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
658 666 CONTINUE
659*
660* Set final_* when cnt hits ithresh.
661*
662 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
663 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
664*
665* Compute error bounds.
666*
667 IF ( n_norms .GE. 1 ) THEN
668 err_bnds_norm( j, la_linrx_err_i ) =
669 $ final_dx_x / (1 - dxratmax)
670 END IF
671 IF (n_norms .GE. 2) THEN
672 err_bnds_comp( j, la_linrx_err_i ) =
673 $ final_dz_z / (1 - dzratmax)
674 END IF
675*
676* Compute componentwise relative backward error from formula
677* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
678* where abs(Z) is the componentwise absolute value of the matrix
679* or vector Z.
680*
681* Compute residual RES = B_s - op(A_s) * Y,
682* op(A) = A, A**T, or A**H depending on TRANS (and type).
683*
684 CALL dcopy( n, b( 1, j ), 1, res, 1 )
685 CALL dgbmv(trans, n, n, kl, ku, -1.0d+0, ab, ldab, y(1,j),
686 $ 1, 1.0d+0, res, 1 )
687
688 DO i = 1, n
689 ayb( i ) = abs( b( i, j ) )
690 END DO
691*
692* Compute abs(op(A_s))*abs(Y) + abs(B_s).
693*
694 CALL dla_gbamv( trans_type, n, n, kl, ku, 1.0d+0,
695 $ ab, ldab, y(1, j), 1, 1.0d+0, ayb, 1 )
696
697 CALL dla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
698*
699* End of loop for each RHS
700*
701 END DO
702*
703 RETURN
704*
705* End of DLA_GBRFSX_EXTENDED
706*
707 END
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
character *1 function chla_transtype(TRANS)
CHLA_TRANSTYPE
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGBMV
Definition: dgbmv.f:185
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:138
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:185
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:101
subroutine dla_wwaddw(N, X, Y, W)
DLA_WWADDW adds a vector into a doubled-single vector.
Definition: dla_wwaddw.f:81