Contents
dgbsvx - use the LU factorization to compute the solution to
a real system of linear equations A * X = B, A**T * X = B,
or A**H * X = B,
SUBROUTINE DGBSVX(FACT, TRANSA, N, KL, KU, NRHS, A, LDA, AF,
LDAF, IPIVOT, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR,
BERR, WORK, WORK2, INFO)
CHARACTER * 1 FACT, TRANSA, EQUED
INTEGER N, KL, KU, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER IPIVOT(*), WORK2(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION A(LDA,*), AF(LDAF,*), R(*), C(*), B(LDB,*),
X(LDX,*), FERR(*), BERR(*), WORK(*)
SUBROUTINE DGBSVX_64(FACT, TRANSA, N, KL, KU, NRHS, A, LDA, AF,
LDAF, IPIVOT, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR,
BERR, WORK, WORK2, INFO)
CHARACTER * 1 FACT, TRANSA, EQUED
INTEGER*8 N, KL, KU, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER*8 IPIVOT(*), WORK2(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION A(LDA,*), AF(LDAF,*), R(*), C(*), B(LDB,*),
X(LDX,*), FERR(*), BERR(*), WORK(*)
F95 INTERFACE
SUBROUTINE GBSVX(FACT, [TRANSA], [N], KL, KU, [NRHS], A, [LDA],
AF, [LDAF], IPIVOT, EQUED, R, C, B, [LDB], X, [LDX],
RCOND, FERR, BERR, [WORK], [WORK2], [INFO])
CHARACTER(LEN=1) :: FACT, TRANSA, EQUED
INTEGER :: N, KL, KU, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER, DIMENSION(:) :: IPIVOT, WORK2
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: R, C, FERR, BERR, WORK
REAL(8), DIMENSION(:,:) :: A, AF, B, X
SUBROUTINE GBSVX_64(FACT, [TRANSA], [N], KL, KU, [NRHS], A,
[LDA], AF, [LDAF], IPIVOT, EQUED, R, C, B, [LDB], X, [LDX],
RCOND, FERR, BERR, [WORK], [WORK2], [INFO])
CHARACTER(LEN=1) :: FACT, TRANSA, EQUED
INTEGER(8) :: N, KL, KU, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER(8), DIMENSION(:) :: IPIVOT, WORK2
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: R, C, FERR, BERR, WORK
REAL(8), DIMENSION(:,:) :: A, AF, B, X
C INTERFACE
#include <sunperf.h>
void dgbsvx(char fact, char transa, int n, int kl, int ku,
int nrhs, double *a, int lda, double *af, int
ldaf, int *ipivot, char equed, double *r, double
*c, double *b, int ldb, double *x, int ldx, double
*rcond, double *ferr, double *berr, int *info);
void dgbsvx_64(char fact, char transa, long n, long kl, long
ku, long nrhs, double *a, long lda, double *af,
long ldaf, long *ipivot, char equed, double *r,
double *c, double *b, long ldb, double *x, long
ldx, double *rcond, double *ferr, double *berr,
long *info);
dgbsvx uses the LU factorization to compute the solution to
a real system of linear equations A * X = B, A**T * X = B,
or A**H * X = B, where A is a band matrix of order N with KL
subdiagonals and KU superdiagonals, and X and B are N-by-
NRHS matrices.
Error bounds on the solution and a condition estimate are
also provided.
The following steps are performed by this subroutine:
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 = L * U,
where L is a product of permutation and unit lower tri-
angular
matrices with KL subdiagonals, and U is upper triangular
with
KL+KU superdiagonals.
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 fac-
tored form
of A.
5. Iterative refinement is applied to improve the computed
solution
matrix and calculate error bounds and backward error
estimates
for it.
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.
FACT (input)
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 IPIVOT contain
the factored form of A. If EQUED is not 'N', the
matrix A has been equilibrated with scaling fac-
tors given by R and C. A, AF, and IPIVOT 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.
TRANSA (input)
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)
TRANSA is defaulted to 'N' for F95 INTERFACE.
N (input) The number of linear equations, i.e., the order of
the matrix A. N >= 0.
KL (input)
The number of subdiagonals within the band of A.
KL >= 0.
KU (input)
The number of superdiagonals within the band of A.
KU >= 0.
NRHS (input)
The number of right hand sides, i.e., the number
of columns of the matrices B and X. NRHS >= 0.
A (input/output)
On entry, the 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 A as follows: A(KU+1+i-
j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
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 fol-
lows: EQUED = 'R': A := diag(R) * A
EQUED = 'C': A := A * diag(C)
EQUED = 'B': A := diag(R) * A * diag(C).
LDA (input)
The leading dimension of the array A. LDA >=
KL+KU+1.
AF (input or output)
If FACT = 'F', then AF is an input argument and on
entry contains details of the LU factorization of
the band matrix A, as computed by SGBTRF. U is
stored as an upper triangular band matrix with
KL+KU superdiagonals in rows 1 to KL+KU+1, and the
multipliers used during the factorization are
stored in rows KL+KU+2 to 2*KL+KU+1. 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 details of the LU factorization of
A.
If FACT = 'E', then AF is an output argument and
on exit returns details of the LU factorization of
the equilibrated matrix A (see the description of
A for the form of the equilibrated matrix).
LDAF (input)
The leading dimension of the array AF. LDAF >=
2*KL+KU+1.
IPIVOT (input)
If FACT = 'F', then IPIVOT is an input argument
and on entry contains the pivot indices from the
factorization A = L*U as computed by SGBTRF; row i
of the matrix was interchanged with row IPIVOT(i).
If FACT = 'N', then IPIVOT is an output argument
and on exit contains the pivot indices from the
factorization A = L*U of the original matrix A.
If FACT = 'E', then IPIVOT is an output argument
and on exit contains the pivot indices from the
factorization A = L*U of the equilibrated matrix
A.
EQUED (input or output)
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 equili-
bration, i.e., A has been postmultiplied by
diag(C). = 'B': Both row and column equilibra-
tion, 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)
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.
C (input or output)
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.
B (input/output)
On entry, the right hand side matrix B. On exit,
if EQUED = 'N', B is not modified; if TRANSA = 'N'
and EQUED = 'R' or 'B', B is overwritten by
diag(R)*B; if TRANSA = 'T' or 'C' and EQUED = 'C'
or 'B', B is overwritten by diag(C)*B.
LDB (input)
The leading dimension of the array B. LDB >=
max(1,N).
X (output)
If INFO = 0 or INFO = N+1, the N-by-NRHS 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 TRANSA = 'N' and EQUED
= 'C' or 'B', or inv(diag(R))*X if TRANSA = 'T' or
'C' and EQUED = 'R' or 'B'.
LDX (input)
The leading dimension of the array X. LDX >=
max(1,N).
RCOND (output)
The estimate of the reciprocal condition number of
the matrix A after equilibration (if done). If
RCOND is less than the machine precision (in par-
ticular, if RCOND = 0), the matrix is singular to
working precision. This condition is indicated by
a return code of INFO > 0.
FERR (output)
The estimated forward error bound for each solu-
tion vector X(j) (the j-th column of the solution
matrix X). If XTRUE is the true solution
corresponding to X(j), FERR(j) is an estimated
upper bound for the magnitude of the largest ele-
ment in (X(j) - XTRUE) divided by the magnitude of
the largest element in X(j). The estimate is as
reliable as the estimate for RCOND, and is almost
always a slight overestimate of the true error.
BERR (output)
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).
WORK (workspace)
dimension(3*N) On exit, WORK(1) contains the
reciprocal pivot growth factor norm(A)/norm(U).
The "max absolute element" norm is used. If
WORK(1) 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, condition estimator RCOND, and forward error
bound FERR could be unreliable. If factorization
fails with 0<INFO<=N, then WORK(1) contains the
reciprocal pivot growth factor for the leading
INFO columns of A.
WORK2 (workspace)
dimension(N)
INFO (output)
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an ille-
gal value
> 0: if INFO = i, and i is
<= N: U(i,i) 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.