LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 .LT. 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 .LE. 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 .LT. 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 *> \date April 2012
433 *
434 *> \ingroup doubleGBcomputational
435 *
436 * =====================================================================
437  SUBROUTINE dgbrfsx( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB,
438  $ ldafb, ipiv, r, c, b, ldb, x, ldx, rcond,
439  $ berr, n_err_bnds, err_bnds_norm,
440  $ err_bnds_comp, nparams, params, work, iwork,
441  $ info )
442 *
443 * -- LAPACK computational routine (version 3.6.1) --
444 * -- LAPACK is a software package provided by Univ. of Tennessee, --
445 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
446 * April 2012
447 *
448 * .. Scalar Arguments ..
449  CHARACTER TRANS, EQUED
450  INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
451  $ nparams, n_err_bnds
452  DOUBLE PRECISION RCOND
453 * ..
454 * .. Array Arguments ..
455  INTEGER IPIV( * ), IWORK( * )
456  DOUBLE PRECISION AB( ldab, * ), AFB( ldafb, * ), B( ldb, * ),
457  $ x( ldx , * ),work( * )
458  DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
459  $ err_bnds_norm( nrhs, * ),
460  $ err_bnds_comp( nrhs, * )
461 * ..
462 *
463 * ==================================================================
464 *
465 * .. Parameters ..
466  DOUBLE PRECISION ZERO, ONE
467  parameter ( zero = 0.0d+0, one = 1.0d+0 )
468  DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
469  DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
470  DOUBLE PRECISION DZTHRESH_DEFAULT
471  parameter ( itref_default = 1.0d+0 )
472  parameter ( ithresh_default = 10.0d+0 )
473  parameter ( componentwise_default = 1.0d+0 )
474  parameter ( rthresh_default = 0.5d+0 )
475  parameter ( dzthresh_default = 0.25d+0 )
476  INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
477  $ la_linrx_cwise_i
478  parameter ( la_linrx_itref_i = 1,
479  $ la_linrx_ithresh_i = 2 )
480  parameter ( la_linrx_cwise_i = 3 )
481  INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
482  $ la_linrx_rcond_i
483  parameter ( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
484  parameter ( la_linrx_rcond_i = 3 )
485 * ..
486 * .. Local Scalars ..
487  CHARACTER(1) NORM
488  LOGICAL ROWEQU, COLEQU, NOTRAN
489  INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
490  INTEGER N_NORMS
491  DOUBLE PRECISION ANORM, RCOND_TMP
492  DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
493  LOGICAL IGNORE_CWISE
494  INTEGER ITHRESH
495  DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
496 * ..
497 * .. External Subroutines ..
498  EXTERNAL xerbla, dgbcon
499  EXTERNAL dla_gbrfsx_extended
500 * ..
501 * .. Intrinsic Functions ..
502  INTRINSIC max, sqrt
503 * ..
504 * .. External Functions ..
505  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
506  EXTERNAL dlamch, dlangb, dla_gbrcond
507  DOUBLE PRECISION DLAMCH, DLANGB, DLA_GBRCOND
508  LOGICAL LSAME
509  INTEGER BLAS_FPINFO_X
510  INTEGER ILATRANS, ILAPREC
511 * ..
512 * .. Executable Statements ..
513 *
514 * Check the input parameters.
515 *
516  info = 0
517  trans_type = ilatrans( trans )
518  ref_type = int( itref_default )
519  IF ( nparams .GE. la_linrx_itref_i ) THEN
520  IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
521  params( la_linrx_itref_i ) = itref_default
522  ELSE
523  ref_type = params( la_linrx_itref_i )
524  END IF
525  END IF
526 *
527 * Set default parameters.
528 *
529  illrcond_thresh = dble( n ) * dlamch( 'Epsilon' )
530  ithresh = int( ithresh_default )
531  rthresh = rthresh_default
532  unstable_thresh = dzthresh_default
533  ignore_cwise = componentwise_default .EQ. 0.0d+0
534 *
535  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
536  IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
537  params( la_linrx_ithresh_i ) = ithresh
538  ELSE
539  ithresh = int( params( la_linrx_ithresh_i ) )
540  END IF
541  END IF
542  IF ( nparams.GE.la_linrx_cwise_i ) THEN
543  IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
544  IF ( ignore_cwise ) THEN
545  params( la_linrx_cwise_i ) = 0.0d+0
546  ELSE
547  params( la_linrx_cwise_i ) = 1.0d+0
548  END IF
549  ELSE
550  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
551  END IF
552  END IF
553  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
554  n_norms = 0
555  ELSE IF ( ignore_cwise ) THEN
556  n_norms = 1
557  ELSE
558  n_norms = 2
559  END IF
560 *
561  notran = lsame( trans, 'N' )
562  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
563  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
564 *
565 * Test input parameters.
566 *
567  IF( trans_type.EQ.-1 ) THEN
568  info = -1
569  ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
570  $ .NOT.lsame( equed, 'N' ) ) THEN
571  info = -2
572  ELSE IF( n.LT.0 ) THEN
573  info = -3
574  ELSE IF( kl.LT.0 ) THEN
575  info = -4
576  ELSE IF( ku.LT.0 ) THEN
577  info = -5
578  ELSE IF( nrhs.LT.0 ) THEN
579  info = -6
580  ELSE IF( ldab.LT.kl+ku+1 ) THEN
581  info = -8
582  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
583  info = -10
584  ELSE IF( ldb.LT.max( 1, n ) ) THEN
585  info = -13
586  ELSE IF( ldx.LT.max( 1, n ) ) THEN
587  info = -15
588  END IF
589  IF( info.NE.0 ) THEN
590  CALL xerbla( 'DGBRFSX', -info )
591  RETURN
592  END IF
593 *
594 * Quick return if possible.
595 *
596  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
597  rcond = 1.0d+0
598  DO j = 1, nrhs
599  berr( j ) = 0.0d+0
600  IF ( n_err_bnds .GE. 1 ) THEN
601  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
602  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
603  END IF
604  IF ( n_err_bnds .GE. 2 ) THEN
605  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
606  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
607  END IF
608  IF ( n_err_bnds .GE. 3 ) THEN
609  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
610  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
611  END IF
612  END DO
613  RETURN
614  END IF
615 *
616 * Default to failure.
617 *
618  rcond = 0.0d+0
619  DO j = 1, nrhs
620  berr( j ) = 1.0d+0
621  IF ( n_err_bnds .GE. 1 ) THEN
622  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
623  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
624  END IF
625  IF ( n_err_bnds .GE. 2 ) THEN
626  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
627  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
628  END IF
629  IF ( n_err_bnds .GE. 3 ) THEN
630  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
631  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
632  END IF
633  END DO
634 *
635 * Compute the norm of A and the reciprocal of the condition
636 * number of A.
637 *
638  IF( notran ) THEN
639  norm = 'I'
640  ELSE
641  norm = '1'
642  END IF
643  anorm = dlangb( norm, n, kl, ku, ab, ldab, work )
644  CALL dgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
645  $ work, iwork, info )
646 *
647 * Perform refinement on each right-hand side
648 *
649  IF ( ref_type .NE. 0 .AND. info .EQ. 0 ) THEN
650 
651  prec_type = ilaprec( 'E' )
652 
653  IF ( notran ) THEN
654  CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
655  $ nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b,
656  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
657  $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
658  $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
659  $ ignore_cwise, info )
660  ELSE
661  CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
662  $ nrhs, ab, ldab, afb, ldafb, ipiv, rowequ, r, b,
663  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
664  $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
665  $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
666  $ ignore_cwise, info )
667  END IF
668  END IF
669 
670  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
671  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
672 *
673 * Compute scaled normwise condition number cond(A*C).
674 *
675  IF ( colequ .AND. notran ) THEN
676  rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
677  $ ldafb, ipiv, -1, c, info, work, iwork )
678  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
679  rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
680  $ ldafb, ipiv, -1, r, info, work, iwork )
681  ELSE
682  rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
683  $ ldafb, ipiv, 0, r, info, work, iwork )
684  END IF
685  DO j = 1, nrhs
686 *
687 * Cap the error at 1.0.
688 *
689  IF ( n_err_bnds .GE. la_linrx_err_i
690  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
691  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
692 *
693 * Threshold the error (see LAWN).
694 *
695  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
696  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
697  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
698  IF ( info .LE. n ) info = n + j
699  ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
700  $ THEN
701  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
702  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
703  END IF
704 *
705 * Save the condition number.
706 *
707  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
708  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
709  END IF
710 
711  END DO
712  END IF
713 
714  IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 2) THEN
715 *
716 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
717 * each right-hand side using the current solution as an estimate of
718 * the true solution. If the componentwise error estimate is too
719 * large, then the solution is a lousy estimate of truth and the
720 * estimated RCOND may be too optimistic. To avoid misleading users,
721 * the inverse condition number is set to 0.0 when the estimated
722 * cwise error is at least CWISE_WRONG.
723 *
724  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
725  DO j = 1, nrhs
726  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
727  $ THEN
728  rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
729  $ ldafb, ipiv, 1, x( 1, j ), info, work, iwork )
730  ELSE
731  rcond_tmp = 0.0d+0
732  END IF
733 *
734 * Cap the error at 1.0.
735 *
736  IF ( n_err_bnds .GE. la_linrx_err_i
737  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
738  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
739 *
740 * Threshold the error (see LAWN).
741 *
742  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
743  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
744  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
745  IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
746  $ .AND. info.LT.n + j ) info = n + j
747  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
748  $ .LT. err_lbnd ) THEN
749  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
750  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
751  END IF
752 *
753 * Save the condition number.
754 *
755  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
756  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
757  END IF
758 
759  END DO
760  END IF
761 *
762  RETURN
763 *
764 * End of DGBRFSX
765 *
766  END
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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 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:442
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
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:126
subroutine dgbcon(NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
DGBCON
Definition: dgbcon.f:148
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.
Definition: dla_gbrcond.f:172
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55