\documentclass[landscape,dvips]{foils}

\usepackage{amsmath,amsthm,amssymb,amscd}
\usepackage{alltt}
\usepackage{epsf}
\usepackage{epic}
\def\w{{\cal W}}
\usepackage{latexsym}
\usepackage{times}
\usepackage{bookman}
\usepackage{graphics}
\usepackage{color}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage[ps2pdf]{hyperref}

\newcommand{\img}[2]{{\\epsfxsize=#1truecm\epsfbox{Images/PS/#2}}}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}
\newtheorem{definition}{Definition}
\newtheorem{proposition}{Proposition}

\def\basswidt{4.0in}
\def\gutter{0.2in}

\newenvironment{comment}{\setbox0=\vbox\bgroup}{\egroup}
\newenvironment{singlespace}{\def\baselinestretch{1.0}\large\normalsize}{\par}

\setlength{\textheight}{8in}
\setlength{\textwidth}{9in}
\voffset=-0.9in
\renewcommand\arraystretch{1.15}

\newcommand{\ie}{{\sl i.e.}}

\newcommand{\blue}[1]{{\color{blue} #1}}
\definecolor{mygreen}{rgb}{0.1,0.6,0.1} 
\newcommand{\green}[1]{{\color{mygreen} #1}}
\newcommand{\red}[1]{{\color{red} #1}}
\newcommand{\stitle}[1]{\foilhead{\red{\noindent{#1}}}\vspace{-0.5in}}
\newcommand{\bigitem}[2]{\item{#1\begin{small}\begin{itemize}#2\end{itemize}\end{small}}}
\newcommand{\mcite}[2]{{\small{\green{[#1, #2]}}}}
\newcommand{\eg}{{\sl e.g.}}
\newcommand{\etal}{{\sl et al.}}

\begin{document}

\title{\red{
Parallel banded preconditioners for non-symmetric linear system solvers
}}

\author{\\\\
{Murat Manguoglu, Ahmed Sameh, and Ananth Grama}\\
{Department of Computer Science, Purdue University}\\\\\\
{\tt \{mmanguog,sameh,ayg\}@cs.purdue.edu}
}

\date{June 12, 2008}

\maketitle

\vspace{1.0in}
\noindent
{\small{\blue{Acknowledgment:} Funding for this work was provided by 
the US National Science Foundation}}

\thispagestyle{empty}

\stitle{Theses of Talk}

\begin{itemize}
{\em
\item
Building on our prior work demonstrating the effectiveness of Spike
as a solver and preconditioner, we derive highly scalable parallel
formulations and demonstrate performance on large parallel platforms.
\item
Using a pseudo-analytical performance model, we derive estimates
of parallel time, show these estimates to be highly accurate, and
use these estimates to argue scaling of Spike to much larger machine
configurations.
}
\end{itemize}


\stitle{Layout of Presentation}

\begin{itemize}
\item
Summary of Spike Performance (Serial)
\item
Spike Parallelization Strategy
\item
Parallel Performance Results
\item
Pseudoanalytical Performance Model for Spike
\item
Validation of Performance Model
\item
Limitations of Approach
\end{itemize}


\stitle{Summary of Spike Solver}

\begin{itemize}
\item
Targeted to banded, or low-rank perturbations of banded systems
(dense or sparse within the band).
\item
Banded approximations used as preconditioners.
\item
Spike is designed to optimize memory system as well as parallel performance.
\end{itemize}


\stitle{Summary of Spike Solver}

\noindent
Solving $Ax = F$ with four partitions:

\centerline{\psfig{file=spike1.ps,width=7.0in}}
Partitioning of the matrix $A$ and the RHS $F$, with $p=4$.


\stitle{Summary of Spike Solver}

\begin{itemize}
\item
Spike factorization is: $A= D \ S$, where $D$
is a block-diagonal matrix consisting only of the diagonal blocks
$A_j$,

$$D = \mathrm{diag}(A_1,..., A_p),$$

\item
Matrix $S$ has the following structure:

\centerline{\psfig{file=spike2.ps,width=7.0in}}

\end{itemize}


\stitle{Summary of Spike Solver}

\begin{itemize}
\item
The spikes $V_j$ and $W_j$ are given by
\begin{equation}
\begin{array}{l}
V_j=(A_j)^{-1}\left[\begin{array}{ll}0 \\
I_m\end{array}\right]B_j, %\quad \mbox{for} %\quad
% j>1
 \quad \mbox{and} \quad %\\[15pt]
W_j=(A_j)^{-1}\left[\begin{array}{ll}I_m \\ 0\end{array}\right]C_j. %\quad \mbox{for} \quad
%j<p.
\end{array}
\label{eq_spikes}
\end{equation}

\item
Spikes $V_j$ and $W_j$, $j=2,...,p-1$, are generated by
solving,
\begin{gather}
A_j
\begin{bmatrix}V_j,W_j\end{bmatrix}
=
\begin{bmatrix}0 & C_j    \\ \vdots  &  0  \\0 & \vdots \\B_j &0 \end{bmatrix}.
%\begin{bmatrix}0 & \dots & 0  & B_j    \\
% C_j  & 0 & \dots & 0  \end{bmatrix}^T.
\label{eq_solveA}
\end{gather}
\item
Solving the system $AX=F$ now reduces to two
steps: $D G = F$ and $S X = G$.

\end{itemize}



\stitle{Spike Performance (Serial)}

\begin{itemize}
\item
Spike used as a preconditioner with BiCGStab as the iterative solver.
\item
Comparison with:
\begin{itemize}
\item
ILUPACK: Multilevel ILU [Bollhofer] \\
\url{http://www.math.tu-berlin.de/ilupack/}
\item
ILUT: Incomplete LU Factorization from Sparsekit [Saad] 
\url{http://www-users.cs.umn.edu/~saad/software/SPARSKIT/sparskit.html}
\item
ILUT-I: Improved ILUT [Benzi] 
(reorder using HSL-MC64 to maximize the product of the diagonals and
scale the matrix, apply symmetric RCM reordering,
incomplete factorization via ILUT)
\end{itemize}
\end{itemize}


\stitle{Spike-Based Preconditioning: Preprocessing}

\noindent
Extracting a banded preconditioner:

\begin{itemize}
\item
reorder using HSL-MC64 to make the diagonal zero free
\item
reorder $\|A\| + \|A^T\|$ using HSL-MC73 to place larger elements
closest to the main diagonal
\item
extract a banded preconditioner, such that 99.9\% of the weight is
inside the band 
\item
factorize the banded preconditioner
\end{itemize}


\stitle{Spike Performance}

\noindent
Test Problems: \\ \\

\centerline{\psfig{file=spike3.ps,width=8.0in}}


\stitle{Spike Performance}

\noindent
Comparison to ILUPACK AMF/PQ preconditoners on an uniprocessor \\ \\ \\


\centerline{\psfig{file=spike4.ps,width=8.0in}}


\noindent
Outer Iterative Solver: unrestarted GMRES, ILUPACK Parameters:
droptol : 1e-1 , bound for inv(L), inv(U): 10 , elbow space: 100




\stitle{Spike Performance}

\noindent
Comparison to ILUT and Improved-ILUT Preconditioners on an uniprocessor \\ \\ \\




\centerline{\psfig{file=spike5.ps,width=8.0in}}

\noindent
Outer Iterative Solver : BiCGStab


\stitle{Spike Performance}

\noindent
WSO solver speedup, SGI Altix: \\ \\ \\

\centerline{\psfig{file=spike6.ps,width=8.0in}}


\stitle{Spike Performance}

\noindent
Reservoir Simulation (SPE10 benchmarks): Problem 1: $N=2,244,000$ 
Problem 2 : $N=2,462,265$. \\ 

\centerline{\psfig{file=spike7.ps,width=8.0in}}


\stitle{Spike Performance}

\noindent
Spike Factorization Performance: \\ \\

\centerline{\psfig{file=spike8.ps,width=8.0in}}


\stitle{Spike Performance}

\noindent
Reservoir Simulation Benchmark 2: \\ \\

\centerline{\psfig{file=spike9.ps,width=8.0in}}










\stitle{Spike: Parallel Steps}

\begin{enumerate}
 \item Compute the $LU$ (via DDBTRF)
 \begin{itemize}
    \item $L_j U_j \leftarrow A_j$ for $j=1,2,3,4$
 \end{itemize}


 \item Compute spikes (via DTBTRS)
  \begin{itemize}
   \item   Solve for $V_j$:
      $L_j U_j V_j =
       \left[ \begin{matrix}
       0 \; \cdots \;
     0 \;
         B_j^T
        \end{matrix}\right]^T
     $ for $j=1,2,3$
  \item  Solve for $W_j$:
           $L_j U_j W_j=
          \left[ \begin{matrix}
               C_j^T \;
           0 \; \cdots \;
         0 \end{matrix}\right]^T$ for $j=2,3,4$.
   \end{itemize}
 \item Communicate spike tips $W^{(t)}_j$ (processor $j$ sends to processor
$j-1$, for $j=2,3,4$). $W^{(t)}_j$ is
 the top $k \times k$ part of $W_j$.


 \item Factorize the reduced system (via DGETRF)
     \begin{equation}
       \begin{pmatrix}
           I & V^{(b)}_j \\
           W^{(t)}_j & I
      \end{pmatrix}
     \end{equation}

\newpage

 \item Modify the right hand side  by solving (via DTBTRS)  :
       $L_j U_j g_j = f_j$ ($j=1,2,3,4$)
 \item Communicate the modified right hand side tips $g^{(t)}_j$ (processor $j$
sends to processor $j-1$ for $j=2,3,4$).
 \item Solve the reduced system
        \begin{equation}
          \begin{pmatrix}
              I & V^{(b)}_j \\
              W^{(t)}_{j+1} & I
          \end{pmatrix}
          \begin{pmatrix}
             x^{(b)}_j \\
             x^{(t)}_{j+1}
          \end{pmatrix} =
          \begin{pmatrix}
             g^{(b)}_j \\
             g^{(t)}_{j+1}
          \end{pmatrix}
         \end{equation}
\item Communicate the reduced system solution  $g^{(t)}_{j+1}$ (processor $j$
sends to processor $j+1$ for $j=1,2,3
$).
\item Retrieve $x_j$ ($j=1,2,3,4$) (via DGEMV)
   $x_j =f_j -V_j x^{(t)}_{j+1} - W_j x^{(b)}_{j-1} $  ($V_4 = 0$ and $W_1 = 0$)
\end{enumerate}





\stitle{Analytical Characterization of Spike Steps}

\noindent
Notation:
\begin{itemize}
\item
$k=kl=ku$ (upper and lower bandwidths identical)
\item the dimension of the matrix is $N$
\item the dimension of the partitioned blocks are $n = N /p$.
\end{itemize}


\begin{table}[h!b!p!]
%\caption{Truncated Spike Algorithm Steps}
\small
\begin{tabular}{lll}
\hline
Stage & Description & Cost            \\\hline
1    &   Factorize the Diagonal Blocks & $\alpha_1  n k^{2} + \beta_1  n  k$   \\
2    &   Compute  Spikes               & $\alpha_2  n k^{2} + \beta_2  n  k$         \\
3    &   Communicate Tips of Spikes    & $\alpha_3  k^{2} (p-1) + \beta_3 k^{2} + \gamma_3 $         \\
4    &   Factorize The Reduced System  & $\alpha_4 k^{3} + \beta_4 k^{2} $             \\
5    &   Modify the Right Hand Side    & $\alpha_5 n k $    \\
6    &   Communicate Tips of MRHS      & $\alpha_6 k (p-1) + \beta_6 k + \gamma_6  $         \\
7    &   Solve the Reduced System      &$\alpha_7 k^{2}$  \\
8    &   Comm. Soln. of Reduced System & $\alpha_8 k (p-1) + \beta_8 k + \gamma_8  $ \\
9    &   Retrieve the Solution         & $\alpha_9 n k + \beta_9 n $ \\
\hline
\end{tabular}
\label{table1}
\end{table}


\stitle{Analytical Characterization of Spike Steps: Observations}

\begin{itemize}
\item
Stage 1 and Stage 2  cost $O(n k^{2})$ computation and $O(n  k)$ memory
references.
\item
Stage 4 has $O(k^{3})$ computation an $O(k^{2})$ memory references.
\item
Stage 5 has $O(n k)$ computation and $O(n k)$ memory references.
\item
Stage 7 has $O(k^{2})$ computation and $O(k^{2})$ memory references.
\item
Stage 9 has $O(n k)$ computation for matrix vector product and $O(n)$
computation for vector addition.
\item
In Stages $3,6,8$ we model the cost by $O(data)$ for the communication
of data, $O(data \times (p-1))$ for network saturation
\end{itemize}


\stitle{Training the Analytical Model}

\begin{itemize}
\item
Digonally dominant toeplitz system of dimension $5,000,000$.
\item
Training platform: Ranger Sun Constellation Linux Cluster at TACC ($3,936$
nodes, each node has 16 cores of AMD Barcelona Processor. Interconnect
is Infiniband.
\item
Train using $16, 32, 64$ processors ( $1,2,4$ nodes respectively) for
bandwidths $k=15,25, 35$.
\item
Using a least squares approximation we find the following parameters
$\alpha_i,\beta_i , \gamma_i$ for $i=1,..,9$.
\end{itemize}

\stitle{Cost Model Parameters}

\noindent
Training yields the following parameters:

\begin{table}[h!b!p!]
\centering
\begin{tabular}{llll}
\hline
i    & $\alpha_i$ & $\beta_i$   & $\gamma_i$           \\\hline
1    & $8.62\times 10^{-10}$ &  $2.61\times10^{-8}$&   -         \\
2    & $3.96\times 10^{-8}$  &  $3.66\times10^{-7}$&   -         \\
3    & $9.07\times 10^{-14}$ &  $1.21\times10^{-6}$&  $9.48\times 10^{-4}$\\
4    & $1.65\times 10^{-29}$ &  $7.43\times10^{-6}$&   -         \\
5    & $2.85\times 10^{-8}$  &     -       &   -         \\
6    & $4.31\times 10^{-9}$  &  $4.10\times10^{-7}$&  $3.59\times 10^{-5}$\\
7    & $1.19\times 10^{-7}$  &      -              &   -         \\
8    & $2.57\times 10^{-7}$  &  $3.39\times10^{-31}$& $6.42\times 10^{-4}$\\
9    & $8.86\times 10^{-8}$  &  $9.51\times10^{-9}$&  -           \\
\hline
\end{tabular}
\label{table2}
\end{table}

\noindent
{\bf Note:} Parameters also tell us what regime our algorithm operates in: is
the computation memory or processor bound? is communication latency or
bandwidth bound, etc.


\stitle{Cost Model Verification}

\noindent
Validation on two systems -- one using the same toeplitz system used in
training, the other, a toeplitz system of dimension $10,000,000$.
We verify the model for $128$, $246$, $512$, and $1,024$ processors.

\begin{table}[h!b!p!]
\centering
\begin{tabular}{llllll}
\hline
N               & k     &  p      & Observed  & Model   & Error   \\\hline \hline
$5,000,000$     & $35$  & $128$   & $2.78$    & $2.64$  & $0.13$  \\
$5,000,000$     & $25$  & $128$   & $1.55$    & $1.49$  & $0.06$  \\
$5,000,000$     & $15$  & $128$   & $0.70$    & $0.66$  & $0.04$  \\\hline
$5,000,000$     & $35$  & $256$   & $1.49$    & $1.33$  & $0.16$  \\
$5,000,000$     & $25$  & $256$   & $0.79$    & $0.75$  & $0.04$  \\
$5,000,000$     & $15$  & $256$   & $0.35$    & $0.33$  & $0.02$  \\\hline
$5,000,000$     & $35$  & $512$   & $0.67$    & $0.67$  & $0.00$  \\
$5,000,000$     & $25$  & $512$   & $0.38$    & $0.38$  & $0.00$  \\
$5,000,000$     & $15$  & $512$   & $0.20$    & $0.17$  & $0.03$  \\\hline
$5,000,000$     & $35$  & $1,024$ & $0.37$    & $0.35$  & $0.02$  \\
$5,000,000$     & $25$  & $1,024$ & $0.21$    & $0.20$  & $0.01$  \\
$5,000,000$     & $15$  & $1,024$ & $0.10$    & $0.09$  & $0.01$  \\\hline \hline
\end{tabular}
\label{table3}
\end{table}


\stitle{Cost Model Verification (contd)}

\begin{table}[h!b!p!]
\centering
\begin{tabular}{llllll}
\hline
N               & k     &  p      & Observed  & Model   & Error   \\\hline \hline
$10,000,000$    & $35$  & $128$   & $5.36$    & $5.27$  & $0.09$  \\
$10,000,000$    & $25$  & $128$   & $3.03$    & $2.98$  & $0.06$  \\
$10,000,000$    & $15$  & $128$   & $1.36$    & $1.31$  & $0.05$  \\\hline
$10,000,000$    & $35$  & $256$   & $2.78$    & $2.64$  & $0.13$ \\
$10,000,000$    & $25$  & $256$   & $1.55$    & $1.49$  & $0.06$ \\
$10,000,000$    & $15$  & $256$   & $0.68$    & $0.66$  & $0.02$ \\\hline
$10,000,000$    & $35$  & $512$   & $1.51$    & $1.33$  & $0.18$ \\
$10,000,000$    & $25$  & $512$   & $0.79$    & $0.75$  & $0.04$ \\
$10,000,000$    & $15$  & $512$   & $0.35$    & $0.33$  & $0.02$ \\\hline
$10,000,000$    & $35$  & $1,024$ & $0.69$    & $0.68$  & $0.01$ \\
$10,000,000$    & $25$  & $1,024$ & $0.45$    & $0.38$  & $0.07$ \\
$10,000,000$    & $15$  & $1,024$ & $0.19$    & $0.17$  & $0.02$ \\
\hline
\end{tabular}
\label{table3}
\end{table}


\stitle{Cost Model Verification}

\begin{itemize}
\item
Model yields outstanding accuracy in predicting performance well beyond
training machine size.
\item
Model predicts linear scaling (linearly increasing speedup when keeping
problem size per processor constant)!
\end{itemize}

\noindent
{\em Spike is an excellent candidate for emerging ultrascale platforms.}




\stitle{Limitations of Approach}

\begin{itemize}
\item
Serial performance is harder to characterize than parallel performance!
  \begin{itemize}
  \item
	Operate in stable region of serial performance
  \end{itemize}
\item
Analytical models may vary across platforms
  \begin{itemize}
  \item
	Separate models for message passing and shared address space machines
  \end{itemize}
\item
Granularity of aggregates
  \begin{itemize}
  \item
	It may not always be possible to clearly identify and analytically
	characterize various steps in a parallel algorithm.
  \end{itemize}
\end{itemize}



\stitle{Thank You!}



\end{document}


\stitle{Banded Preconditioners for Iterative Methods}

\begin{itemize}
\item
Derive banded approximations to the matrix, which can act as good
preconditioners.
\item
Bands must be narrow and capture much of the matrix norm.
\item
Use banded solvers with high FLOP counts and concurrency.
\end{itemize}

\noindent
{\bf Key questions:}

\begin{itemize}
\item
How do we derive such narrow-banded preconditioners (non-symmetric
permutations)?
\item
Can such simple preconditioners be competitive against
traditional preconditioners in terms of iteration counts, FLOPS,
FLOP counts, and parallelism?
\end{itemize}



\stitle{Contributions and Results}

\begin{itemize}
\item
Banded preconditioners (with suitable reordering) can significantly
outperform ILU preconditioners in terms of iteration counts, FLOP
counts, as well as concurrency for large classes of matrices!
\item
Reordering schemes based on weighted spectral methods are highly
effective in deriving narrow banded preconditioners.
\item
The overhead of such reordering schemes is easily offset by the
lower solution cost for the system.
\item
A number of banded solvers (LAPACK, Spike) can be used for the
inner solve.
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Bandwidth Reduction}

\begin{itemize}
\item{Traditional algorithms (\eg, Cuthill-McKee \mcite{Cuthill \& McKee}{1969}, 
      Spectral reordering \mcite{Barnard \etal}{1995}) are aimed at minimizing the
      bandwidth
      \begin{displaymath}
       BW(A) = \max_{i,j: A(i,j) > 0}{|i-j|}
      \end{displaymath}
}

\bigitem {\red{Heavy} (high-magnitude) nonzeros that are \blue{distant} from the 
diagonal may significantly \red{degrade} the performance (convergence 
rate)}{
\item Particularly, for ill conditioned matrices.
}

%\item{
%We define {\bf \blue{weighted bandwidth}} as
%\begin{displaymath}
%WBW(A) = \frac{\sum_{i,j}{|A(i,j)|\times|i-j|}}{||A||_F}
%\end{displaymath}
%}

\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\stitle{Accounting for Heavy Entries}

\begin{itemize}
\item{We \green{generalize} the definition of bandwith}

\item{For given $b$, we define \blue{bandweight} as
\begin{displaymath}
w_b(A) = \sum_{i,j: |i-j| < b}{|A(i,j)|}
\end{displaymath}
}
\bigitem{Then, for given $\alpha$, we define \blue{$\alpha$-bandwidth} as
the \green{smallest bandwidth} that encapsulates an $\alpha$ fraction of \red{total
matrix weight}
\begin{displaymath}
BW_{\alpha}(A) = \min{b}~\text{such that}~{w_b(A) \geq \alpha \times \sum_{i,j}{|A(i,j)|}} 
\end{displaymath}
}
{\item Observe that this is a generalization of bandwidth, such that $BW_1(A) = BW(A)$}

\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Spectral Ordering}

\begin{itemize}

\item{Commonly used in graph-theoretic applications and matrix algorithms}

\bigitem{
Find $x$ that \red{minimizes}
\begin{displaymath}
\sum_{i,j: A(i,j) > 0}{(x(i)-x(j))^2}
\end{displaymath}}
{\item Reorder rows and columns of $A$ accordingly
}

\bigitem{The eigenvector that corresponds to the smallest non-zero eigenvalue of the \blue{Laplacian} 
\begin{displaymath}
\begin{array}{l}
L(i,j) = -1~~~~~~~~~~~~~~~~~~~~~~~~~
{\bf if}~i \neq j \wedge A(i,j) > 0\\
L(i,i) = |\{j: A(i,j) > 0 \}|,
\end{array}
\end{displaymath}
a.k.a, \blue{Fiedler vector}, minimizes this  cost function \mcite{Fiedler}{1973}}{
\item{Can be computed effectively using iterative techniques (\eg, CG \mcite{Kruyt}{1995})}
}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Weighted Spectral Ordering}

\begin{itemize}
\item{Fiedler's result generalizes to \green{weighted} graphs (matrices) as well}

\item{Define \blue{Weighted Laplacian} as
\begin{displaymath}
\begin{array}{l}
\bar{L}(i,j) = -|A(i,j)|~~~~~~~~~~~
{\bf if}~i \neq j \\
\bar{L}(i,i) = \sum_j{|A(i,j)|}
\end{array}
\end{displaymath}
}
\item{The eigenvector that corresponds to the smallest non-zero eigenvalue 
of the weighted Laplacian \red{minimizes}
\begin{displaymath}
x^T\bar{L}x = \sum_{i,j}{|A(i,j)|(x(i)-x(j))^2},
\end{displaymath}
}
\item{Observe that $x^TLx$ is closely related to $\sum_{i,j}{|A(i,j)|}-w_b(A)$, with
proper quantization}

\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Weighted Bandwidth Reduction}

\begin{itemize}

\bigitem{For \green{large enough $\alpha$}, use Weighted Spectral Ordering as a heuristic
to minimize $\alpha$-bandwidth}
{\item Find $\hat{x}$, the eigenvector corresponding to smallest non-zero eigenvalue 
       of $\bar{L}$
 \item Find permutation $\Pi=\{i_1,i_2,...,i_n:~\text{if}~j<k,~x(i_j) < x(i_k)\}$
 \item Reorder rows and columns of $A$ accordingly to obtain $\bar{A}=A(\Pi,\Pi)$
} 

\item{Observe that the heavy entries of the reordered matrix, $\bar{A}$ are close
to its diagonal, \ie, $\bar{A}$ has a smaller $\alpha$-bandwidth compared to $A$}

\item{\red{Drop} all entries that are outside $\alpha$-bandwidth of $\bar{A}$}
\begin{displaymath}
\tilde{A} = \{\tilde{A}(i,j):\tilde{A}(i,j)=\bar{A}(i,j)~\text{if}~|i-j|\leq 
              BW_\alpha(\bar{A}),~0~\text{else}\} 
\end{displaymath}

\item{Use $\tilde{A}$ as a \blue{banded preconditioner} to solve the system $A$}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Experimental Results}

\begin{itemize}
\bigitem{Matrices gathered from \green{UF Sparse Matrix Collection}}
{
 \item All software is implemented in \green{Fortran}
 \item Sequential timings were done on a clovertown machine
 \item \green{BICGSTAB} is used as the iterative solver
 \item All matrices are first reordered using \green{MC64}, to move heaviest
       entries to the diagonal
}



\stitle{Experimental Results}
\bigitem{Application of \blue{Weighted Spectral Ordering (WSO)}}
{
\item Reorder $|A|+|A^T|$ using \blue{MC73}
\item Find the bandwidth that encapsulates \blue{99\%} of overall matrix norm ($\alpha=0.99$)
\item Drop entries that fall out of this bandwidth to obtain the \blue{Weighted Spectral
Preconditioner}
}



\stitle{Experimental Results}
\bigitem{Comparison with \red{no preconditioner} and \red{ILU}}
{
\item \red{ILUT} \mcite{Saad}{1994} is used as a basis for comparison.
\item Fill-in is set to match the storage required for dense storage
of WSO's required bandwidth (\red{ILUBW}) 
\item Drop tolerance is set to $10^{-1}$ (\red{ILU1}), $10^{-3}$ (\red{ILU3})
}
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{epb0 Matrix}

\begin{itemize}
\bigitem{Plate-fin heat exchanger w/ simple model}
{\item $1794\times1794$, $7764$ non-zeros
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/epb0.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/epb0_p.eps, width=0.475\textwidth}\\
Bandwidth: $3587$  & Bandwidth: $76$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $19$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{epb0 Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/epb0_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/epb0_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB converges in 14 iterations with WSO preconditioner}
{\item No convergence after 300 iterations with no preconditioner or ILU with
drop tolerance $10^{-1}$
}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{ASIC\_680k Matrix}

\begin{itemize}
\bigitem{Sandia, Xyce circuit simulation matrix (stripped)}
{\item $682862\times682862$, $2638997$ non-zeros
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/ASIC_680k.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/ASIC_680k_p.eps, width=0.475\textwidth}\\
Bandwidth: $1328578$  & Bandwidth: $1328726$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $5$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{ASIC\_680k Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/ASIC_680k_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/ASIC_680k_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB converges in 9 iterations with WSO preconditioner}
{
\item ILU with fill-in equivalent to bandwith of WSO preconditioner
converges faster (4 iterations), but factorization takes too much time
\item ILU factorization unsucessful for drop tolerance $10^{-1}$,$10^{-3}$ 
\item No convergence after 300 iterations with no preconditioner 
}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{lhr01 Matrix}

\begin{itemize}
\bigitem{Light hydrocarbon recovery}
{\item $14777\times14777$, $18427$ non-zeros
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/lhr01.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/lhr01_p.eps, width=0.475\textwidth}\\
Bandwidth: $2537$  & Bandwidth: $1341$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $601$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{lhr01 Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/lhr01_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/lhr01_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB converges in only 4 iterations with WSO preconditioner!}{
 \item No convergence with no preconditioner in 300 iterations
 \item ILU factorization unsucessful for all variants
}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{west0479 Matrix}

\begin{itemize}
\bigitem{U8 stage column section, all sections rigorous (chemical engineering)}
{\item $479\times479$, $1888$ non-zeros
 \item Known as a \red{horror matrix}
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/west0479.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/west0479_p.eps, width=0.475\textwidth}\\
Bandwidth: $726$  & Bandwidth: $872$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $221$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{west0479 Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/west0479_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/west0479_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB converges in 293 iterations with WSO preconditioner}{
 \item No convergence with no preconditioner in 300 iterations
 \item ILU factorization unsucessful for all variants
}
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{fp Matrix}

\begin{itemize}
\bigitem{2-D Fokker Planck equation, electron dynamics in external field}
{\item $7548\times7548$, $834222$ non-zeros
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/fp.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/fp_p.eps, width=0.475\textwidth}\\
Bandwidth: $15086$  & Bandwidth: $12368$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $5$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{fp Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/fp_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/fp_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB converges in 2 iterations with WSO preconditioner, as well as
 ILU with $10^{-3}$ drop tolerance, $2$ fill-in}
{\item Convergence in 7 iterations with no preconditioner
 \item Reordering, factorization take too much time as compared to
preconditioner savings
}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{matrix\_9 Matrix}

\begin{itemize}
\bigitem{Semiconductor device problem}
{\item $103430\times103430$, $1205518$ non-zeros
}
\end{itemize}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Original Matrix} &  \green{Reordered Matrix}\\
 \epsfig{file=figures/matrix_9.eps, width=0.475\textwidth} & 
 \epsfig{file=figures/matrix_9_p.eps, width=0.475\textwidth}\\
Bandwidth: $107410$  & Bandwidth: $198951$\\
 & \\
\multicolumn{2}{l}{$99\%$ of matrix norm lies within a bandwidth of $10673$ after WSO}\\
\end{tabular*}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{fp Results}

\begin{center}
\begin{tabular*}{\textwidth}%
     {@{\extracolsep{\fill}}cc}
\green{Residual Plot} & \green{Runtime Performance}\\
 \epsfig{file=figures/matrix_9_residue.eps,width=0.475\textwidth} &
 \epsfig{file=figures/matrix_9_performance.eps, width=0.475\textwidth}\\
\end{tabular*}
\end{center}

\begin{itemize}
\bigitem{BICGSTAB runs out of memory with WSO preconditioner}
{ \item BICGSTAB does not converge with no preconditioner
 \item ILU with $10^{-3}$ fill-in tolerance converges in 16 iterations
}
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\stitle{Summary}

\begin{center}
\begin{tabular}{lccccc}
\hline
& \multicolumn{5}{c}{\red{Preconditioner}}\\
\red{Matrix} & \blue{None} & \blue{ILUBW} & \blue{ILU1} & \blue{ILU3} & \blue{WSO}\\
\hline
\green{epb0} & $>0.019$ & $0.018$ & $>0.029$ & $0.059$ & $0.009$ \\
              & $>300$    & $33$     & $>300$    & $135$    & $14$ \\
\hline
\green{ASIC\_680k} & $>38.2$ & $329.8$ & $\infty$ & $\infty$ & $10.1$\\
                   & $>300$ & $4$ & $\infty$ & $\infty$ &  $9$\\ 
\hline
\green{lhr01} & $>0.02$ & $\infty$ & $\infty$ & $\infty$ & $0.11$\\
              & $>300$ &  $\infty$ & $\infty$ & $\infty$ & $4$\\
\hline
\green{west0479} & $>0.009$ & $\infty$ & $\infty$ & $\infty$ & $0.094$\\
                 & $>300$ & $\infty$ & $\infty$ & $\infty$ & $293$\\
\hline
\green{fp} & $0.05$ & $21.97$ & $\infty$ & $2.09$ & $0.68$\\
           & $7$ & $2$ & $\infty$ & $2$ & $2$ \\
\hline
\green{matrix\_9} & $>8.2$ & $39.2$ & $\infty$ & $3.5$ & $\infty$\\
                  & $>300$ & $91$ & $\infty$ & $16$ & $\infty$\\
\hline
\end{tabular}
\end{center}

\vspace{0.2in}
\centerline{Total runtime (reordering+factorization+bicgstab) is reported in seconds.}
\vspace{0.1in}
\centerline{Number of iterations are reported on the row below.}


\stitle{Remarks}

\begin{itemize}
\item
Banded preconditioners with suitable reordering techniques can be
very powerful for diverse classes of applications.
\item
Banded preconditioners typically yield much better CPU performance and
parallel performance.
\item
Due to memory reuse associated with dense kernels, they are well-suited
to conventional chip multiprocessor architectures.
\end{itemize}

\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\stitle{References}

\begin{thebibliography}{1}

%\bibitem{chan94}
%P.~K. Chan, M.~D.~F. Schlag, and J.Y. Zien.
%\newblock Spectral k-way ratio-cut partitioning and clustering.
%\newblock {\em IEEE Transactions on Computer-Aided Design of Integrated
%  Circuits and Systems}, 13(9):1088--96, 1994.

\bibitem{cuthill69}
E.~Cuthill and J.~McKee.
\newblock Reducing the bandwidth of sparse symmetric matrices.
\newblock In {\em In Proc. 24th Nat. Conf. ACM}, pages 157--172, 1969.

\bibitem{fiedler73}
M.~Fiedler.
\newblock Algebraic connectivity of graphs.
\newblock {\em Czechoslovak Mathematical Journal}, 1973.

\bibitem{saad94}
Y.~Saad. I
\newblock {ILUT}: A dual threshold incomplete {ILU} factorization. 
\newblock {\em Numerical Linear Algebra with Applications}, 1:387--402, 
1994.

\bibitem{barnard95}
S.~T. Barnard, A. Pothen, and H. Simon.
\newblock A spectral algorithm for envelope reduction of sparse matrices.
\newblock {\em Numerical Linear Algebra with Applications}, 2(4):317--334,
  1995.

\bibitem{kruyt96}
N.~P.~Kruyt.
\newblock A conjugate gradient method for the spectral partitioning of 
graphs.
\newblock {\em Parallel Computing}, 22(11):1493--1502, 1996.



%\bibitem{koren03}
%Yehuda Koren.
%\newblock On spectral graph drawing.
%\newblock In {\em Computing and Combinatorics: 9th Annual International
%  Conference, COCOON 2003}, pages 496 -- 508, 2003.

\end{thebibliography}

\end{document}
