3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZGBRFSX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbrfsx.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbrfsx.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbrfsx.f">
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,
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
33 * .. Array Arguments ..
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( * )
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
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.
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.
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)
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
90 *> = 'C': Column equilibration, i.e., A has been postmultiplied
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.
100 *> The order of the matrix A. N >= 0.
106 *> The number of subdiagonals within the band of A. KL >= 0.
112 *> The number of superdiagonals within the band of A. KU >= 0.
118 *> The number of right hand sides, i.e., the number of columns
119 *> of the matrices B and X. NRHS >= 0.
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).
134 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
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.
150 *> The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
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).
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
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
198 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
199 *> The right hand side matrix B.
205 *> The leading dimension of the array B. LDB >= max(1,N).
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.
218 *> The leading dimension of the array X. LDX >= max(1,N).
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-
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).
242 *> \param[in] N_ERR_BNDS
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.
250 *> \param[out] ERR_BNDS_NORM
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:
257 *> Normwise relative error in the ith solution vector:
258 *> max_j (abs(XTRUE(j,i) - X(j,i)))
259 *> ------------------------------
262 *> The array is indexed by the type of error information as described
263 *> below. There currently are up to three pieces of information
266 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
269 *> The second index in ERR_BNDS_NORM(:,err) contains the following
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').
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.
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.
290 *> See Lapack Working Note 165 for further details and extra
294 *> \param[out] ERR_BNDS_COMP
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:
301 *> Componentwise relative error in the ith solution vector:
302 *> abs(XTRUE(j,i) - X(j,i))
303 *> max_j ----------------------
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.
314 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
317 *> The second index in ERR_BNDS_COMP(:,err) contains the following
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').
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.
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.
340 *> See Lapack Working Note 165 for further details and extra
344 *> \param[in] NPARAMS
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.
351 *> \param[in,out] PARAMS
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.
359 *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
360 *> refinement or not.
362 *> = 0.0 : No refinement is performed, and no error bounds are
364 *> = 1.0 : Use the double-precision refinement algorithm,
365 *> possibly with doubled-single computations if the
366 *> compilation environment does not support DOUBLE
368 *> (other values are reserved for future use)
370 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
371 *> computations allowed for refinement.
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
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)
389 *> WORK is COMPLEX*16 array, dimension (2*N)
394 *> RWORK is DOUBLE PRECISION array, dimension (2*N)
400 *> = 0: Successful exit. The solution to every right-hand side is
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
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
427 *> \author Univ. of Tennessee
428 *> \author Univ. of California Berkeley
429 *> \author Univ. of Colorado Denver
434 *> \ingroup complex16GBcomputational
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,
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..--
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
454 * .. Array Arguments ..
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( * )
463 * ==================================================================
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,
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,
483 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
484 PARAMETER ( LA_LINRX_RCOND_I = 3 )
486 * .. Local Scalars ..
488 LOGICAL ROWEQU, COLEQU, NOTRAN, IGNORE_CWISE
489 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS,
491 DOUBLE PRECISION ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND,
492 $ CWISE_WRONG, RTHRESH, UNSTABLE_THRESH
494 * .. External Subroutines ..
495 EXTERNAL XERBLA, ZGBCON, ZLA_GBRFSX_EXTENDED
497 * .. Intrinsic Functions ..
498 INTRINSIC MAX, SQRT, TRANSFER
500 * .. External Functions ..
501 EXTERNAL LSAME, ILAPREC
502 EXTERNAL DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
503 DOUBLE PRECISION DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
505 INTEGER ILATRANS, ILAPREC
507 * .. Executable Statements ..
509 * Check the input parameters.
512 TRANS_TYPE = ILATRANS( TRANS )
513 REF_TYPE = INT( ITREF_DEFAULT )
514 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
515 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
516 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
518 REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
522 * Set default parameters.
524 ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
525 ITHRESH = INT( ITHRESH_DEFAULT )
526 RTHRESH = RTHRESH_DEFAULT
527 UNSTABLE_THRESH = DZTHRESH_DEFAULT
528 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
530 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
531 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
532 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
534 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
537 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
538 IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
539 IF ( IGNORE_CWISE ) THEN
540 PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
542 PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
545 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
548 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
550 ELSE IF ( IGNORE_CWISE ) THEN
556 NOTRAN = LSAME( TRANS, 'N' )
557 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
558 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
560 * Test input parameters.
562 IF( TRANS_TYPE.EQ.-1 ) THEN
564 ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
565 $ .NOT.LSAME( EQUED, 'N' ) ) THEN
567 ELSE IF( N.LT.0 ) THEN
569 ELSE IF( KL.LT.0 ) THEN
571 ELSE IF( KU.LT.0 ) THEN
573 ELSE IF( NRHS.LT.0 ) THEN
575 ELSE IF( LDAB.LT.KL+KU+1 ) THEN
577 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
579 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
581 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
585 CALL XERBLA( 'ZGBRFSX', -INFO )
589 * Quick return if possible.
591 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
595 IF ( N_ERR_BNDS .GE. 1 ) THEN
596 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
597 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
599 IF ( N_ERR_BNDS .GE. 2 ) THEN
600 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
601 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
603 IF ( N_ERR_BNDS .GE. 3 ) THEN
604 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
605 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
611 * Default to failure.
616 IF ( N_ERR_BNDS .GE. 1 ) THEN
617 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
618 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
620 IF ( N_ERR_BNDS .GE. 2 ) THEN
621 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
622 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
624 IF ( N_ERR_BNDS .GE. 3 ) THEN
625 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
626 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
630 * Compute the norm of A and the reciprocal of the condition
638 ANORM = ZLANGB( NORM, N, KL, KU, AB, LDAB, RWORK )
639 CALL ZGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND,
640 $ WORK, RWORK, INFO )
642 * Perform refinement on each right-hand side
644 IF ( REF_TYPE .NE. 0 .AND. INFO .EQ. 0 ) THEN
646 PREC_TYPE = ILAPREC( 'E' )
649 CALL ZLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
650 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B,
651 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
652 $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
653 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
654 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
657 CALL ZLA_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, RWORK, WORK(N+1),
661 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
662 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
667 ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
668 IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN
670 * Compute scaled normwise condition number cond(A*C).
672 IF ( COLEQU .AND. NOTRAN ) THEN
673 RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
674 $ LDAFB, IPIV, C, .TRUE., INFO, WORK, RWORK )
675 ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
676 RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
677 $ LDAFB, IPIV, R, .TRUE., INFO, WORK, RWORK )
679 RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
680 $ LDAFB, IPIV, C, .FALSE., INFO, WORK, RWORK )
684 * Cap the error at 1.0.
686 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
687 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0)
688 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
690 * Threshold the error (see LAWN).
692 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
693 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
694 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
695 IF ( INFO .LE. N ) INFO = N + J
696 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
698 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
699 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
702 * Save the condition number.
704 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
705 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
711 IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
713 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
714 * each right-hand side using the current solution as an estimate of
715 * the true solution. If the componentwise error estimate is too
716 * large, then the solution is a lousy estimate of truth and the
717 * estimated RCOND may be too optimistic. To avoid misleading users,
718 * the inverse condition number is set to 0.0 when the estimated
719 * cwise error is at least CWISE_WRONG.
721 CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
723 IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
725 RCOND_TMP = ZLA_GBRCOND_X( TRANS, N, KL, KU, AB, LDAB,
726 $ AFB, LDAFB, IPIV, X( 1, J ), INFO, WORK, RWORK )
731 * Cap the error at 1.0.
733 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
734 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
735 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
737 * Threshold the error (see LAWN).
739 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
740 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
741 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
742 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0
743 $ .AND. INFO.LT.N + J ) INFO = N + J
744 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
745 $ .LT. ERR_LBND ) THEN
746 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
747 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
750 * Save the condition number.
752 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
753 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP