NAG Library Routine Document
F07HBF (DPBSVX)
1 Purpose
F07HBF (DPBSVX) uses the Cholesky factorization
to compute the solution to a real system of linear equations
where
$A$ is an
$n$ by
$n$ symmetric positive definite band matrix of bandwidth
$\left(2{k}_{d}+1\right)$ and
$X$ and
$B$ are
$n$ by
$r$ matrices. Error bounds on the solution and a condition estimate are also provided.
2 Specification
SUBROUTINE F07HBF ( 
FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO) 
INTEGER 
N, KD, NRHS, LDAB, LDAFB, LDB, LDX, IWORK(N), INFO 
REAL (KIND=nag_wp) 
AB(LDAB,*), AFB(LDAFB,*), S(*), B(LDB,*), X(LDX,*), RCOND, FERR(NRHS), BERR(NRHS), WORK(3*N) 
CHARACTER(1) 
FACT, UPLO, EQUED 

The routine may be called by its
LAPACK
name dpbsvx.
3 Description
F07HBF (DPBSVX) performs the following steps:
 If ${\mathbf{FACT}}=\text{'E'}$, real diagonal scaling factors, ${D}_{S}$, are computed to equilibrate the system:
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 ${D}_{S}A{D}_{S}$ and $B$ by ${D}_{S}B$.
 If ${\mathbf{FACT}}=\text{'N'}$ or $\text{'E'}$, the Cholesky decomposition is used to factor the matrix $A$ (after equilibration if ${\mathbf{FACT}}=\text{'E'}$) as $A={U}^{\mathrm{T}}U$ if ${\mathbf{UPLO}}=\text{'U'}$ or $A=L{L}^{\mathrm{T}}$ if ${\mathbf{UPLO}}=\text{'L'}$, where $U$ is an upper triangular matrix and $L$ is a lower triangular matrix.
 If the leading $i$ by $i$ principal minor of $A$ is not positive definite, then the routine returns with ${\mathbf{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, ${\mathbf{INFO}}={\mathbf{N}+{\mathbf{1}}}$ is returned as a warning, but the routine still goes on to solve for $X$ and compute error bounds as described below.
 The system of equations is solved for $X$ using the factored form of $A$.
 Iterative refinement is applied to improve the computed solution matrix and to calculate error bounds and backward error estimates for it.
 If equilibration was used, the matrix $X$ is premultiplied by ${D}_{S}$ so that it solves the original system before equilibration.
4 References
Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999)
LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia
http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia
5 Parameters
 1: FACT – CHARACTER(1)Input
On entry: specifies whether or not the factorized form of the matrix
$A$ is supplied on entry, and if not, whether the matrix
$A$ should be equilibrated before it is factorized.
 ${\mathbf{FACT}}=\text{'F'}$
 AFB contains the factorized form of $A$. If ${\mathbf{EQUED}}=\text{'Y'}$, the matrix $A$ has been equilibrated with scaling factors given by S. AB and AFB will not be modified.
 ${\mathbf{FACT}}=\text{'N'}$
 The matrix $A$ will be copied to AFB and factorized.
 ${\mathbf{FACT}}=\text{'E'}$
 The matrix $A$ will be equilibrated if necessary, then copied to AFB and factorized.
Constraint:
${\mathbf{FACT}}=\text{'F'}$, $\text{'N'}$ or $\text{'E'}$.
 2: UPLO – CHARACTER(1)Input
On entry: if
${\mathbf{UPLO}}=\text{'U'}$, the upper triangle of
$A$ is stored.
If ${\mathbf{UPLO}}=\text{'L'}$, the lower triangle of $A$ is stored.
Constraint:
${\mathbf{UPLO}}=\text{'U'}$ or $\text{'L'}$.
 3: N – INTEGERInput
On entry: $n$, the number of linear equations, i.e., the order of the matrix $A$.
Constraint:
${\mathbf{N}}\ge 0$.
 4: KD – INTEGERInput
On entry: ${k}_{d}$, the number of superdiagonals of the matrix $A$ if ${\mathbf{UPLO}}=\text{'U'}$, or the number of subdiagonals if ${\mathbf{UPLO}}=\text{'L'}$.
Constraint:
${\mathbf{KD}}\ge 0$.
 5: NRHS – INTEGERInput
On entry: $r$, the number of righthand sides, i.e., the number of columns of the matrix $B$.
Constraint:
${\mathbf{NRHS}}\ge 0$.
 6: AB(LDAB,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
AB
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: the upper or lower triangle of the symmetric band matrix
$A$, except if
${\mathbf{FACT}}=\text{'F'}$ and
${\mathbf{EQUED}}=\text{'Y'}$, in which case
AB must contain the equilibrated matrix
${D}_{S}A{D}_{S}$.
The matrix is stored in rows
$1$ to
${k}_{d}+1$, more precisely,
 if ${\mathbf{UPLO}}=\text{'U'}$, the elements of the upper triangle of $A$ within the band must be stored with element ${A}_{ij}$ in ${\mathbf{AB}}\left({k}_{d}+1+ij,j\right)\text{ for}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j{k}_{d}\right)\le i\le j$;
 if ${\mathbf{UPLO}}=\text{'L'}$, the elements of the lower triangle of $A$ within the band must be stored with element ${A}_{ij}$ in ${\mathbf{AB}}\left(1+ij,j\right)\text{ for}j\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+{k}_{d}\right)\text{.}$
On exit: if
${\mathbf{FACT}}=\text{'E'}$ and
${\mathbf{EQUED}}=\text{'Y'}$,
AB is overwritten by
${D}_{S}A{D}_{S}$.
 7: LDAB – INTEGERInput
On entry: the first dimension of the array
AB as declared in the (sub)program from which F07HBF (DPBSVX) is called.
Constraint:
${\mathbf{LDAB}}\ge {\mathbf{KD}}+1$.
 8: AFB(LDAFB,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
AFB
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'F'}$,
AFB contains the triangular factor
$U$ or
$L$ from the Cholesky factorization
$A={U}^{\mathrm{T}}U$ or
$A=L{L}^{\mathrm{T}}$ of the band matrix
$A$, in the same storage format as
$A$. If
${\mathbf{EQUED}}=\text{'Y'}$,
AFB is the factorized form of the equilibrated matrix
$A$.
On exit: if
${\mathbf{FACT}}=\text{'N'}$,
AFB returns the triangular factor
$U$ or
$L$ from the Cholesky factorization
$A={U}^{\mathrm{T}}U$ or
$A=L{L}^{\mathrm{T}}$.
If
${\mathbf{FACT}}=\text{'E'}$,
AFB returns the triangular factor
$U$ or
$L$ from the Cholesky factorization
$A={U}^{\mathrm{T}}U$ or
$A=L{L}^{\mathrm{T}}$ of the equilibrated matrix
$A$ (see the description of
AB for the form of the equilibrated matrix).
 9: LDAFB – INTEGERInput
On entry: the first dimension of the array
AFB as declared in the (sub)program from which F07HBF (DPBSVX) is called.
Constraint:
${\mathbf{LDAFB}}\ge {\mathbf{KD}}+1$.
 10: EQUED – CHARACTER(1)Input/Output
On entry: if
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
EQUED need not be set.
If
${\mathbf{FACT}}=\text{'F'}$,
EQUED must specify the form of the equilibration that was performed as follows:
 if ${\mathbf{EQUED}}=\text{'N'}$, no equilibration;
 if ${\mathbf{EQUED}}=\text{'Y'}$, equilibration was performed, i.e., $A$ has been replaced by ${D}_{S}A{D}_{S}$.
On exit: if
${\mathbf{FACT}}=\text{'F'}$,
EQUED is unchanged from entry.
Otherwise, if no constraints are violated,
EQUED specifies the form of the equilibration that was performed as specified above.
Constraint:
if ${\mathbf{FACT}}=\text{'F'}$, ${\mathbf{EQUED}}=\text{'N'}$ or $\text{'Y'}$.
 11: S($*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the dimension of the array
S
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
S need not be set.
If
${\mathbf{FACT}}=\text{'F'}$ and
${\mathbf{EQUED}}=\text{'Y'}$,
S must contain the scale factors,
${D}_{S}$, for
$A$; each element of
S must be positive.
On exit: if
${\mathbf{FACT}}=\text{'F'}$,
S is unchanged from entry.
Otherwise, if no constraints are violated and
${\mathbf{EQUED}}=\text{'Y'}$,
S contains the scale factors,
${D}_{S}$, for
$A$; each element of
S is positive.
 12: B(LDB,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
B
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On entry: the $n$ by $r$ righthand side matrix $B$.
On exit: if
${\mathbf{EQUED}}=\text{'N'}$,
B is not modified.
If
${\mathbf{EQUED}}=\text{'Y'}$,
B is overwritten by
${D}_{S}B$.
 13: LDB – INTEGERInput
On entry: the first dimension of the array
B as declared in the (sub)program from which F07HBF (DPBSVX) is called.
Constraint:
${\mathbf{LDB}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 14: X(LDX,$*$) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array
X
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On exit: if ${\mathbf{INFO}}={\mathbf{0}}$ or ${\mathbf{N}+{\mathbf{1}}}$, the $n$ by $r$ solution matrix $X$ to the original system of equations. Note that the arrays $A$ and $B$ are modified on exit if ${\mathbf{EQUED}}=\text{'Y'}$, and the solution to the equilibrated system is ${D}_{S}^{1}X$.
 15: LDX – INTEGERInput
On entry: the first dimension of the array
X as declared in the (sub)program from which F07HBF (DPBSVX) is called.
Constraint:
${\mathbf{LDX}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 16: RCOND – REAL (KIND=nag_wp)Output
On exit: if no constraints are violated, an estimate of the reciprocal condition number of the matrix $A$ (after equilibration if that is performed), computed as ${\mathbf{RCOND}}=1.0/\left({\Vert A\Vert}_{1}{\Vert {A}^{1}\Vert}_{1}\right)$.
 17: FERR(NRHS) – REAL (KIND=nag_wp) arrayOutput
On exit: if
${\mathbf{INFO}}={\mathbf{0}}$ or
${\mathbf{N}+{\mathbf{1}}}$, an estimate of the forward error bound for each computed solution vector, such that
${\Vert {\hat{x}}_{j}{x}_{j}\Vert}_{\infty}/{\Vert {x}_{j}\Vert}_{\infty}\le {\mathbf{FERR}}\left(j\right)$ where
${\hat{x}}_{j}$ is the
$j$th column of the computed solution returned in the array
X and
${x}_{j}$ is the corresponding column of the exact solution
$X$. The estimate is as reliable as the estimate for
RCOND, and is almost always a slight overestimate of the true error.
 18: BERR(NRHS) – REAL (KIND=nag_wp) arrayOutput
On exit: if ${\mathbf{INFO}}={\mathbf{0}}$ or ${\mathbf{N}+{\mathbf{1}}}$, an estimate of the componentwise relative backward error of each computed solution vector ${\hat{x}}_{j}$ (i.e., the smallest relative change in any element of $A$ or $B$ that makes ${\hat{x}}_{j}$ an exact solution).
 19: WORK($3\times {\mathbf{N}}$) – REAL (KIND=nag_wp) arrayWorkspace
 20: IWORK(N) – INTEGER arrayWorkspace
 21: INFO – INTEGEROutput
On exit:
${\mathbf{INFO}}=0$ unless the routine detects an error (see
Section 6).
6 Error Indicators and Warnings
Errors or warnings detected by the routine:
 ${\mathbf{INFO}}<0$
If ${\mathbf{INFO}}=i$, the $i$th argument had an illegal value. An explanatory message is output, and execution of the program is terminated.
 ${\mathbf{INFO}}>0\text{and}{\mathbf{INFO}}\le {\mathbf{N}}$
If ${\mathbf{INFO}}=i$ and $i\le {\mathbf{N}}$, the leading minor of order $i$ of $A$ is not positive definite, so the factorization could not be completed, and the solution has not been computed. ${\mathbf{RCOND}}=0.0$ is returned.
 ${\mathbf{INFO}}={\mathbf{N}}+1$
The triangular matrix
$U$ (or
$L$) 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.
7 Accuracy
For each righthand side vector
$b$, the computed solution
$x$ is the exact solution of a perturbed system of equations
$\left(A+E\right)x=b$, where
 if ${\mathbf{UPLO}}=\text{'U'}$, $\leftE\right\le c\left(n\right)\epsilon \left{U}^{\mathrm{T}}\right\leftU\right$;
 if ${\mathbf{UPLO}}=\text{'L'}$, $\leftE\right\le c\left(n\right)\epsilon \leftL\right\left{L}^{\mathrm{T}}\right$,
$c\left(n\right)$ is a modest linear function of
$n$, and
$\epsilon $ is the
machine precision. See Section 10.1 of
Higham (2002) for further details.
If
$\hat{x}$ is the true solution, then the computed solution
$x$ satisfies a forward error bound of the form
where
$\mathrm{cond}\left(A,\hat{x},b\right)={\Vert \left{A}^{1}\right\left(\leftA\right\left\hat{x}\right+\leftb\right\right)\Vert}_{\infty}/{\Vert \hat{x}\Vert}_{\infty}\le \mathrm{cond}\left(A\right)={\Vert \left{A}^{1}\right\leftA\right\Vert}_{\infty}\le {\kappa}_{\infty}\left(A\right)$.
If
$\hat{x}$ is the
$j$th column of
$X$, then
${w}_{c}$ is returned in
${\mathbf{BERR}}\left(j\right)$ and a bound on
${\Vert x\hat{x}\Vert}_{\infty}/{\Vert \hat{x}\Vert}_{\infty}$ is returned in
${\mathbf{FERR}}\left(j\right)$. See Section 4.4 of
Anderson et al. (1999) for further details.
When $n\gg k$, the factorization of $A$ requires approximately $n{\left(k+1\right)}^{2}$ floating point operations, where $k$ is the number of superdiagonals.
For each righthand side, computation of the backward error involves a minimum of $8nk$ floating point operations. Each step of iterative refinement involves an additional $12nk$ operations. At most five steps of iterative refinement are performed, but usually only one or two steps are required. Estimating the forward error involves solving a number of systems of equations of the form $Ax=b$; the number is usually $4$ or $5$ and never more than $11$. Each solution involves approximately $4nk$ operations.
The complex analogue of this routine is
F07HPF (ZPBSVX).
9 Example
This example solves the equations
where
$A$ is the symmetric positive definite band matrix
and
Error estimates for the solutions, information on equilibration and an estimate of the reciprocal of the condition number of the scaled matrix $A$ are also output.
9.1 Program Text
Program Text (f07hbfe.f90)
9.2 Program Data
Program Data (f07hbfe.d)
9.3 Program Results
Program Results (f07hbfe.r)