LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgbrfsx.f
Go to the documentation of this file.
1*> \brief \b DGBRFSX
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGBRFSX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbrfsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbrfsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbrfsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGBRFSX( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB,
22* LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND,
23* BERR, N_ERR_BNDS, ERR_BNDS_NORM,
24* ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK,
25* INFO )
26*
27* .. Scalar Arguments ..
28* CHARACTER TRANS, EQUED
29* INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
30* $ NPARAMS, N_ERR_BNDS
31* DOUBLE PRECISION RCOND
32* ..
33* .. Array Arguments ..
34* INTEGER IPIV( * ), IWORK( * )
35* DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
36* $ X( LDX , * ),WORK( * )
37* DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
38* $ ERR_BNDS_NORM( NRHS, * ),
39* $ ERR_BNDS_COMP( NRHS, * )
40* ..
41*
42*
43*> \par Purpose:
44* =============
45*>
46*> \verbatim
47*>
48*> DGBRFSX improves the computed solution to a system of linear
49*> equations and provides error bounds and backward error estimates
50*> for the solution. In addition to normwise error bound, the code
51*> provides maximum componentwise error bound if possible. See
52*> comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
53*> error bounds.
54*>
55*> The original system of linear equations may have been equilibrated
56*> before calling this routine, as described by arguments EQUED, R
57*> and C below. In this case, the solution and error bounds returned
58*> are for the original unequilibrated system.
59*> \endverbatim
60*
61* Arguments:
62* ==========
63*
64*> \verbatim
65*> Some optional parameters are bundled in the PARAMS array. These
66*> settings determine how refinement is performed, but often the
67*> defaults are acceptable. If the defaults are acceptable, users
68*> can pass NPARAMS = 0 which prevents the source code from accessing
69*> the PARAMS argument.
70*> \endverbatim
71*>
72*> \param[in] TRANS
73*> \verbatim
74*> TRANS is CHARACTER*1
75*> Specifies the form of the system of equations:
76*> = 'N': A * X = B (No transpose)
77*> = 'T': A**T * X = B (Transpose)
78*> = 'C': A**H * X = B (Conjugate transpose = Transpose)
79*> \endverbatim
80*>
81*> \param[in] EQUED
82*> \verbatim
83*> EQUED is CHARACTER*1
84*> Specifies the form of equilibration that was done to A
85*> before calling this routine. This is needed to compute
86*> the solution and error bounds correctly.
87*> = 'N': No equilibration
88*> = 'R': Row equilibration, i.e., A has been premultiplied by
89*> diag(R).
90*> = 'C': Column equilibration, i.e., A has been postmultiplied
91*> by diag(C).
92*> = 'B': Both row and column equilibration, i.e., A has been
93*> replaced by diag(R) * A * diag(C).
94*> The right hand side B has been changed accordingly.
95*> \endverbatim
96*>
97*> \param[in] N
98*> \verbatim
99*> N is INTEGER
100*> The order of the matrix A. N >= 0.
101*> \endverbatim
102*>
103*> \param[in] KL
104*> \verbatim
105*> KL is INTEGER
106*> The number of subdiagonals within the band of A. KL >= 0.
107*> \endverbatim
108*>
109*> \param[in] KU
110*> \verbatim
111*> KU is INTEGER
112*> The number of superdiagonals within the band of A. KU >= 0.
113*> \endverbatim
114*>
115*> \param[in] NRHS
116*> \verbatim
117*> NRHS is INTEGER
118*> The number of right hand sides, i.e., the number of columns
119*> of the matrices B and X. NRHS >= 0.
120*> \endverbatim
121*>
122*> \param[in] AB
123*> \verbatim
124*> AB is DOUBLE PRECISION array, dimension (LDAB,N)
125*> The original band matrix A, stored in rows 1 to KL+KU+1.
126*> The j-th column of A is stored in the j-th column of the
127*> array AB as follows:
128*> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
129*> \endverbatim
130*>
131*> \param[in] LDAB
132*> \verbatim
133*> LDAB is INTEGER
134*> The leading dimension of the array AB. LDAB >= KL+KU+1.
135*> \endverbatim
136*>
137*> \param[in] AFB
138*> \verbatim
139*> AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
140*> Details of the LU factorization of the band matrix A, as
141*> computed by DGBTRF. U is stored as an upper triangular band
142*> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
143*> the multipliers used during the factorization are stored in
144*> rows KL+KU+2 to 2*KL+KU+1.
145*> \endverbatim
146*>
147*> \param[in] LDAFB
148*> \verbatim
149*> LDAFB is INTEGER
150*> The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
151*> \endverbatim
152*>
153*> \param[in] IPIV
154*> \verbatim
155*> IPIV is INTEGER array, dimension (N)
156*> The pivot indices from DGETRF; for 1<=i<=N, row i of the
157*> matrix was interchanged with row IPIV(i).
158*> \endverbatim
159*>
160*> \param[in,out] R
161*> \verbatim
162*> R is DOUBLE PRECISION array, dimension (N)
163*> The row scale factors for A. If EQUED = 'R' or 'B', A is
164*> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
165*> is not accessed. R is an input argument if FACT = 'F';
166*> otherwise, R is an output argument. If FACT = 'F' and
167*> EQUED = 'R' or 'B', each element of R must be positive.
168*> If R is output, each element of R is a power of the radix.
169*> If R is input, each element of R should be a power of the radix
170*> to ensure a reliable solution and error estimates. Scaling by
171*> powers of the radix does not cause rounding errors unless the
172*> result underflows or overflows. Rounding errors during scaling
173*> lead to refining with a matrix that is not equivalent to the
174*> input matrix, producing error estimates that may not be
175*> reliable.
176*> \endverbatim
177*>
178*> \param[in,out] C
179*> \verbatim
180*> C is DOUBLE PRECISION array, dimension (N)
181*> The column scale factors for A. If EQUED = 'C' or 'B', A is
182*> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
183*> is not accessed. C is an input argument if FACT = 'F';
184*> otherwise, C is an output argument. If FACT = 'F' and
185*> EQUED = 'C' or 'B', each element of C must be positive.
186*> If C is output, each element of C is a power of the radix.
187*> If C is input, each element of C should be a power of the radix
188*> to ensure a reliable solution and error estimates. Scaling by
189*> powers of the radix does not cause rounding errors unless the
190*> result underflows or overflows. Rounding errors during scaling
191*> lead to refining with a matrix that is not equivalent to the
192*> input matrix, producing error estimates that may not be
193*> reliable.
194*> \endverbatim
195*>
196*> \param[in] B
197*> \verbatim
198*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
199*> The right hand side matrix B.
200*> \endverbatim
201*>
202*> \param[in] LDB
203*> \verbatim
204*> LDB is INTEGER
205*> The leading dimension of the array B. LDB >= max(1,N).
206*> \endverbatim
207*>
208*> \param[in,out] X
209*> \verbatim
210*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
211*> On entry, the solution matrix X, as computed by DGETRS.
212*> On exit, the improved solution matrix X.
213*> \endverbatim
214*>
215*> \param[in] LDX
216*> \verbatim
217*> LDX is INTEGER
218*> The leading dimension of the array X. LDX >= max(1,N).
219*> \endverbatim
220*>
221*> \param[out] RCOND
222*> \verbatim
223*> RCOND is DOUBLE PRECISION
224*> Reciprocal scaled condition number. This is an estimate of the
225*> reciprocal Skeel condition number of the matrix A after
226*> equilibration (if done). If this is less than the machine
227*> precision (in particular, if it is zero), the matrix is singular
228*> to working precision. Note that the error may still be small even
229*> if this number is very small and the matrix appears ill-
230*> conditioned.
231*> \endverbatim
232*>
233*> \param[out] BERR
234*> \verbatim
235*> BERR is DOUBLE PRECISION array, dimension (NRHS)
236*> Componentwise relative backward error. This is the
237*> componentwise relative backward error of each solution vector X(j)
238*> (i.e., the smallest relative change in any element of A or B that
239*> makes X(j) an exact solution).
240*> \endverbatim
241*>
242*> \param[in] N_ERR_BNDS
243*> \verbatim
244*> N_ERR_BNDS is INTEGER
245*> Number of error bounds to return for each right hand side
246*> and each type (normwise or componentwise). See ERR_BNDS_NORM and
247*> ERR_BNDS_COMP below.
248*> \endverbatim
249*>
250*> \param[out] ERR_BNDS_NORM
251*> \verbatim
252*> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
253*> For each right-hand side, this array contains information about
254*> various error bounds and condition numbers corresponding to the
255*> normwise relative error, which is defined as follows:
256*>
257*> Normwise relative error in the ith solution vector:
258*> max_j (abs(XTRUE(j,i) - X(j,i)))
259*> ------------------------------
260*> max_j abs(X(j,i))
261*>
262*> The array is indexed by the type of error information as described
263*> below. There currently are up to three pieces of information
264*> returned.
265*>
266*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
267*> right-hand side.
268*>
269*> The second index in ERR_BNDS_NORM(:,err) contains the following
270*> three fields:
271*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
272*> reciprocal condition number is less than the threshold
273*> sqrt(n) * dlamch('Epsilon').
274*>
275*> err = 2 "Guaranteed" error bound: The estimated forward error,
276*> almost certainly within a factor of 10 of the true error
277*> so long as the next entry is greater than the threshold
278*> sqrt(n) * dlamch('Epsilon'). This error bound should only
279*> be trusted if the previous boolean is true.
280*>
281*> err = 3 Reciprocal condition number: Estimated normwise
282*> reciprocal condition number. Compared with the threshold
283*> sqrt(n) * dlamch('Epsilon') to determine if the error
284*> estimate is "guaranteed". These reciprocal condition
285*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
286*> appropriately scaled matrix Z.
287*> Let Z = S*A, where S scales each row by a power of the
288*> radix so all absolute row sums of Z are approximately 1.
289*>
290*> See Lapack Working Note 165 for further details and extra
291*> cautions.
292*> \endverbatim
293*>
294*> \param[out] ERR_BNDS_COMP
295*> \verbatim
296*> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
297*> For each right-hand side, this array contains information about
298*> various error bounds and condition numbers corresponding to the
299*> componentwise relative error, which is defined as follows:
300*>
301*> Componentwise relative error in the ith solution vector:
302*> abs(XTRUE(j,i) - X(j,i))
303*> max_j ----------------------
304*> abs(X(j,i))
305*>
306*> The array is indexed by the right-hand side i (on which the
307*> componentwise relative error depends), and the type of error
308*> information as described below. There currently are up to three
309*> pieces of information returned for each right-hand side. If
310*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
311*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
312*> the first (:,N_ERR_BNDS) entries are returned.
313*>
314*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
315*> right-hand side.
316*>
317*> The second index in ERR_BNDS_COMP(:,err) contains the following
318*> three fields:
319*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
320*> reciprocal condition number is less than the threshold
321*> sqrt(n) * dlamch('Epsilon').
322*>
323*> err = 2 "Guaranteed" error bound: The estimated forward error,
324*> almost certainly within a factor of 10 of the true error
325*> so long as the next entry is greater than the threshold
326*> sqrt(n) * dlamch('Epsilon'). This error bound should only
327*> be trusted if the previous boolean is true.
328*>
329*> err = 3 Reciprocal condition number: Estimated componentwise
330*> reciprocal condition number. Compared with the threshold
331*> sqrt(n) * dlamch('Epsilon') to determine if the error
332*> estimate is "guaranteed". These reciprocal condition
333*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
334*> appropriately scaled matrix Z.
335*> Let Z = S*(A*diag(x)), where x is the solution for the
336*> current right-hand side and S scales each row of
337*> A*diag(x) by a power of the radix so all absolute row
338*> sums of Z are approximately 1.
339*>
340*> See Lapack Working Note 165 for further details and extra
341*> cautions.
342*> \endverbatim
343*>
344*> \param[in] NPARAMS
345*> \verbatim
346*> NPARAMS is INTEGER
347*> Specifies the number of parameters set in PARAMS. If <= 0, the
348*> PARAMS array is never referenced and default values are used.
349*> \endverbatim
350*>
351*> \param[in,out] PARAMS
352*> \verbatim
353*> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
354*> Specifies algorithm parameters. If an entry is < 0.0, then
355*> that entry will be filled with default value used for that
356*> parameter. Only positions up to NPARAMS are accessed; defaults
357*> are used for higher-numbered parameters.
358*>
359*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
360*> refinement or not.
361*> Default: 1.0D+0
362*> = 0.0: No refinement is performed, and no error bounds are
363*> computed.
364*> = 1.0: Use the double-precision refinement algorithm,
365*> possibly with doubled-single computations if the
366*> compilation environment does not support DOUBLE
367*> PRECISION.
368*> (other values are reserved for future use)
369*>
370*> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
371*> computations allowed for refinement.
372*> Default: 10
373*> Aggressive: Set to 100 to permit convergence using approximate
374*> factorizations or factorizations other than LU. If
375*> the factorization uses a technique other than
376*> Gaussian elimination, the guarantees in
377*> err_bnds_norm and err_bnds_comp may no longer be
378*> trustworthy.
379*>
380*> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
381*> will attempt to find a solution with small componentwise
382*> relative error in the double-precision algorithm. Positive
383*> is true, 0.0 is false.
384*> Default: 1.0 (attempt componentwise convergence)
385*> \endverbatim
386*>
387*> \param[out] WORK
388*> \verbatim
389*> WORK is DOUBLE PRECISION array, dimension (4*N)
390*> \endverbatim
391*>
392*> \param[out] IWORK
393*> \verbatim
394*> IWORK is INTEGER array, dimension (N)
395*> \endverbatim
396*>
397*> \param[out] INFO
398*> \verbatim
399*> INFO is INTEGER
400*> = 0: Successful exit. The solution to every right-hand side is
401*> guaranteed.
402*> < 0: If INFO = -i, the i-th argument had an illegal value
403*> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
404*> has been completed, but the factor U is exactly singular, so
405*> the solution and error bounds could not be computed. RCOND = 0
406*> is returned.
407*> = N+J: The solution corresponding to the Jth right-hand side is
408*> not guaranteed. The solutions corresponding to other right-
409*> hand sides K with K > J may not be guaranteed as well, but
410*> only the first such right-hand side is reported. If a small
411*> componentwise error is not requested (PARAMS(3) = 0.0) then
412*> the Jth right-hand side is the first with a normwise error
413*> bound that is not guaranteed (the smallest J such
414*> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
415*> the Jth right-hand side is the first with either a normwise or
416*> componentwise error bound that is not guaranteed (the smallest
417*> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
418*> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
419*> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
420*> about all of the right-hand sides check ERR_BNDS_NORM or
421*> ERR_BNDS_COMP.
422*> \endverbatim
423*
424* Authors:
425* ========
426*
427*> \author Univ. of Tennessee
428*> \author Univ. of California Berkeley
429*> \author Univ. of Colorado Denver
430*> \author NAG Ltd.
431*
432*> \ingroup gbrfsx
433*
434* =====================================================================
435 SUBROUTINE dgbrfsx( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB,
436 $ LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND,
437 $ BERR, N_ERR_BNDS, ERR_BNDS_NORM,
438 $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK,
439 $ INFO )
440*
441* -- LAPACK computational routine --
442* -- LAPACK is a software package provided by Univ. of Tennessee, --
443* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
444*
445* .. Scalar Arguments ..
446 CHARACTER TRANS, EQUED
447 INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
448 $ NPARAMS, N_ERR_BNDS
449 DOUBLE PRECISION RCOND
450* ..
451* .. Array Arguments ..
452 INTEGER IPIV( * ), IWORK( * )
453 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
454 $ X( LDX , * ),WORK( * )
455 DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
456 $ err_bnds_norm( nrhs, * ),
457 $ err_bnds_comp( nrhs, * )
458* ..
459*
460* ==================================================================
461*
462* .. Parameters ..
463 DOUBLE PRECISION ZERO, ONE
464 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
465 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
466 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
467 DOUBLE PRECISION DZTHRESH_DEFAULT
468 parameter( itref_default = 1.0d+0 )
469 parameter( ithresh_default = 10.0d+0 )
470 parameter( componentwise_default = 1.0d+0 )
471 parameter( rthresh_default = 0.5d+0 )
472 parameter( dzthresh_default = 0.25d+0 )
473 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
474 $ la_linrx_cwise_i
475 parameter( la_linrx_itref_i = 1,
476 $ la_linrx_ithresh_i = 2 )
477 parameter( la_linrx_cwise_i = 3 )
478 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
479 $ la_linrx_rcond_i
480 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
481 parameter( la_linrx_rcond_i = 3 )
482* ..
483* .. Local Scalars ..
484 CHARACTER(1) NORM
485 LOGICAL ROWEQU, COLEQU, NOTRAN
486 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
487 INTEGER N_NORMS
488 DOUBLE PRECISION ANORM, RCOND_TMP
489 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
490 LOGICAL IGNORE_CWISE
491 INTEGER ITHRESH
492 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
493* ..
494* .. External Subroutines ..
495 EXTERNAL xerbla, dgbcon
496 EXTERNAL dla_gbrfsx_extended
497* ..
498* .. Intrinsic Functions ..
499 INTRINSIC max, sqrt
500* ..
501* .. External Functions ..
502 EXTERNAL lsame, ilatrans, ilaprec
503 EXTERNAL dlamch, dlangb, dla_gbrcond
504 DOUBLE PRECISION DLAMCH, DLANGB, DLA_GBRCOND
505 LOGICAL LSAME
506 INTEGER ILATRANS, ILAPREC
507* ..
508* .. Executable Statements ..
509*
510* Check the input parameters.
511*
512 info = 0
513 trans_type = ilatrans( trans )
514 ref_type = int( itref_default )
515 IF ( nparams .GE. la_linrx_itref_i ) THEN
516 IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
517 params( la_linrx_itref_i ) = itref_default
518 ELSE
519 ref_type = params( la_linrx_itref_i )
520 END IF
521 END IF
522*
523* Set default parameters.
524*
525 illrcond_thresh = dble( n ) * dlamch( 'Epsilon' )
526 ithresh = int( ithresh_default )
527 rthresh = rthresh_default
528 unstable_thresh = dzthresh_default
529 ignore_cwise = componentwise_default .EQ. 0.0d+0
530*
531 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
532 IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
533 params( la_linrx_ithresh_i ) = ithresh
534 ELSE
535 ithresh = int( params( la_linrx_ithresh_i ) )
536 END IF
537 END IF
538 IF ( nparams.GE.la_linrx_cwise_i ) THEN
539 IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
540 IF ( ignore_cwise ) THEN
541 params( la_linrx_cwise_i ) = 0.0d+0
542 ELSE
543 params( la_linrx_cwise_i ) = 1.0d+0
544 END IF
545 ELSE
546 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
547 END IF
548 END IF
549 IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
550 n_norms = 0
551 ELSE IF ( ignore_cwise ) THEN
552 n_norms = 1
553 ELSE
554 n_norms = 2
555 END IF
556*
557 notran = lsame( trans, 'N' )
558 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
559 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
560*
561* Test input parameters.
562*
563 IF( trans_type.EQ.-1 ) THEN
564 info = -1
565 ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
566 $ .NOT.lsame( equed, 'N' ) ) THEN
567 info = -2
568 ELSE IF( n.LT.0 ) THEN
569 info = -3
570 ELSE IF( kl.LT.0 ) THEN
571 info = -4
572 ELSE IF( ku.LT.0 ) THEN
573 info = -5
574 ELSE IF( nrhs.LT.0 ) THEN
575 info = -6
576 ELSE IF( ldab.LT.kl+ku+1 ) THEN
577 info = -8
578 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
579 info = -10
580 ELSE IF( ldb.LT.max( 1, n ) ) THEN
581 info = -13
582 ELSE IF( ldx.LT.max( 1, n ) ) THEN
583 info = -15
584 END IF
585 IF( info.NE.0 ) THEN
586 CALL xerbla( 'DGBRFSX', -info )
587 RETURN
588 END IF
589*
590* Quick return if possible.
591*
592 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
593 rcond = 1.0d+0
594 DO j = 1, nrhs
595 berr( j ) = 0.0d+0
596 IF ( n_err_bnds .GE. 1 ) THEN
597 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
598 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
599 END IF
600 IF ( n_err_bnds .GE. 2 ) THEN
601 err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
602 err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
603 END IF
604 IF ( n_err_bnds .GE. 3 ) THEN
605 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
606 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
607 END IF
608 END DO
609 RETURN
610 END IF
611*
612* Default to failure.
613*
614 rcond = 0.0d+0
615 DO j = 1, nrhs
616 berr( j ) = 1.0d+0
617 IF ( n_err_bnds .GE. 1 ) THEN
618 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
619 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
620 END IF
621 IF ( n_err_bnds .GE. 2 ) THEN
622 err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
623 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
624 END IF
625 IF ( n_err_bnds .GE. 3 ) THEN
626 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
627 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
628 END IF
629 END DO
630*
631* Compute the norm of A and the reciprocal of the condition
632* number of A.
633*
634 IF( notran ) THEN
635 norm = 'I'
636 ELSE
637 norm = '1'
638 END IF
639 anorm = dlangb( norm, n, kl, ku, ab, ldab, work )
640 CALL dgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
641 $ work, iwork, info )
642*
643* Perform refinement on each right-hand side
644*
645 IF ( ref_type .NE. 0 .AND. info .EQ. 0 ) THEN
646
647 prec_type = ilaprec( 'E' )
648
649 IF ( notran ) THEN
650 CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
651 $ nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b,
652 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
653 $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
654 $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
655 $ ignore_cwise, info )
656 ELSE
657 CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
658 $ nrhs, ab, ldab, afb, ldafb, ipiv, rowequ, r, b,
659 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
660 $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
661 $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
662 $ ignore_cwise, info )
663 END IF
664 END IF
665
666 err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
667 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
668*
669* Compute scaled normwise condition number cond(A*C).
670*
671 IF ( colequ .AND. notran ) THEN
672 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
673 $ ldafb, ipiv, -1, c, info, work, iwork )
674 ELSE IF ( rowequ .AND. .NOT. notran ) THEN
675 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
676 $ ldafb, ipiv, -1, r, info, work, iwork )
677 ELSE
678 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
679 $ ldafb, ipiv, 0, r, info, work, iwork )
680 END IF
681 DO j = 1, nrhs
682*
683* Cap the error at 1.0.
684*
685 IF ( n_err_bnds .GE. la_linrx_err_i
686 $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
687 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
688*
689* Threshold the error (see LAWN).
690*
691 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
692 err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
693 err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
694 IF ( info .LE. n ) info = n + j
695 ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
696 $ THEN
697 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
698 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
699 END IF
700*
701* Save the condition number.
702*
703 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
704 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
705 END IF
706
707 END DO
708 END IF
709
710 IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 2) THEN
711*
712* Compute componentwise condition number cond(A*diag(Y(:,J))) for
713* each right-hand side using the current solution as an estimate of
714* the true solution. If the componentwise error estimate is too
715* large, then the solution is a lousy estimate of truth and the
716* estimated RCOND may be too optimistic. To avoid misleading users,
717* the inverse condition number is set to 0.0 when the estimated
718* cwise error is at least CWISE_WRONG.
719*
720 cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
721 DO j = 1, nrhs
722 IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
723 $ THEN
724 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
725 $ ldafb, ipiv, 1, x( 1, j ), info, work, iwork )
726 ELSE
727 rcond_tmp = 0.0d+0
728 END IF
729*
730* Cap the error at 1.0.
731*
732 IF ( n_err_bnds .GE. la_linrx_err_i
733 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
734 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
735*
736* Threshold the error (see LAWN).
737*
738 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
739 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
740 err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
741 IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
742 $ .AND. info.LT.n + j ) info = n + j
743 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
744 $ .LT. err_lbnd ) THEN
745 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
746 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
747 END IF
748*
749* Save the condition number.
750*
751 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
752 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
753 END IF
754
755 END DO
756 END IF
757*
758 RETURN
759*
760* End of DGBRFSX
761*
762 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
DGBCON
Definition dgbcon.f:146
subroutine dgbrfsx(trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DGBRFSX
Definition dgbrfsx.f:440
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
integer function ilatrans(trans)
ILATRANS
Definition ilatrans.f:58
double precision function dla_gbrcond(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
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...
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlangb(norm, n, kl, ku, ab, ldab, work)
DLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlangb.f:124
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48