\documentstyle [11pt] {article}
%
\parindent 0in
\topmargin=-0.60in
\hoffset=-0.75in
\textheight=8.75in \vsize\textheight
\textwidth=6.5in \hsize\textwidth
%
\begin {document}
%
{\small
\begin{center}
{\hfill \bf PHYS 410/555 \hfill SOLVING GENERAL LINEAR BANDED SYSTEMS \hfill October 2001 \hfill}
\end  {center}

\medskip
In the following, we will assume that we are solving a system of $n$
equations in $n$ unknowns so that the number of rows and columns 
in our matrix {\bf A} are equal. 

\medskip
A general banded system is one in which the non-zero elements of {\bf A}
are confined to a certain number of diagonals of the matrix, usually 
clustered about the main diagonal.  Tridiagonal systems clearly 
fall into this category. Systems with more than 3 non-zero diagonals 
also often arise in finite difference calculations (and other 
applications) and can also be efficiently LU-decomposed by taking 
advantage of the non-zero structure of the array. 

\medskip
\medskip
{\bf Example: Pentadiagonal (5-diagonal) system}

\medskip
Recall the BVP we have considered previously 
$$
u''(x) = f(x) \quad\quad
0 \le x \le 1 \quad\quad u(0), u(1)\,\,\,\, \hbox{\rm specified}
$$
Let us approximate $u''(x)$ using a fourth-order ($O(h^4)$) 
finite-difference approximation:

\medskip
$$
\frac{-u_{i-2} + 16 u_{i-1} - 30 u_i + 16 u_{i+1} - u_{i+2}}{12h^2}
 = u''(x_j) + O(h^4) \quad \quad i = 3,\ldots,n-2
$$

\medskip
(Part of Homework 4 involves using {\tt maple} to verify that this 
{\sl is} an $O(h^4)$ approximation to the second derivative.)
Because this approximation couples the $i$th unknown, $u_i$ to 
its nearest, and next-nearest neighbors, $u_{i-1}, u_{i+1}, u_{i-2}$
and $u_{i+2}$, the linear system which results from applying the 
approximation at mesh points $x_i, \quad i = 3,\ldots,n-2$ is 
(ignoring momentarily the equations for $u_1, u_2, u_{n-1}$ and $u_n$)
is clearly pentadiagonal:

\medskip
We define the {\sl total bandwidth}, $w$, of a matrix as the total number 
of diagonals which span the non-zero portion of the matrix. In
our current example, the bandwidth is $w = 5$, and all 5 diagonals 
are generally comprised of non-0 elements.  However, particularly 
for finite-difference approximations to problems in multiple dimensions,
some of the diagonals within the bandwidth of the matrix may have 
elements which are all 0.

\medskip
\medskip
As already stated, LU decomposition can be optimized for banded 
systems, both in terms of execution time:
$$
O(c_w n) \quad {\rm vs.} \quad O(n^3)
$$
and storage requirements
$$
O(n) \approx (2 {\tt kl} + {\tt ku} + 1) n \quad {\rm vs.} \quad O(n^2)
$$
%
where

\medskip
$$
c_w~\equiv~\hbox{\rm $w$-dependent constant}
$$

\medskip
$$
{\tt kl}~\equiv~\hbox{\rm number of lower diagonals}
$$

\medskip
$$
{\tt ku}~\equiv~\hbox{\rm number of upper diagonals}
$$

\medskip
Note that although $c_w$ is a fairly rapidly varying function of $w$
(must go to $n^2$ as $w$ goes to $n$), for any fixed bandwidth, the 
number of operations needed to solve the banded system is proportional
to the size, $n$,  of the system.

\newpage
{\bf LAPACK STORAGE SCHEME FOR BANDED MATRICES}

\medskip
\medskip
First note that the length of the $k$th sub- or super-diagonal is 
$n_k = n - k$.  Thus, particularly for $n \gg w$, $n_k \approx n$,
and we can efficiently store the matrix in a 2-d array in which 
the {\sl diagonals} of the matrix are stored in either rows or 
columns of the array. {\tt LAPACK} adopts the former strategy:
\begin{itemize}
\item Columns of matrix are stored in corresponding columns of 
      array (but elements are shifted)
\item Diagonals of matrix are stored in rows of array
\end{itemize}

\medskip
{\bf Example:} $n=5,~~{\tt kl}=2,~~{\tt ku}=1,~~w=4$

\medskip
\medskip
\medskip
$$
\left[
\begin{array}{ccccc}
 a_{11}  &  a_{12}  &          &          &          \\
 a_{21}  &  a_{22}  &  a_{23}  &          &          \\
 a_{31}  &  a_{32}  &  a_{33}  &  a_{34}  &          \\
         &  a_{42}  &  a_{43}  &  a_{44}  &  a_{45}  \\
         &          &  a_{53}  &  a_{54}  &  a_{55}  \\
\end{array}
\right]
\hbox{\rm~~~~is stored as~~~~}
\left[
\begin{array}{ccccc}
  \star  &  a_{12}  &  a_{23}  &  a_{34}  &  a_{45}  \\
 a_{11}  &  a_{22}  &  a_{33}  &  a_{44}  &  a_{55}  \\
 a_{21}  &  a_{32}  &  a_{43}  &  a_{54}  &  \star   \\
 a_{31}  &  a_{42}  &  a_{53}  &  \star   &  \star   \\
\end{array}
\right]
$$

\medskip
\medskip
\medskip
Here, $\star$, denotes an ``undefined'' element which will not 
be referenced by {\tt LAPACK} routines.

\medskip
\medskip
{\bf Additional Wrinkle:} To perform the LU decomposition with partial
pivoting and row interchanges (for numerical stability) {\tt LAPACK}
routines require {\tt kl} additional rows of storage. (Elements of 
these rows need not be assigned values prior to calling the solver.)
Thus, in the above example, we would actually need to set-up our 
banded array {\tt ab} as follows:

\medskip
\medskip
\medskip
$$
\left[
\begin{array}{ccccc}
    -    &    -     &    -     &    -     &    -     \\
    -    &    -     &    -     &    -     &    -     \\
  \star  &  a_{12}  &  a_{23}  &  a_{34}  &  a_{45}  \\
 a_{11}  &  a_{22}  &  a_{33}  &  a_{44}  &  a_{55}  \\
 a_{21}  &  a_{32}  &  a_{43}  &  a_{54}  &  \star   \\
 a_{31}  &  a_{42}  &  a_{53}  &  \star   &  \star   \\
\end{array}
\right]
$$
 
\medskip
\medskip
\medskip
Note that the leading dimension of this array (call it {\tt ab}) is thus
$$
{\tt ldab}~~=~~2\,{\tt kl} + {\tt ku} + 1
$$
%

\newpage
\medskip
{\bf FORTRAN IMPLEMENTATION: (illustrative purposes only)}

\medskip
\begin{verbatim}
      integer        maxn
      parameter    ( maxn = 1 000 )

      real*8         a(maxn,maxn)

      integer        kl,       ku,       ldab
      parameter    ( kl = 2,   ku = 1,   ldab = 6 )
      real*8         ab(ldab,maxn)
\end{verbatim}

\medskip
Assume that we have defined the (full) matrix {\tt a} (i.e. set 
{\tt a(i,j) for {\tt i, j~}$ = 1,\ldots,n$}), then the following 
code will pack the coefficients using the {\tt LAPACK} band-storage 
scheme.

\medskip
\begin{verbatim}
      do j = 1 , n
         do i = max( 1 , j - ku ) , min( n , j + kl )

            ab( kl + ku + 1 + i - j , j ) = a( i , j )

         end do
      end do
\end{verbatim}

\medskip
Note that the above indexing guarantees that (a) elements remain 
in their original columns and (b) diagonals ({\tt i} - {\tt j} = {\tt constant})
are stored in rows of {\tt ab}.  As we will see, in a real 
application, we will generally compute $A_{ij}$ and load the value into
{\tt ab(kl+ku+1+i-j,j)} without explicitly storing {\tt a(i,j)}.
However, the ``index translation'' remains the same.

\medskip
\medskip
\medskip
{\bf LAPACK GENERAL BANDED SOLVER: DGBSV}

\smallskip
\begin{verbatim}
      subroutine dgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )

         integer  n            ! size of system
         integer  kl, ku       ! # of lower and upper diagonals 
         integer  nrhs         ! Number of right hand sides (nominally 1)
         integer  ldab         ! Leading dimension of ab 
                               ! Note: ldab must be .ge. 2 * kl + ku + 1
         real*8   ab(ldab,*)   ! Matrix in LAPACK band-storage form
         integer  ldb          ! Leading dimension of rhs array
         real*8   b(ldb,*)     ! Right hand side(s) on input
                               ! Solution(s) on output
         integer  ipiv(n)      ! Storage for pivot vector
         integer  info         ! Usual LAPACK return code
\end{verbatim}

\newpage
{\bf EXAMPLE: 1D BVP}

\medskip
\medskip
\medskip
We again consider

\medskip
$$
u''(x) = f(x) \quad\quad
0 \le x \le 1 \quad\quad u(0), u(1)\,\,\,\, \hbox{\rm specified}
$$

\medskip
and adopt the same uniform discrete domain as previously, but now 
use the $O(h^4)$ (five-term) approximation to $u''(x)$ where possible, and the 
$O(h^2)$ (three-term) approximation elsewhere.  Then we have the 
following system of equations for the $n$ unknowns, $u_i$:

\medskip
\medskip
$$
\begin{array}{cccclccr}
u_1 & = & u(0) & & i = 1 & & (1)\\[0.10in]
h^{-2}(u_{i-1} - 2u_i + u_{i+1}) & = & f_i & &i = 2  & & (2)\\[0.10in]
(12h^2)^{-1}(-u_{i-2} + 16 u_{i-1} - 30 u_i + 16 u_{i+1} - u_{i+2}) &
 = & f_i & &i = 3, \ldots n - 2  & & (3)\\[0.10in]
h^{-2}(u_{i-1} - 2u_i + u_{i+1}) & = & f_i & &i = n-1  & & (4)\\[0.10in]
u_n & = & u(1) & & i = n  & & (5)\\[0.10in]
\end{array}
$$

\medskip
\medskip
Clearly, this set of equations is of the form 
$$
{\bf A} {\bf u} = {\bf b}
$$
where {\bf A} is an $n \times n$ pentadiagonal matrix.

\medskip
\medskip
Having derived the linear system we must solve, the only slightly 
challenging part of the exercise is correctly setting up the 
array {\tt ab} in {\tt LAPACK} banded form. Organizationally,
this process can 
be expedited by noting that all of the equations (1-5) can be written
in the form:

\medskip
\medskip
$$
	\sum_{k=-2}^{k=2} c_k u_{i+k} = b_i
$$

\medskip
\medskip
where, for each of the 5 sub-cases, the $c_k$ are constant coefficients,
and it is to be understood that the limits on the $k$-sum are to be changed 
if they result in references to $u_{i+k}$ such that $i+k < 1$ or $i+k>n$.

\medskip
\medskip
Specifically, we have

\medskip
\medskip
$$ c = [~0~,~0~,~1~,~0~,~0~] \eqno (1,5)$$
$$ c = h^{-2}\,\,[~0~,~1~,~-2~,~1~,~0~] \eqno (2,4)$$
$$ c = (12h^2)^{-1}\,\,[~-1~,~16~,~-30~,~16~,~-1~] \eqno (3)$$

\medskip
\medskip
Refer to the program source for {\tt bvp1d4.f} (available on-line) 
for full implementation details.
}
\end {document}
