LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zgbrfsx.f
Go to the documentation of this file.
1 *> \brief \b ZGBRFSX
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGBRFSX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbrfsx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbrfsx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbrfsx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGBRFSX( 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 * DOUBLE PRECISION RCOND
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IPIV( * )
35 * COMPLEX*16 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, * ), RWORK( * )
40 * ..
41 *
42 *
43 *> \par Purpose:
44 * =============
45 *>
46 *> \verbatim
47 *>
48 *> ZGBRFSX 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 COMPLEX*16 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*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (2*N)
390 *> \endverbatim
391 *>
392 *> \param[out] RWORK
393 *> \verbatim
394 *> RWORK is DOUBLE PRECISION 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 *> \date April 2012
433 *
434 *> \ingroup complex16GBcomputational
435 *
436 * =====================================================================
437  SUBROUTINE zgbrfsx( 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, rwork,
441  $ info )
442 *
443 * -- LAPACK computational routine (version 3.4.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( * )
456  COMPLEX*16 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, * ), rwork( * )
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, ignore_cwise
489  INTEGER j, trans_type, prec_type, ref_type, n_norms,
490  $ ithresh
491  DOUBLE PRECISION anorm, rcond_tmp, illrcond_thresh, err_lbnd,
492  $ cwise_wrong, rthresh, unstable_thresh
493 * ..
494 * .. External Subroutines ..
496 * ..
497 * .. Intrinsic Functions ..
498  INTRINSIC max, sqrt, transfer
499 * ..
500 * .. External Functions ..
501  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
503  DOUBLE PRECISION dlamch, zlangb, zla_gbrcond_x, zla_gbrcond_c
504  LOGICAL lsame
505  INTEGER blas_fpinfo_x
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( 'ZGBRFSX', -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 = zlangb( norm, n, kl, ku, ab, ldab, rwork )
640  CALL zgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
641  $ work, rwork, info )
642 *
643 * Perform refinement on each right-hand side
644 *
645  IF ( ref_type .NE. 0 ) THEN
646 
647  prec_type = ilaprec( 'E' )
648 
649  IF ( notran ) THEN
650  CALL zla_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, rwork, work(n+1),
654  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
655  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
656  $ info )
657  ELSE
658  CALL zla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
659  $ nrhs, ab, ldab, afb, ldafb, ipiv, rowequ, r, b,
660  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
661  $ err_bnds_comp, work, rwork, work(n+1),
662  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
663  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
664  $ info )
665  END IF
666  END IF
667 
668  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
669  IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 1) THEN
670 *
671 * Compute scaled normwise condition number cond(A*C).
672 *
673  IF ( colequ .AND. notran ) THEN
674  rcond_tmp = zla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
675  $ ldafb, ipiv, c, .true., info, work, rwork )
676  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
677  rcond_tmp = zla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
678  $ ldafb, ipiv, r, .true., info, work, rwork )
679  ELSE
680  rcond_tmp = zla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
681  $ ldafb, ipiv, c, .false., info, work, rwork )
682  END IF
683  DO j = 1, nrhs
684 *
685 * Cap the error at 1.0.
686 *
687  IF ( n_err_bnds .GE. la_linrx_err_i
688  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0)
689  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
690 *
691 * Threshold the error (see LAWN).
692 *
693  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
694  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
695  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
696  IF ( info .LE. n ) info = n + j
697  ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
698  $ THEN
699  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
700  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
701  END IF
702 *
703 * Save the condition number.
704 *
705  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
706  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
707  END IF
708 
709  END DO
710  END IF
711 
712  IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 2) THEN
713 *
714 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
715 * each right-hand side using the current solution as an estimate of
716 * the true solution. If the componentwise error estimate is too
717 * large, then the solution is a lousy estimate of truth and the
718 * estimated RCOND may be too optimistic. To avoid misleading users,
719 * the inverse condition number is set to 0.0 when the estimated
720 * cwise error is at least CWISE_WRONG.
721 *
722  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
723  DO j = 1, nrhs
724  IF (err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
725  $ THEN
726  rcond_tmp = zla_gbrcond_x( trans, n, kl, ku, ab, ldab,
727  $ afb, ldafb, ipiv, x( 1, j ), info, work, rwork )
728  ELSE
729  rcond_tmp = 0.0d+0
730  END IF
731 *
732 * Cap the error at 1.0.
733 *
734  IF ( n_err_bnds .GE. la_linrx_err_i
735  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
736  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
737 *
738 * Threshold the error (see LAWN).
739 *
740  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
741  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
742  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
743  IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
744  $ .AND. info.LT.n + j ) info = n + j
745  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
746  $ .LT. err_lbnd ) THEN
747  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
748  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
749  END IF
750 *
751 * Save the condition number.
752 *
753  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
754  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
755  END IF
756 
757  END DO
758  END IF
759 *
760  return
761 *
762 * End of ZGBRFSX
763 *
764  END