Uses extra precise iterative refinement to improve the solution to the system of linear equations with a banded matrix A and provides error bounds and backward error estimates.
Syntax
FORTRAN 77:
call sgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info )
call dgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info )
call cgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, rwork, info )
call zgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, rwork, info )
C:
lapack_int LAPACKE_sgbrfsx( int matrix_order, char trans, char equed, lapack_int n, lapack_int kl, lapack_int ku, lapack_int nrhs, const float* ab, lapack_int ldab, const float* afb, lapack_int ldafb, const lapack_int* ipiv, const float* r, const float* c, const float* b, lapack_int ldb, float* x, lapack_int ldx, float* rcond, float* berr, lapack_int n_err_bnds, float* err_bnds_norm, float* err_bnds_comp, lapack_int nparams, float* params );
lapack_int LAPACKE_dgbrfsx( int matrix_order, char trans, char equed, lapack_int n, lapack_int kl, lapack_int ku, lapack_int nrhs, const double* ab, lapack_int ldab, const double* afb, lapack_int ldafb, const lapack_int* ipiv, const double* r, const double* c, const double* b, lapack_int ldb, double* x, lapack_int ldx, double* rcond, double* berr, lapack_int n_err_bnds, double* err_bnds_norm, double* err_bnds_comp, lapack_int nparams, double* params );
lapack_int LAPACKE_cgbrfsx( int matrix_order, char trans, char equed, lapack_int n, lapack_int kl, lapack_int ku, lapack_int nrhs, const lapack_complex_float* ab, lapack_int ldab, const lapack_complex_float* afb, lapack_int ldafb, const lapack_int* ipiv, const float* r, const float* c, const lapack_complex_float* b, lapack_int ldb, lapack_complex_float* x, lapack_int ldx, float* rcond, float* berr, lapack_int n_err_bnds, float* err_bnds_norm, float* err_bnds_comp, lapack_int nparams, float* params );
lapack_int LAPACKE_zgbrfsx( int matrix_order, char trans, char equed, lapack_int n, lapack_int kl, lapack_int ku, lapack_int nrhs, const lapack_complex_double* ab, lapack_int ldab, const lapack_complex_double* afb, lapack_int ldafb, const lapack_int* ipiv, const double* r, const double* c, const lapack_complex_double* b, lapack_int ldb, lapack_complex_double* x, lapack_int ldx, double* rcond, double* berr, lapack_int n_err_bnds, double* err_bnds_norm, double* err_bnds_comp, lapack_int nparams, double* params );
Description
The routine improves the computed solution to a system of linear equations and provides error bounds and backward error estimates for the solution. In addition to a normwise error bound, the code provides a maximum componentwise error bound, if possible. See comments for err_bnds_norm and err_bnds_comp for details of the error bounds.
The original system of linear equations may have been equilibrated before calling this routine, as described by the parameters equed, r, and c below. In this case, the solution and error bounds returned are for the original unequilibrated system.
Input Parameters
The data types are given for the Fortran interface. A <datatype> placeholder, if present, is used for the C interface data types in the C interface section above. See C Interface Conventions for the C interface principal conventions and type definitions.
- trans
CHARACTER*1. Must be 'N', 'T', or 'C'.
Specifies the form of the system of equations:
If trans = 'N', the system has the form A*X = B (No transpose).
If trans = 'T', the system has the form A**T*X = B (Transpose).
If trans = 'C', the system has the form A**H*X = B (Conjugate transpose = Transpose).
- equed
CHARACTER*1. Must be 'N', 'R', 'C', or 'B'.
Specifies the form of equilibration that was done to A before calling this routine.
If equed = 'N', no equilibration was done.
If equed = 'R', row equilibration was done, that is, A has been premultiplied by diag(r).
If equed = 'C', column equilibration was done, that is, A has been postmultiplied by diag(c).
If equed = 'B', both row and column equilibration was done, that is, A has been replaced by diag(r)*A*diag(c). The right-hand side B has been changed accordingly.
- n
INTEGER. The number of linear equations; the order of the matrix A; n ≥ 0.
- kl
INTEGER. The number of subdiagonals within the band of A; kl ≥ 0.
- ku
INTEGER. The number of superdiagonals within the band of A; ku ≥ 0.
- nrhs
INTEGER. The number of right-hand sides; the number of columns of the matrices B and X; nrhs ≥ 0.
- ab, afb, b, work
REAL for sgbrfsx
DOUBLE PRECISION for dgbrfsx
COMPLEX for cgbrfsx
DOUBLE COMPLEX for zgbrfsx.
Arrays: ab(ldab,*), afb(ldafb,*), b(ldb,*), work(*).
The array ab contains the original matrix A in band storage, in rows 1 to kl+ku+1. The j-th column of A is stored in the j-th column of the array ab as follows:
ab(ku+1+i-j,j) = A(i,j) for max(1,j-ku) ≤ i ≤ min(n,j+kl).
The array afb contains details of the LU factorization of the banded matrix A as computed by ?gbtrf. U is stored as an upper triangular banded matrix with kl + ku superdiagonals in rows 1 to kl + ku + 1. The multipliers used during the factorization are stored in rows kl + ku + 2 to 2*kl + ku + 1.
The array b contains the matrix B whose columns are the right-hand sides for the systems of equations. The second dimension of b must be at least max(1,nrhs).
work(*) is a workspace array. The dimension of work must be at least max(1,4*n) for real flavors, and at least max(1,2*n) for complex flavors.
- ldab
INTEGER. The leading dimension of the array ab; ldab ≥ kl+ku+1.
- ldafb
INTEGER. The leading dimension of the array afb; ldafb ≥ 2*kl+ku+1.
- ipiv
INTEGER.
Array, DIMENSION at least max(1, n). Contains the pivot indices as computed by ?gbtrf; for row 1 ≤ i ≤ n, row i of the matrix was interchanged with row ipiv(i).
- r, c
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Arrays: r(n), c(n). The array r contains the row scale factors for A, and the array c contains the column 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 equed = 'R' or 'B', each element of r must be positive.
If equed = 'C' or 'B', A is multiplied on the right by diag(c); if equed = 'N' or 'R', c is not accessed.
If equed = 'C' or 'B', each element of c must be positive.
Each element of r or c should be a power of the radix to ensure a reliable solution and error estimates. Scaling by powers of the radix does not cause rounding errors unless the result underflows or overflows. Rounding errors during scaling lead to refining with a matrix that is not equivalent to the input matrix, producing error estimates that may not be reliable.
- ldb
INTEGER. The leading dimension of the array b; ldb ≥ max(1, n).
- x
REAL for sgbrfsx
DOUBLE PRECISION for dgbrfsx
COMPLEX for cgbrfsx
DOUBLE COMPLEX for zgbrfsx.
Array, DIMENSION (ldx,*).
The solution matrix X as computed by sgbtrs/dgbtrs for real flavors or cgbtrs/zgbtrs for complex flavors.
- ldx
INTEGER. The leading dimension of the output array x; ldx ≥ max(1, n).
- n_err_bnds
INTEGER. Number of error bounds to return for each right-hand side and each type (normwise or componentwise). See err_bnds_norm and err_bnds_comp descriptions in Output Arguments section below.
- nparams
INTEGER. Specifies the number of parameters set in params. If ≤ 0, the params array is never referenced and default values are used.
- params
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Array, DIMENSION nparams. Specifies algorithm parameters. If an entry is less than 0.0, that entry will be filled with the default value used for that parameter. Only positions up to nparams are accessed; defaults are used for higher-numbered parameters. If defaults are acceptable, you can pass nparams = 0, which prevents the source code from accessing the params argument.
params(la_linrx_itref_i = 1) : Whether to perform iterative refinement or not. Default: 1.0 (for single precision flavors), 1.0D+0 (for double precision flavors).
=0.0
No refinement is performed and no error bounds are computed.
=1.0
Use the double-precision refinement algorithm, possibly with doubled-single computations if the compilation environment does not support DOUBLE PRECISION.
(Other values are reserved for futute use.)
params(la_linrx_ithresh_i = 2) : Maximum number of resudual computations allowed for refinement.
Default
10
Aggressive
Set to 100 to permit convergence using approximate factorizations or factorizations other than LU. If the factorization uses a technique other than Gaussian elimination, the quarantees in err_bnds_norm and err_bnds_comp may no longer be trustworthy.
params(la_linrx_cwise_i = 3) : Flag determining if the code will attempt to find a solution with a small componentwise relative error in the double-precision algorithm. Positive is true, 0.0 is false. Default: 1.0 (attempt componentwise convergence).
- iwork
INTEGER. Workspace array, DIMENSION at least max(1, n); used in real flavors only.
- rwork
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Workspace array, DIMENSION at least max(1, 3*n); used in complex flavors only.
Output Parameters
- x
REAL for sgbrfsx
DOUBLE PRECISION for dgbrfsx
COMPLEX for cgbrfsx
DOUBLE COMPLEX for zgbrfsx.
The improved solution matrix X.
- rcond
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Reciprocal scaled condition number. An estimate of the reciprocal Skeel condition number of the matrix A after equilibration (if done). If rcond is less than the machine precision, in particular, if rcond = 0, the matrix is singular to working precision. Note that the error may still be small even if this number is very small and the matrix appears ill-conditioned.
- berr
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Array, DIMENSION at least max(1, nrhs). Contains the componentwise relative backward error for each solution vector x(j), that is, the smallest relative change in any element of A or B that makes x(j) an exact solution.
- err_bnds_norm
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Array, DIMENSION (nrhs,n_err_bnds). For each right-hand side, contains information about various error bounds and condition numbers corresponding to the normwise relative error, which is defined as follows:
Normwise relative error in the i-th solution vector
The array is indexed by the type of error information as described below. There are currently up to three pieces of information returned.
The first index in err_bnds_norm(i,:) corresponds to the i-th right-hand side.
The second index in err_bnds_norm(:,err) contains the following three fields:
err=1
"Trust/don't trust" boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors.
err=2
"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 sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors. This error bound should only be trusted if the previous boolean is true.
err=3
Reciprocal condition number. Estimated normwise reciprocal condition number. Compared with the threshold sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors to determine if the error estimate is "guaranteed". These reciprocal condition numbers are 1/(norm(1/z,inf)*norm(z,inf)) for some appropriately scaled matrix Z.
Let z=s*a, where s scales each row by a power of the radix so all absolute row sums of z are approximately 1.
- err_bnds_comp
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Array, DIMENSION (nrhs,n_err_bnds). For each right-hand side, contains information about various error bounds and condition numbers corresponding to the componentwise relative error, which is defined as follows:
Componentwise relative error in the i-th solution vector:
The array is indexed by the right-hand side i, on which the componentwise relative error depends, and by the type of error information as described below. There are currently up to three pieces of information returned for each right-hand side. If componentwise accuracy is nit requested (params(3) = 0.0), then err_bnds_comp is not accessed. If n_err_bnds < 3, then at most the first (:,n_err_bnds) entries are returned.
The first index in err_bnds_comp(i,:) corresponds to the i-th right-hand side.
The second index in err_bnds_comp(:,err) contains the follwoing three fields:
err=1
"Trust/don't trust" boolean. Trust the answer if the reciprocal condition number is less than the threshold sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors.
err=2
"Guaranteed" error bpound. 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 sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors. This error bound should only be trusted if the previous boolean is true.
err=3
Reciprocal condition number. Estimated componentwise reciprocal condition number. Compared with the threshold sqrt(n)*slamch(ε) for single precision flavors and sqrt(n)*dlamch(ε) for double precision flavors to determine if the error estimate is "guaranteed". These reciprocal condition numbers are 1/(norm(1/z,inf)*norm(z,inf)) for some appropriately scaled matrix Z.
Let z=s*(a*diag(x)), where x is the solution for the current right-hand side and s scales each row of a*diag(x) by a power of the radix so all absolute row sums of z are approximately 1.
- params
REAL for single precision flavors
DOUBLE PRECISION for double precision flavors.
Output parameter only if the input contains erroneous values, namely, in params(1), params(2), params(3). In such a case, the corresponding elements of params are filled with default values on output.
- info
INTEGER. If info = 0, the execution is successful. The solution to every right-hand side is guaranteed.
If info = -i, the i-th parameter had an illegal value.
If 0 < info ≤ 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.
If info = n+j: The solution corresponding to the j-th right-hand side is not guaranteed. The solutions corresponding to other right-hand sides k with k > j may not be guaranteed as well, but only the first such right-hand side is reported. If a small componentwise error is not requested params(3) = 0.0, then the j-th right-hand side is the first with a normwise error bound that is not guaranteed (the smallest j such that err_bnds_norm(j,1) = 0.0 or err_bnds_comp(j,1) = 0.0. See the definition of err_bnds_norm(;,1) and err_bnds_comp(;,1). To get information about all of the right-hand sides, check err_bnds_norm or err_bnds_comp.