\documentclass[11pt]{article}
\usepackage{makeidx}
%\usepackage{showidx}
\usepackage{latexsym}
%\usepackage{german}
%\usepackage{mypage}
%\usepackage{fancyhdr}
\usepackage{epsfig}
\usepackage{verbatim}
%\usepackage{ncode}
%\usepackage{myexercise}
\usepackage{amsfonts}
%\usepackage{rotating}
%\usepackage{fancybox}



%\input pagesizeLetter 
%%%% pagesizeLetter.tex %%%%%%%%%
%\setlength{\oddsidemargin}{-5mm}
%\setlength{\evensidemargin}{-5mm}
%\setlength{\headheight}{2em}
%\setlength{\headsep}{1cm}
%\setlength{\topmargin}{-15mm} %letter
%\setlength{\textheight}{235mm}
%\setlength{\textwidth}{173mm}


\input pagesizeA4

\input macros

%\selectlanguage{english}
\makeindex


\newcounter{twoyear}
\setcounter{twoyear}{\number\year}
\addtocounter{twoyear}{-2000}

\pagestyle{plain}
%  \pagestyle{fancy}
%  \lhead{\fancyplain{}{\sl\leftmark}}
%  \rhead{\rm\thepage}
%  \chead{}
%  \cfoot{}
%  %\rfoot{{\tiny SHA \number\day-\number\month-\number\year}}
%  \rfoot{{\tiny SHA \number\day-\number\month-0\thetwoyear}}
%  \lfoot{}

\title{Symmetric Banded Matrices}
\author{Version 0.1 December, 2001\\
Andreas Stahel
%\\\texttt{Andreas.Stahel@hta-bi.bfh.ch}\\
%\texttt{http://www.hta-bi.bfh.ch/\~{}sha}
}
\date{\today}
%\includeonly{}


\begin{document}
\newcommand{\currdir}{./}

\maketitle
%\input elements
%%%%%%%%%%%%%%%%%%%%%%%% Document %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\currdir{./}
\newcommand\ID[1]{\index{#1}}
\setlength{\marginparwidth}{25mm}

%\newcommand\eqref[1]{(\ref{#1})}


% to put a translation of a given word in the margin
\newcommand\DF[1]{$^*$\marginpar{#1}}
%\renewcommand\DF[1]{}

% \thispagestyle{empty}
% \dummy
% %\vspace{18cm}
% \vfill
% {\small
% \noindent
% \copyright Andreas Stahel, 2000\\
% All rights reserved. This work may not be translated or copied in whole or in
% part without the written permission by the author, except for brief excerpts in
% connection with reviews or scholarly analysis. Use in connection with any form
% of information storage and retrieval, electronic adaptation, computer software
% is forbidden.
% }
% \newpage

%\pagenumbering{roman}
\tableofcontents
%\addcontentsline{toc}{section}{List of Figures}
%\listoffigures

%\addcontentsline{toc}{chapter}{List of Tables}
%\listoftables

%\newpage
%\pagenumbering{arabic}


\section{Basic description}
Many  matrices used to solve PDE (using FEM) are symmetric. It the
nodes are numbered properly then the matrix will show
a band structure, i.e. all nonzero elements are located close to the main
diagonal. The algorithm of Cholesky or the $LDL^T$ factorization can take
advantage of this structure, see~\cite{GoluVanLoan96}. For a symmetric
matrix $A$ of size $n\times n$ with semi-bandwidth $b$ the approximate
computational cost to solve one system of equations is 
given by
\[ \mbox{Gauss}\approx \frac{1}{3}\;n^3 \text{and}
 \mbox{Band Cholesky}\approx \frac{1}{2}\,n\,b^2\]
Obviously for $b\ll n$ is is advantageous to
use a banded solver. A more detailed analysis and an implementation is
given in~\cite{VarFem}.

To take advantage of the symmetry and the band structure the matrices will
be stored in a modified format, as illustrated below.
\[\Matrix{5}{10&2&3&0&0\\2&20&4&5&0\\3&4&30&6&7\\0&5&6&40&8\\0&0&7&8&50}
\longrightarrow
\Matrix{3}{10&2&3\\20&4&5\\30&6&7\\40&8&0\\50&0&0}\]
A banded version of the $LDL^T$ factorization in~\cite{GoluVanLoan96}
can be implemented. If the matrix $A$ is strictly positive definite, then
the algorithm is known to be stable. If $A$ is not positive definite, then
problems might occur, since no pivoting is done. The matrix $A$ is
positive definite if and only if the diagonal matrix $D$ is positive.

For a given matrix some of its smallest eigenvalues can be computed with an
algorithm based on inverse power iteration. Precise information on the
numerical errors is provided. The code is capable of finding eigenvalues of
medium size matrices, where the standard command \texttt{eig()} is either very
slow of will fail.

%This author prefers the notation of a $R^TDR$ factorization, the
%difference is in notation only, i.e. $R^T=L$. 
%The presented code can be obtained form authors home
%page\footnote{\texttt{http://www.hta-bi.bfh.ch/\~{}sha}}.

\begin{table}[htbp]
  \begin{center}
\begin{tabular}{|l|l|}\hline
\multicolumn{2}{|l|}{Operations for symmetric, banded matrices}\\\hline
\texttt{SBSolve()} & solve a system of linear equations \\
\texttt{SBFactor()} & find the $R^TDR$ factorization \\
\texttt{SBBacksub()} & use back-substitution to solve system of equations\\
\texttt{SBEig()} & find a few of the smallest eigenvalues and eigenvectors\\
\texttt{SBProd()} & multiply symmetric banded matrix with full matrix\\
\texttt{FullToBand()} & convert a symmetric matrix to a banded matrix\\
\texttt{BandToFull()} & convert a banded matrix to a symmetric matrix\\
\texttt{BandToSparse()} & convert a banded matrix to a sparse matrix\\
\hline\end{tabular}
    \caption{List of commands}
    \label{tab:commands}
  \end{center}
\end{table}
\begin{center}
\end{center}


\section{Description of the commands}
 
\subsection{\texttt{SBSolve}}\ID{SBSolve}
The basic factorization algorithm is implemented in \texttt{SBSolve}. The
function can return the solution of the system of linear equations, or the
solution and the factorization of the original matrix.
Multiple sets of equations can be solved.

\begin{verbatim}
[...] = SBSolve(...)
  solve a system of linear equations with a symmetric banded matrix

  X=SBSolve(A,B)
  [R,X]=SBSolve(A,B)

   solves A X = B

   A is mxt where t-1 is number of non-zero super-diagonals
   B is mxn
   X is mxn
   R is mxt

  if A would be ! 11000 ! then A= ! 11 ! 
                ! 14300 !         ! 43 ! 
                ! 03520 !         ! 52 ! 
                ! 00285 !         ! 85 ! 
                ! 00059 !         ! 90 ! 

  B is a full matrix

  The code is based on a LDL' decomposition (use L=R'), without pivoting.
  If A is positive definite, then it reduces to the Cholesky algorithm.

  R is an upper right band matrix
  The first column of R contains the entries of a diagonal matrix D. 
  If the first column of R is filled by 1's, then we have R'*D*R = A
\end{verbatim}

To determine the \ID{matrix, inverse}inverse matrix $\mathtt{A}^{-1}$ one can
use the command \texttt{invA = SBSolve(A,eye(n));}.  Be aware that calculating
the inverse matrix is rarely a wise thing to do. Most often the inverse of a
banded matrix will loose the band structure.

If the matrix \textbf{A} is strictly positive definite, then the algorithm is
stable and one can expect the solution to be as accurate as the conditions
number of \textbf{A} permits. If \textbf{A} is semidefinite, then large errors
might occur, since \textbf{not pivoting} is implemented in the code. The
matrix is positive definite iff all eigenvalues are positive, this can be
verified by inspection the sign of the numbers in the first column of
\textbf{R}.  The matrix is positive definite if the first column of the
factorization matrix \texttt{R} (use \texttt{SBFactor()}) contains positive
numbers only. A description of the algorithm can be found
in~\cite{GoluVanLoan96} or~\cite{VarFem}.

\subsection{\texttt{SBFactor} and \texttt{SBBacksub}}
\ID{SBFactor}\ID{SBBacksub}
Instead of calling \texttt{X=SBSolve(A,B)} one can first call
\texttt{R=SBFactor(A)} to determine the factorization $A=R^TDR$ and
then \texttt{B=SBBacksub(R,X)} to solve the system(s) $A\cdot X=B$~.
Since most of the computational effort is in the factorization, this can be
useful if many system of linear equations have to be solved sequentially.
If multiple system are to be solved simultaneously it is preferable to use
\texttt{SBSolve(A,B)} with a matrix \texttt{B}~.

\begin{verbatim}
[...] = SBFactor(...)
  find the R'DR factorization of a symmetric banded matrix

  R=SBFactor(A)

   A is mxt where t-1 is number of non-zero super diagonals
   R is mxt

  if A would be ! 11000 ! then A= ! 11 !
                ! 14300 !         ! 43 !
                ! 03520 !         ! 52 !
                ! 00285 !         ! 85 !
                ! 00059 !         ! 90 !


  The code is based on a LDL' decomposition (use L=R'), without pivoting.
  If A is positive definite, then it reduces to the Cholesky algorithm.

  R is an upper right band matrix
  The first column of R contains the entries of a diagonal matrix D.
  If the first column of R is filled by 1's, then we have R'*D*R = A
\end{verbatim}

\begin{verbatim}
[...] = SBBacksub(...)
  using backsubstitution  to return the solution of a system of linear equations

  X=SBBacksub(R,B)

   B is mxn
   X is mxn
   R is mxt

   R is produced by a call of [X,R] = SBSolve(A,B) or R = SBFactor(A)
   It is an upper right band matrix
   The first column of R contains the entries of a diagonal matrix D.
   If the first column of R is filled by 1's, then we have R'*D*R = A
\end{verbatim}


If there is interest in the classical Cholesky decomposition
\ID{Cholesky decomposition}  of the matrix \texttt{A} 
(i.e. $\mathtt{A}=\mathtt{R}^\prime\cdot \mathtt{R}$) then \texttt{R} can be
computed by
\begin{verbatim}
rBand=SBFactor(A);
d=sqrt(rBand(:,1));
rBand(:,1)=ones(n,1);
r=triu(diag(d)*rBand)
\end{verbatim}

The number of positive/negative numbers in the first column of \textbf{R}
equals the number of positive/negative eigenvalues of \textbf{A}.


\subsection{\texttt{SBEig}}\ID{SBEig}
For given symmetric matrices \textbf{A} and \textbf{B} the standard
(resp. generalized) eigenvalue problem will be solved, i.e.
\[ \mathbf{A}\,\vec v=\lambda\,\vec v \text{resp.}
   \mathbf{A}\,\vec v=\lambda\,\mathbf{B}\,\vec v  \]

Using inverse power iteration a given number of the smallest (absolute value)
eigenvalues if a symmetric matrix \textbf{A} are computed. If needed the
eigenvectors are also generated.  A set of initial vectors \textbf{V} have to
be given. If those are already close to the eigenvectors, then the algorithm
will converge rather quickly. For a precise description and analysis
consult~\cite{GoluVanLoan96}.
\begin{verbatim}
[...] = SBEig(...)
  find a few eigenvalues of the symmetric, banded matrix
  inverse power iteration is used for the standard and generalized
  eigenvalue problem

  [Lambda,{Ev,err}] = SBEig(A,V,tol)     solve A*Ev = Ev*diag(Lambda)
                    standard eigenvalue problem

  [Lambda,{Ev,err}] = SBEig(A,B,V,tol)   solve A*Ev = B*Ev*diag(Lambda)
                    generalized eigenvalue problem

   A   is mxt, where t-1 is number of non-zero superdiagonals
   B   is mxs, where s-1 is number of non-zero superdiagonals
   V   is mxn, where n is the number of eigenvalues desired
       contains the initial eigenvectors for the iteration
   tol is the relative error, used as the stopping criterion

   X   is a column vector with the eigenvalues
   Ev  is a matrix whose columns represent normalized eigenvectors
   err is a vector with the a posteriori error estimates for the eigenvalues
\end{verbatim}

The algorithm is based on inverse power iteration 
\ID{iteration, inverse power}  with $n$ independent vectors.
The iteration will proceed until the relative change of all eigenvalues is
smaller than the given value of \texttt{tol}. This does not guarantee that the
relative error is smaller than \texttt{tol}.  The initial guesses \textbf{V}
for the eigenvectors have to be linearly independent. The closer the initial
guess is to the actual eigenvector,the faster the algorithm will converge. 
The algorithm returns $n$ eigenvalues closest to~0~.


For the standard eigenvalue problem $\mathbf{A}\;\vec v_i = \lambda_i\, \vec
v_i$ the eigenvectors $\vec v_i$ will be orthonormal with respect to the
standard scalar product, i.e, $\langle \vec v_i\,,\,\vec v_j\rangle
=\delta_{i,j}$. For the generalized eigenvalue problem 
$\mathbf{A}\;\vec v_i = \lambda_i\, \mathbf{B}\,\vec v_i$ this translates to
$\langle \vec v_i\,,\,\mathbf{B}\,\vec v_j\rangle=\delta_{i,j}$. The symmetric
matrix \textbf{B} should be positive definite. The columns of \texttt{Ev} can
be used to restart the algorithm if higher accuracy is required.

The algorithm will return reliable estimates for the errors in the eigenvalues.
The \`a posteriori\ID{a posteriori estimate} error estimate is based on the
residual $\vec r = \mathbf{A}\,\vec v -\lambda\,\vec v$ and
\[ \min_{\lambda_i\in\sigma(\mathbf{A})}|\lambda-\lambda_i|\leq 
  \langle \vec r\,,\,\vec r\rangle=\|\vec r\|\]
where we use the normalization $\langle \vec v\,,\,\vec v\rangle =1$.
If one of the eigenvalues has to be computed with high accuracy, the
approximate value $\lambda$ may be subtracted from the diagonal of the matrix.
Then the eigenvalue closest to zero of the modified matrix
$\mathbf{A}-\lambda\id$ can be computed, using the already computed
eigenvector. If the eigenvalue is isolated the algorithm will converge very
quickly. This algorithm is similar to the Rayleigh 
\ID{Rayleigh quotient iteration} quotient iteration. A good description is
given in~\cite{GoluVanLoan96}.
\medskip

If the eigenvalue closest to $\lambda$ is denoted by $\lambda_i$ we have
the improved estimate
\[ |\lambda-\lambda_i|\leq  \dfrac{\|\vec r\|^2}{\mbox{gap}}\text{where}
 \mbox{gap}=\min\{|\lambda-\lambda_j|\;:\;
  \lambda_j\in\sigma(\mathbf{A}),j\neq i\} \]
It is very easy to implement this test in \Octave{}. If the estimate is
based on approximate values of the eigenvalues, then the result is not as
reliable as the previous one. Since the value of \texttt{gap} will carry an
approximation error. The situation is particularly bad if some eigenvalues are
clustered. A code  sample is provided.

\bigskip

For the generalized eigenvalue problem \ID{eigenvalue problem, generalized}
we use the residual
$\vec r = \mathbf{A}\,\vec v -\lambda\,\mathbf{B}\,\vec v$
and the estimate
\[ \min_{\lambda_i\in\sigma(\mathbf{A})}|\lambda-\lambda_i|\leq 
  \sqrt{\langle \vec r\,,\,\mathbf{B}^{-1}\vec r\rangle}\text{and}
|\lambda-\lambda_i|\leq  \dfrac{\|\vec r\|^2}{\mbox{gap}}\]
where we use the normalization 
$\langle \vec v\,,\,\mathbf{B}\,\vec v\rangle =1$.
The precise algorithm and proof of the above estimate is given
in~\cite{VarFem}.

\subsection{\texttt{SBProd}}
With this command a symmetric banded matrix can be multiplied with a full
matrix.
\begin{verbatim}
[...] = SBProd(...)
  multiplies a symmetric banded matrix with a matrix

  X=SBProd(A,B)

   A is mxt where t-1 is number of non-zero super diagonals
   B is mxn
   X is mxn

  if A would be ! 11000 ! then A= ! 11 ! 
                ! 14300 !         ! 43 ! 
                ! 03520 !         ! 52 ! 
                ! 00285 !         ! 85 ! 
                ! 00059 !         ! 90 ! 

  B is full matrix Ax=B
\end{verbatim}
\subsection{\texttt{BandToFull}, \texttt{FullToBand} and \texttt{BandToSparse}}
With these commands conversion between full, symmetric matrices and banded
symmetric matrices is possible. A conversion to a sparse format is also
included.

%\newpage

\bibliographystyle{keysort}
%\bibliographystyle{apalike}
%\bibliographystyle{siam}
\addcontentsline{toc}{section}{Bibliography}
\bibliography{school,book,fem,new}


%\clearpage



%\printindex
%\addcontentsline{toc}{section}{Index}

\end{document}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
