1*> \brief \b CGBRFSX
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGBRFSX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbrfsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbrfsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbrfsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGBRFSX( 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, RWORK,
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*       REAL               RCOND
32*       ..
33*       .. Array Arguments ..
34*       INTEGER            IPIV( * )
35*       COMPLEX            AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
36*      $                   X( LDX , * ),WORK( * )
37*       REAL               R( * ), C( * ), PARAMS( * ), BERR( * ),
38*      $                   ERR_BNDS_NORM( NRHS, * ),
39*      $                   ERR_BNDS_COMP( NRHS, * ), RWORK( * )
40*       ..
41*
42*
43*> \par Purpose:
44*  =============
45*>
46*> \verbatim
47*>
48*>    CGBRFSX 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)
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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
140*>     Details of the LU factorization of the band matrix A, as
141*>     computed by CGBTRF.  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 CGETRF; 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 REAL 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 REAL 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 COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
211*>     On entry, the solution matrix X, as computed by CGETRS.
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 REAL
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 REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.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 COMPLEX array, dimension (2*N)
390*> \endverbatim
391*>
392*> \param[out] RWORK
393*> \verbatim
394*>          RWORK is REAL array, dimension (2*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 complexGBcomputational
433*
434*  =====================================================================
435      SUBROUTINE CGBRFSX( 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, RWORK,
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      REAL               RCOND
450*     ..
451*     .. Array Arguments ..
452      INTEGER            IPIV( * )
453      COMPLEX            AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
454     $                   X( LDX , * ),WORK( * )
455      REAL               R( * ), C( * ), PARAMS( * ), BERR( * ),
456     $                   ERR_BNDS_NORM( NRHS, * ),
457     $                   ERR_BNDS_COMP( NRHS, * ), RWORK( * )
458*     ..
459*
460*  ==================================================================
461*
462*     .. Parameters ..
463      REAL               ZERO, ONE
464      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
465      REAL               ITREF_DEFAULT, ITHRESH_DEFAULT,
466     $                   COMPONENTWISE_DEFAULT
467      REAL               RTHRESH_DEFAULT, DZTHRESH_DEFAULT
468      PARAMETER          ( ITREF_DEFAULT = 1.0 )
469      PARAMETER          ( ITHRESH_DEFAULT = 10.0 )
470      PARAMETER          ( COMPONENTWISE_DEFAULT = 1.0 )
471      PARAMETER          ( RTHRESH_DEFAULT = 0.5 )
472      PARAMETER          ( DZTHRESH_DEFAULT = 0.25 )
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, IGNORE_CWISE
486      INTEGER            J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS,
487     $                   ITHRESH
488      REAL               ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND,
489     $                   CWISE_WRONG, RTHRESH, UNSTABLE_THRESH
490*     ..
491*     .. External Subroutines ..
492      EXTERNAL           XERBLA, CGBCON, CLA_GBRFSX_EXTENDED
493*     ..
494*     .. Intrinsic Functions ..
495      INTRINSIC          MAX, SQRT, TRANSFER
496*     ..
497*     .. External Functions ..
498      EXTERNAL           LSAME, ILATRANS, ILAPREC
499      EXTERNAL           SLAMCH, CLANGB, CLA_GBRCOND_X, CLA_GBRCOND_C
500      REAL               SLAMCH, CLANGB, CLA_GBRCOND_X, CLA_GBRCOND_C
501      LOGICAL            LSAME
502      INTEGER            ILATRANS, ILAPREC
503*     ..
504*     .. Executable Statements ..
505*
506*     Check the input parameters.
507*
508      INFO = 0
509      TRANS_TYPE = ILATRANS( TRANS )
510      REF_TYPE = INT( ITREF_DEFAULT )
511      IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
512         IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0 ) THEN
513            PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
514         ELSE
515            REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
516         END IF
517      END IF
518*
519*     Set default parameters.
520*
521      ILLRCOND_THRESH = REAL( N ) * SLAMCH( 'Epsilon' )
522      ITHRESH = INT( ITHRESH_DEFAULT )
523      RTHRESH = RTHRESH_DEFAULT
524      UNSTABLE_THRESH = DZTHRESH_DEFAULT
525      IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0
526*
527      IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
528         IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0 ) THEN
529            PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
530         ELSE
531            ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
532         END IF
533      END IF
534      IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
535         IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0 ) THEN
536            IF ( IGNORE_CWISE ) THEN
537               PARAMS( LA_LINRX_CWISE_I ) = 0.0
538            ELSE
539               PARAMS( LA_LINRX_CWISE_I ) = 1.0
540            END IF
541         ELSE
542            IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0
543         END IF
544      END IF
545      IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
546         N_NORMS = 0
547      ELSE IF ( IGNORE_CWISE ) THEN
548         N_NORMS = 1
549      ELSE
550         N_NORMS = 2
551      END IF
552*
553      NOTRAN = LSAME( TRANS, 'N' )
554      ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
555      COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
556*
557*     Test input parameters.
558*
559      IF( TRANS_TYPE.EQ.-1 ) THEN
560        INFO = -1
561      ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
562     $         .NOT.LSAME( EQUED, 'N' ) ) THEN
563        INFO = -2
564      ELSE IF( N.LT.0 ) THEN
565        INFO = -3
566      ELSE IF( KL.LT.0 ) THEN
567        INFO = -4
568      ELSE IF( KU.LT.0 ) THEN
569        INFO = -5
570      ELSE IF( NRHS.LT.0 ) THEN
571        INFO = -6
572      ELSE IF( LDAB.LT.KL+KU+1 ) THEN
573        INFO = -8
574      ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
575        INFO = -10
576      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
577        INFO = -13
578      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
579        INFO = -15
580      END IF
581      IF( INFO.NE.0 ) THEN
582        CALL XERBLA( 'CGBRFSX', -INFO )
583        RETURN
584      END IF
585*
586*     Quick return if possible.
587*
588      IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
589         RCOND = 1.0
590         DO J = 1, NRHS
591            BERR( J ) = 0.0
592            IF ( N_ERR_BNDS .GE. 1 ) THEN
593               ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
594               ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
595            END IF
596            IF ( N_ERR_BNDS .GE. 2 ) THEN
597               ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0
598               ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0
599            END IF
600            IF ( N_ERR_BNDS .GE. 3 ) THEN
601               ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0
602               ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0
603            END IF
604         END DO
605         RETURN
606      END IF
607*
608*     Default to failure.
609*
610      RCOND = 0.0
611      DO J = 1, NRHS
612         BERR( J ) = 1.0
613         IF ( N_ERR_BNDS .GE. 1 ) THEN
614            ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
615            ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
616         END IF
617         IF ( N_ERR_BNDS .GE. 2 ) THEN
618            ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
619            ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
620         END IF
621         IF ( N_ERR_BNDS .GE. 3 ) THEN
622            ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0
623            ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0
624         END IF
625      END DO
626*
627*     Compute the norm of A and the reciprocal of the condition
628*     number of A.
629*
630      IF( NOTRAN ) THEN
631         NORM = 'I'
632      ELSE
633         NORM = '1'
634      END IF
635      ANORM = CLANGB( NORM, N, KL, KU, AB, LDAB, RWORK )
636      CALL CGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND,
637     $     WORK, RWORK, INFO )
638*
639*     Perform refinement on each right-hand side
640*
641      IF ( REF_TYPE .NE. 0 .AND. INFO .EQ. 0 ) THEN
642
643         PREC_TYPE = ILAPREC( 'D' )
644
645         IF ( NOTRAN ) THEN
646            CALL CLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE,  N, KL, KU,
647     $           NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B,
648     $           LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
649     $           ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
650     $           TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
651     $           RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
652     $           INFO )
653         ELSE
654            CALL CLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE,  N, KL, KU,
655     $           NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B,
656     $           LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
657     $           ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
658     $           TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
659     $           RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
660     $           INFO )
661         END IF
662      END IF
663
664      ERR_LBND = MAX( 10.0, SQRT( REAL( N ) ) ) * SLAMCH( 'Epsilon' )
665      IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN
666*
667*     Compute scaled normwise condition number cond(A*C).
668*
669         IF ( COLEQU .AND. NOTRAN ) THEN
670            RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
671     $           LDAFB, IPIV, C, .TRUE., INFO, WORK, RWORK )
672         ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
673            RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
674     $           LDAFB, IPIV, R, .TRUE., INFO, WORK, RWORK )
675         ELSE
676            RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
677     $           LDAFB, IPIV, C, .FALSE., INFO, WORK, RWORK )
678         END IF
679         DO J = 1, NRHS
680*
681*     Cap the error at 1.0.
682*
683            IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
684     $           .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0)
685     $           ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
686*
687*     Threshold the error (see LAWN).
688*
689            IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
690               ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
691               ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0
692               IF ( INFO .LE. N ) INFO = N + J
693            ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
694     $              THEN
695               ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
696               ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
697            END IF
698*
699*     Save the condition number.
700*
701            IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
702               ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
703            END IF
704
705         END DO
706      END IF
707
708      IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
709*
710*     Compute componentwise condition number cond(A*diag(Y(:,J))) for
711*     each right-hand side using the current solution as an estimate of
712*     the true solution.  If the componentwise error estimate is too
713*     large, then the solution is a lousy estimate of truth and the
714*     estimated RCOND may be too optimistic.  To avoid misleading users,
715*     the inverse condition number is set to 0.0 when the estimated
716*     cwise error is at least CWISE_WRONG.
717*
718         CWISE_WRONG = SQRT( SLAMCH( 'Epsilon' ) )
719         DO J = 1, NRHS
720            IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
721     $     THEN
722               RCOND_TMP = CLA_GBRCOND_X( TRANS, N, KL, KU, AB, LDAB,
723     $              AFB, LDAFB, IPIV, X( 1, J ), INFO, WORK, RWORK )
724            ELSE
725               RCOND_TMP = 0.0
726            END IF
727*
728*     Cap the error at 1.0.
729*
730            IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
731     $           .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0 )
732     $           ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
733*
734*     Threshold the error (see LAWN).
735*
736            IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
737               ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
738               ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0
739               IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0
740     $              .AND. INFO.LT.N + J ) INFO = N + J
741            ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
742     $              .LT. ERR_LBND ) THEN
743               ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
744               ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
745            END IF
746*
747*     Save the condition number.
748*
749            IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
750               ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
751            END IF
752
753         END DO
754      END IF
755*
756      RETURN
757*
758*     End of CGBRFSX
759*
760      END
761