Appended are the proposed interfaces for extraprecise
refinement in single precision. Other precisions and
matrix types are similar. The interfaces are, alas,
complicated. Comments appreciated.
SGESVXX : Expert driver that factors, solves, and refines.
SGERFSX : Refinement routine itself.
Still in flux:
1) How to query for the workspace size. The complex
version has integer, real, and complex workspace;
adding three LWORKlike parameters feels a bit
strange.
2) The actual workspace size. If we remove the
cautious scaling from the condition number
estimator, we need either N or 2*N less
workspace.
Jason
 next part 
SUBROUTINE SGESVXX(FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
$ EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS,
$ ERR_BNDS, NPARAMS, PARAMS, WORK, IWORK, INFO)
*
*  LAPACK driver routine (version 3.X) 
* Copyright (c) 20052006, Regents of the University of California,
* Univ. of Tennessee, NAG Ltd., Courant Institute, Argonne National
* Lab, and Rice University
* See COPYING file for license.
*
* .. Scalar Arguments ..
CHARACTER EQUED, FACT, TRANS
INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
$ N_ERR_BNDS
REAL RCOND, RPVGRW
* ..
* .. Array Arguments ..
INTEGER IPIV( * ), IWORK( * )
REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
$ R( * ), C( * ), X( LDX, * ), PARAMS(*),
$ ERR_BNDS(NRHS, 2, N_ERR_BNDS), WORK ( * )
*
* ..
*
* Purpose
* =======
*
* SGESVXX uses the LU factorization to compute the solution to a
* real system of linear equations A * X = B, where A is an
* NbyN matrix and X and B are NbyNRHS matrix.
*
* If iterative refinement is performed, both normwise and maximum
* componentwise error bound is returned, as well as a condition
* estimate. SGESVXX will return a solution with O(eps) error (where
* eps is the working machine precision) unless the matrix is very
* illconditioned. Relevant condition numbers also are calculated
* and returned.
*
* Description
* ===========
*
* The following steps are performed:
*
* 1. If FACT = 'E', real scaling factors are computed to equilibrate
* the system:
*
* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
*
* Whether or not the system will be equilibrated depends on the
* scaling of the matrix A, but if equilibration is used, A is
* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
* or diag(C)*B (if TRANS = 'T' or 'C').
*
* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
* matrix A (after equilibration if FACT = 'E') as
*
* A = P * L * U,
*
* where P is a permutation matrix, L is a unit lower triangular
* matrix, and U is upper triangular.
*
* 3. If some U(i,i)=0, so that U is exactly singular, then the routine
* returns with INFO = i. Otherwise, the factored form of A is used
* to estimate the condition number of the matrix A. If the
* reciprocal of the condition number is less than machine precision,
* INFO = N+1 is returned as a warning, but the routine still goes on
* to solve for X and compute error bounds as described below.
*
* 4. The system of equations is solved for X using the factored form
* of A.
*
* 5. PARAMS(ITREF_PARAM) specifies whether we are performing
* iterative refinement. Refinement calculates the residual to
* at least twice the working precision.
*
* 6. If equilibration was used, the matrix X is premultiplied by
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
* Arguments
* =========
*
* Some parameters are bundled in the PARAMS array. These settings
* affect other behavior, so we describe them first. If the defaults
* are acceptable, users can pass NPARAMS = 0 and the code will never
* access the PARAMS argument.
*
* NPARAMS (input) INTEGER
* Specifies the number of parameters set in PARAMS. If .LE. 0, the
* PARAMS array is never referenced.
*
* PARAMS (input / output) REAL array, dimension NPARAMS
* Specifies algorithm parameters. If an entry is .LT. 0.0, then
* that entry will be filled with default value used for that
* parameter. Only positions up to NPARAMS are accessed; defaults
* are used for unavailable parameters. PARAMETER declarations for
* the index parameters below are specified in sgerfsx_params.fh.
*
* PARAMS(ITREF_PARAM = 1) : Whether to perform iterative
* refinement or not.
* Default: 1.0
* = 0.0 : No refinement is perfomed, and no error bounds are
* computed. Error arguments are not touched.
* = 1.0 : Use the doubleprecision refinement algorithm,
* possibly with doubledsingle computations if the
* compilation environment does not support DOUBLE
* PRECISION.
* (other values are reserved for future use)
*
* PARAMS(RCONDTHRESH_PARAM = 2) : Reciprocal condition number
* threshold where the error estimates are no longer
* considered trustworthy.
* Change with extreme caution.
* Default: SQRT(N) * SLAMCH('Epsilon')
*
* PARAMS(ITHRESH_PARAM = 3) : Total number of residual computations
* allowed for refinement.
* Default: 10 for doubleprecision, 5 for singleprecision
* Aggressive: Set to 100 for our "aggresive" settings.
*
* PARAMS(COMPONENTWISE_PARAM = 4) : Flag determining if the code
* will attempt to find a solution with small componentwise
* relative error in the doubleprecision algorithm. Nonzero
* is true, 0.0 is false.
* Default: 1.0 (attempt componentwise convergenge)
*
* PARAMS(RTHRESH_PARAM = 5) : Ratio of consecutive "step sizes"
* required to continue relative normwise or componentwise
* refinement. If the norm of the previous step (dx or dz) is
* SPREV and the norm of the current step is SCURR, then SCURR
* / SPREV must be .LE. PARAMS(RTHRESH_PARAM) to continue.
* Default: 0.5
* Aggressive: Set to 0.9 for our "aggresive" settings.
*
* PARAMS(DZTHRESH_PARAM = 6) : Threshold where the solution is
* considered stable enough for computing componentwise
* measurements.
* Default: 0.25
*
* FACT (input) CHARACTER*1
* Specifies whether or not the factored form of the matrix A is
* supplied on entry, and if not, whether the matrix A should be
* equilibrated before it is factored.
* = 'F': On entry, AF and IPIV contain the factored form of A.
* If EQUED is not 'N', the matrix A has been
* equilibrated with scaling factors given by R and C.
* A, AF, and IPIV are not modified.
* = 'N': The matrix A will be copied to AF and factored.
* = 'E': The matrix A will be equilibrated if necessary, then
* copied to AF and factored.
*
* TRANS (input) CHARACTER*1
* Specifies the form of the system of equations:
* = 'N': A * X = B (No transpose)
* = 'T': A**T * X = B (Transpose)
* = 'C': A**H * X = B (Transpose)
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input/output) REAL array, dimension (LDA,N)
* On entry, the NbyN matrix A. If FACT = 'F' and EQUED is
* not 'N', then A must have been equilibrated by the scaling
* factors in R and/or C. A is not modified if FACT = 'F' or
* 'N', or if FACT = 'E' and EQUED = 'N' on exit.
*
* On exit, if EQUED .ne. 'N', A is scaled as follows:
* EQUED = 'R': A := diag(R) * A
* EQUED = 'C': A := A * diag(C)
* EQUED = 'B': A := diag(R) * A * diag(C).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input or output) REAL array, dimension (LDAF,N)
* If FACT = 'F', then AF is an input argument and on entry
* contains the factors L and U from the factorization
* A = P*L*U as computed by SGETRF. If EQUED .ne. 'N', then
* AF is the factored form of the equilibrated matrix A.
*
* If FACT = 'N', then AF is an output argument and on exit
* returns the factors L and U from the factorization A = P*L*U
* of the original matrix A.
*
* If FACT = 'E', then AF is an output argument and on exit
* returns the factors L and U from the factorization A = P*L*U
* of the equilibrated matrix A (see the description of A for
* the form of the equilibrated matrix).
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input or output) INTEGER array, dimension (N)
* If FACT = 'F', then IPIV is an input argument and on entry
* contains the pivot indices from the factorization A = P*L*U
* as computed by SGETRF; row i of the matrix was interchanged
* with row IPIV(i).
*
* If FACT = 'N', then IPIV is an output argument and on exit
* contains the pivot indices from the factorization A = P*L*U
* of the original matrix A.
*
* If FACT = 'E', then IPIV is an output argument and on exit
* contains the pivot indices from the factorization A = P*L*U
* of the equilibrated matrix A.
*
* EQUED (input or output) CHARACTER*1
* Specifies the form of equilibration that was done.
* = 'N': No equilibration (always true if FACT = 'N').
* = 'R': Row equilibration, i.e., A has been premultiplied by
* diag(R).
* = 'C': Column equilibration, i.e., A has been postmultiplied
* by diag(C).
* = 'B': Both row and column equilibration, i.e., A has been
* replaced by diag(R) * A * diag(C).
* EQUED is an input argument if FACT = 'F'; otherwise, it is an
* output argument.
*
* R (input or output) REAL array, dimension (N)
* The row scale factors for A. If EQUED = 'R' or 'B', A is
* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
* is not accessed. R is an input argument if FACT = 'F';
* otherwise, R is an output argument. If FACT = 'F' and
* EQUED = 'R' or 'B', each element of R must be positive.
* If R is output, each element of R is a power of the radix.
* If R is input, and more than working precision is used,
* then each element of R should be a power of the radix
* for the solution and error bound to be reliable.
*
* C (input or output) REAL array, dimension (N)
* The column scale factors for A. If EQUED = 'C' or 'B', A is
* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
* is not accessed. C is an input argument if FACT = 'F';
* otherwise, C is an output argument. If FACT = 'F' and
* EQUED = 'C' or 'B', each element of C must be positive.
* If C is output, each element of C is a power of the radix.
* If C is input, and more than working precision is used,
* then each element of C should be a power of the radix
* for the solution and error bound to be reliable.
*
* B (input/output) REAL array, dimension (LDB,NRHS)
* On entry, the NbyNRHS right hand side matrix B.
* On exit,
* if EQUED = 'N', B is not modified;
* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
* diag(R)*B;
* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
* overwritten by diag(C)*B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (output) REAL array, dimension (LDX,NRHS)
* If INFO = 0 or INFO = N+1, the NbyNRHS solution matrix X
* to the original system of equations. Note that A and B are
* modified on exit if EQUED .ne. 'N', and the solution to the
* equilibrated system is inv(diag(C))*X if TRANS = 'N' and
* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
* and EQUED = 'R' or 'B'.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* RCOND (output) REAL
* Reciprocal scaled condition number. This is an estimate of the
* reciprocal condition number of the matrix A after equilibration
* (if done). If this is less than the machine precision (in
* particular, if it is zero), the matrix is singular to working
* precision. This condition is indicated by a return code of
* INFO > 0. Note that the error may still be small even if this
* number is very small and the matrix appears illconditioned.
*
* RPVGRW (output) REAL
* Reciprocal pivot growth. On exit, this contains the reciprocal
* pivot growth factor norm(A)/norm(U). The "max absolute element"
* norm is used. If this is much less than 1, then the stability of
* the LU factorization of the (equilibrated) matrix A could be poor.
* This also means that the solution X, estimated condition numbers,
* and error bounds could be unreliable. If factorization fails with
* 0<INFO<=N, then this contains the reciprocal pivot growth factor
* for the leading INFO columns of A. In SGESVX, this quantity is
* returned in WORK(1).
*
* BERR (output) REAL array, dimension (NRHS)
* Componentwise relative backward error. Otherwise, this is the
* componentwise relative backward error of each solution vector X(j)
* (i.e., the smallest relative change in any element of A or B that
* makes X(j) an exact solution).
*
* N_ERR_BNDS (input) INTEGER
* Number of errorbound entries to return for each right hand side
* and each type (normwise or componentwise). See ERR_BNDS below.
*
* ERR_BNDS (output) REAL array, dimension (NRHS, 2, N_ERR_BNDS)
* The array containing information about various error
* bounds and condition numbers. The array is indexed by the right
* hand side j, the error "norm", and the type of error information
* as described below. There currently are up to three pieces of
* information returned for each righthand side and "norm". If
* componentwise accuracy is disabled in PARAMS, the (:,2,:) entries
* are not accessed. If N_ERR_BNDS .LT. 3, then at most the first
* (:,:,N_ERR_BNDS) entries are accessed.
*
* The first index in ERR_BNDS(j,:,:) corresponds to the jth right
* hand side.
*
* The second index in ERR_BNDS(:,nrm,:) corresponds to the "norm".
* If XTRUE(j) is the true solution for X(j) and we let 0/0 = 0, then
* the norms are
*
* = 1 Normwise relative:
* max_j (abs(XTRUE(j)  X(j)))
* 
* max_j abs(X(j))
*
* = 2 Componentwise relative
* abs(XTRUE(j)  X(j))
* max_j 
* abs(X(j))
*
* The final index in ERR_BNDS(:,:,err) is one of:
* = 1 "Guaranteed" error bound: The estimated forward error,
* almost certainly within a factor of 10 of the true error
* so long as the next entry is greater than the threshold
* PARAMS(RCONDTHRESH_PARAM). If the reciprocal condition
* number is less than that threshold, the estimate is
* considered untrustworthy and is set to at least 1.0.
*
* = 2 Reciprocal condition number: Estimated reciprocal condition
* number scaled for the corresponding norm. Compared with
* PARAMS(RCONDTHRESH_PARAM) to determine if the error
* estimate is "guaranteed". (That parameter defaults to
* sqrt(n) * slamch('Epsilon').) These reciprocal condition
* numbers are 1 / (norm(Z^{1},inf) * norm(Z,inf)) for some
* appropriately scaled matrix Z.
* Normwise: Let Z = S*A, where S scales each row by a power
* of the radix so all row sums of Z are
* approximately 1.
* Componentwise: Let Z = S*(A*diag(x)), where x is the
* solution for the current righthand side and S
* scales each row of A*diag(x) by a power of the
* radix so all row sums of Z are approximately 1.
*
* = 3 Raw error estimate: Error estimated as detailed in LAPACK
* Working Note 165. This estimate will never be greater
* than the first, "guaranteed" entry. General systems
* should not use or trust this error estimate. Only use
* this if there is some known structure to your problem that
* is not considered in our condition number.
*
* See Lapack Working Note 165 for further details and extra
* cautions.
*
* WORK (workspace) REAL array, dimension (5*N)
*
* IWORK (workspace) INTEGER array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = i, the ith argument had an illegal value
* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
* has been completed, but the factor U is exactly singular, so
* the solution and error bounds could not be computed. RCOND = 0
* is returned.
* = N+1: U is nonsingular, but RCOND is less than machine precision,
* meaning that the matrix is singular to working precision.
* Nevertheless, the solution and error bounds are computed because
* there are a number of situations where the computed solution can
* be more accurate than the value of RCOND would suggest. However,
* if INFO=N+1+J for some righthand side J, then INFO is never set
* to N+1.
* = N+1+J: The iterative refinement for the Jth righthand side did
* not produce a small error bound. Other righthand sides K with
* K > J may have not converged as well, but only the first
* righthand side to not converge is reported.
*
* =====================================================================
END
 next part 
SUBROUTINE SGERFSX(TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
$ R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS,
$ ERR_BNDS, NPARAMS, PARAMS, WORK, IWORK, INFO )
*
*  LAPACK routine (version 3.X) 
* Copyright (c) 20052006, Regents of the University of California,
* Univ. of Tennessee, NAG Ltd., Courant Institute, Argonne National
* Lab, and Rice University
* See COPYING file for license.
*
* .. Scalar Arguments ..
CHARACTER TRANS, EQUED
INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
$ N_ERR_BNDS
REAL RCOND, RPVGRW
*
* ..
* .. Array Arguments ..
INTEGER IPIV( * ), IWORK( * )
REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
$ R( * ), C( * ), X( LDX, * ), PARAMS(*),
$ ERR_BNDS(NRHS, 2, N_ERR_BNDS), WORK ( * )
*
* Purpose
* =======
*
* SGERFSX improves the computed solution to a system of linear
* equations and provides error bounds and backward error estimates for
* the solution. In addition to normwise error bound, the code provides
* maximum componentwise error bound if possible. See comments for
* ERR_BNDS for details of the error bounds.
*
* The original system of linear equations may have been equilibrated
* before calling this routine, as described by arguments EQUED, R and C
* below. In this case, the solution and error bounds returned are for
* the original unequilibrated system.
*
* Arguments
* =========
*
* Some parameters are bundled in the PARAMS array. These settings
* affect other behavior, so we describe them first. If the defaults
* are acceptable, users can pass NPARAMS = 0 and the code will never
* access the PARAMS argument.
*
* NPARAMS (input) INTEGER
* Specifies the number of parameters set in PARAMS. If .LE. 0, the
* PARAMS array is never referenced.
*
* PARAMS (input / output) REAL array, dimension NPARAMS
* Specifies algorithm parameters. If an entry is .LT. 0.0, then
* that entry will be filled with default value used for that
* parameter. Only positions up to NPARAMS are accessed; defaults
* are used for unavailable parameters. PARAMETER declarations for
* the index parameters below are specified in sgerfsx_params.fh.
*
* PARAMS(ITREF_PARAM = 1) : Whether to perform iterative
* refinement or not.
* Default: 1.0
* = 0.0 : No refinement is perfomed, and no error bounds are
* computed. Error arguments are not touched.
* = 1.0 : Use the doubleprecision refinement algorithm,
* possibly with doubledsingle computations if the
* compilation environment does not support DOUBLE
* PRECISION.
* (other values are reserved for future use)
*
* PARAMS(RCONDTHRESH_PARAM = 2) : Reciprocal condition number
* threshold where the error estimates are no longer
* considered trustworthy.
* Change with extreme caution.
* Default: SQRT(N) * SLAMCH('Epsilon')
*
* PARAMS(ITHRESH_PARAM = 3) : Total number of residual computations
* allowed for refinement.
* Default: 10 for doubleprecision, 5 for singleprecision
* Aggressive: Set to 100 for our "aggresive" settings.
*
* PARAMS(COMPONENTWISE_PARAM = 4) : Flag determining if the code
* will attempt to find a solution with small componentwise
* relative error in the doubleprecision algorithm. Positive
* is true, 0.0 is false.
* Default: 1.0 (attempt componentwise convergenge)
*
* PARAMS(RTHRESH_PARAM = 5) : Ratio of consecutive "step sizes"
* required to continue relative normwise or componentwise
* refinement. If the norm of the previous step (dx or dz) is
* SPREV and the norm of the current step is SCURR, then SCURR
* / SPREV must be .LE. PARAMS(RTHRESH_PARAM) to continue.
* Default: 0.5
* Aggressive: Set to 0.9 for our "aggresive" settings.
*
* PARAMS(DZTHRESH_PARAM = 6) : Threshold where the solution is
* considered stable enough for computing componentwise
* measurements.
* Default: 0.25
*
* TRANS (input) CHARACTER*1
* Specifies the form of the system of equations:
* = 'N': A * X = B (No transpose)
* = 'T': A**T * X = B (Transpose)
* = 'C': A**H * X = B (Conjugate transpose = Transpose)
*
* EQUED (input) CHARACTER*1
* Specifies the form of equilibration that was done to A
* before calling this routine. This is needed to compute
* the solution and error bounds correctly.
* = 'N': No equilibration
* = 'R': Row equilibration, i.e., A has been premultiplied by diag(R).
* = 'C': Column equilibration, i.e., A has been postmultiplied
* by diag(C).
* = 'B': Both row and column equilibration, i.e., A has been
* replaced by diag(R) * A * diag(C).
* The right hand side B has been changed accordingly.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input) REAL array, dimension (LDA,N)
* The original NbyN matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input) REAL array, dimension (LDAF,N)
* The factors L and U from the factorization A = P*L*U
* as computed by SGETRF.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* The pivot indices from SGETRF; for 1<=i<=N, row i of the
* matrix was interchanged with row IPIV(i).
*
* R (input) REAL array, dimension (N)
* The row scale factors for A. If EQUED = 'R' or 'B', A is
* multiplied on the left by diag(R). If EQUED = 'N' or 'C',
* R is not accessed. If R is input, and more than working
* precision is used, then each element of R should be a power
* of the radix for the solution and error bound to be reliable.
*
* C (input) REAL array, dimension (N)
* The column scale factors for A. If EQUED = 'C' or 'B', A is
* multiplied on the right by diag(C). If EQUED = 'N' or 'R',
* C is not accessed. If C is input, and more than working
* precision is used, then each element of C should be a power
* of the radix for the solution and error bound to be reliable.
*
* B (input) REAL array, dimension (LDB,NRHS)
* The right hand side matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (input/output) REAL array, dimension (LDX,NRHS)
* On entry, the solution matrix X, as computed by SGETRS.
* On exit, the improved solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* RCOND (output) REAL
* Reciprocal scaled condition number. This is an estimate of the
* reciprocal condition number of the matrix A after equilibration
* (if done). If this is less than the machine precision (in
* particular, if it is zero), the matrix is singular to working
* precision. This condition is indicated by a return code of
* INFO > 0. Note that the error may still be small even if this
* number is very small and the matrix appears illconditioned.
*
* RPVGRW (output) REAL
* Reciprocal pivot growth. On exit, this contains the reciprocal
* pivot growth factor norm(A)/norm(U). The "max absolute element"
* norm is used. If this is much less than 1, then the stability of
* the LU factorization of the (equilibrated) matrix A could be poor.
* This also means that the solution X, estimated condition numbers,
* and error bounds could be unreliable. If factorization fails with
* 0<INFO<=N, then this contains the reciprocal pivot growth factor
* for the leading INFO columns of A. In SGESVX, this quantity is
* returned in WORK(1).
*
* BERR (output) REAL array, dimension (NRHS)
* Componentwise relative backward error. Otherwise, this is the
* componentwise relative backward error of each solution vector X(j)
* (i.e., the smallest relative change in any element of A or B that
* makes X(j) an exact solution).
*
* N_ERR_BNDS (input) INTEGER
* Number of errorbound entries to return for each right hand side
* and each type (normwise or componentwise). See ERR_BNDS below.
*
* ERR_BNDS (output) REAL array, dimension (NRHS, 2, N_ERR_BNDS)
* The array containing information about various error
* bounds and condition numbers. The array is indexed by the right
* hand side j, the error "norm", and the type of error information
* as described below. There currently are up to three pieces of
* information returned for each righthand side and "norm". If
* componentwise accuracy is disabled in PARAMS, the (:,2,:) entries
* are not accessed. If N_ERR_BNDS .LT. 3, then at most the first
* (:,:,N_ERR_BNDS) entries are accessed.
*
* The first index in ERR_BNDS(j,:,:) corresponds to the jth right
* hand side.
*
* The second index in ERR_BNDS(:,nrm,:) corresponds to the "norm".
* If XTRUE(j) is the true solution for X(j) and we let 0/0 = 0, then
* the norms are
*
* = 1 Normwise relative:
* max_j (abs(XTRUE(j)  X(j)))
* 
* max_j abs(X(j))
*
* = 2 Componentwise relative
* abs(XTRUE(j)  X(j))
* max_j 
* abs(X(j))
*
* The final index in ERR_BNDS(:,:,err) is one of:
* = 1 "Guaranteed" error bound: The estimated forward error,
* almost certainly within a factor of 10 of the true error
* so long as the next entry is greater than the threshold
* PARAMS(RCONDTHRESH_PARAM). If the reciprocal condition
* number is less than that threshold, the estimate is
* considered untrustworthy and is set to at least 1.0.
*
* = 2 Reciprocal condition number: Estimated reciprocal condition
* number scaled for the corresponding norm. Compared with
* PARAMS(RCONDTHRESH_PARAM) to determine if the error
* estimate is "guaranteed". (That parameter defaults to
* sqrt(n) * slamch('Epsilon').) These reciprocal condition
* numbers are 1 / (norm(Z^{1},inf) * norm(Z,inf)) for some
* appropriately scaled matrix Z.
* Normwise: Let Z = S*A, where S scales each row by a power
* of the radix so all row sums of Z are
* approximately 1.
* Componentwise: Let Z = S*(A*diag(x)), where x is the
* solution for the current righthand side and S
* scales each row of A*diag(x) by a power of the
* radix so all row sums of Z are approximately 1.
*
* = 3 Raw error estimate: Error estimated as detailed in LAPACK
* Working Note 165. This estimate will never be greater
* than the first, "guaranteed" entry. General systems
* should not use or trust this error estimate. Only use
* this if there is some known structure to your problem that
* is not considered in our condition number.
*
* See Lapack Working Note 165 for further details and extra
* cautions.
*
* WORK (workspace) REAL array, dimension (3*N)
*
* IWORK (workspace) INTEGER array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = i, the ith argument had an illegal value
* = N+1+J: The iterative refinement for the Jth righthand side
* did not produce small normwise error bound. Other righthand
* sides K with K > J may have not converged as well, but only
* the first righthand side to not converge is reported.
*
* For further details see LAPACK Working Note 165.
*
* =====================================================================
END
