\documentstyle[12pt]{article}
\title{A tricyclic tridiagonal equation solver}
\author{David S. Dodson\thanks{Mathematical
Software Group, Convex Computer Corp.,
P.O. Box 833851, Dallas, Texas 75083-3851.}
\and
Stewart A. Levin\thanks{Mobil Research and Development Corp.,
P.O. Box 819047, Dallas, Texas 75381-9047.}}
\date{\sl draft \rm \today}
\begin{document}
\maketitle
\begin{abstract}
We exhibit an improved method for solving tridiagonal equations on
Convex
computers.  The method employs cyclic reduction by powers of three.
The advantages over conventional power-of-two reduction are several.
The number of divisions per element is cut in half, while the number
of multiplications and additions remain almost exactly the same.
Memory bank conflicts are minimized because all vector strides are
powers of three, i.e. odd.  The number of memory accesses is reduced by
a quarter.  This is especially important for the Convex where cyclic
reduction is memory-limited.  Lastly, the algorithm supports modification
in a recently published manner that maximizes the use of full vector register
segments throughout the reduction.
(R. Reuter, 1988, Solving tridiagonal
systems of linear equations on the IBM 3090 VF, Parallel Comput. 8, 371-376.)
\end{abstract}
\section{Introduction}

\hspace*{\parindent}
Cyclic reduction in various forms is used routinely to improve the speed
of solution of tridiagonal systems on vector and parallel computers.
Variables are grouped into two disjoint index sets and simultaneously
eliminated from the
corresponding tridiagonal equations to yield a descending sequence of
progressively smaller tridiagonal systems to solve.
Because this reduction process and corresponding back substitution are
vectorizable, this algorithm is beneficial despite the doubled arithmetic
load it incurs when compared to the conventional serial algorithm,
e.g. SGTSL from LINPACK~\cite{JJD:Lin}.

In the past, the index sets have been the even and odd variables, resulting in
a halving of the system size at each reduction step.  Problems that researchers
have wrestled with are memory bank conflicts due to the even strides between
vector elements and inefficient use of vector registers due to end effects and
short vectors.

In this report we explore a different choice of index sets,
specifically every third variable in one set and the rest in the other,
and show that it leads to a faster tridiagonal solution algorithm on a
Convex computer.

\section{Tricyclic reduction}
\subsection{Algorithm}
\hspace*{\parindent}
We describe the in-place version of power-of-three cyclic reduction.
Denote the
subdiagonal elements $a_j$, $j=1,...,n-1$, the diagonal elements $b_j$,
$j=0,...,n-1$, the superdiagonal elements $c_j$, $j=0,...,n-2$, and the
input right hand side $d_j$, $j=0,...,n-1$.  Pictorially,
\begin{displaymath}
\begin{array}{lllllllll}
a_{3j-2}&b_{3j-2}&c_{3j-2}& & & & &\;&d_{3j-2} \\
 &a_{3j-1}&b_{3j-1}&c_{3j-1}& & & &&d_{3j-1} \\
 & &a_{3j}&b_{3j}&c_{3j}& & &&d_{3j} \\
 & & &a_{3j+1}&b_{3j+1}&c_{3j+1}& &&d_{3j+1} \\
 & & & &a_{3j+2}&b_{3j+2}&c_{3j+2}&&d_{3j+2}
\end{array}
\end{displaymath}
are the matrix entries participating in the reduction steps.

The first step is to invert the submatrices
\begin{displaymath}
\left(
\begin{array}{ll}
b_{3j+1}&c_{3j+1}\\
a_{3j+2}&b_{3j+2}
\end{array}
\right)
\end{displaymath}
which we do directly using Cramer's Rule to produce
\begin{displaymath}
\left(
\begin{array}{rr}
b_{3j+2}'&-c_{3j+1}'\\
-a_{3j+2}'&b_{3j+1}'
\end{array}
\right)
\; = \;
\frac{1}{b_{3j+1}b_{3j+2}-c_{3j+1}a_{3j+2}}
\left(
\begin{array}{rr}
b_{3j+2}&-c_{3j+1} \\
-a_{3j+2}&b_{3j+1}
\end{array}
\right)
\;.
\end{displaymath}
Here the primed entries may be stored in place of the original unprimed
matrix entries.
This matrix inverse is then multiplied into the corresponding rows
of the original system to bring it into the form
\begin{displaymath}
\begin{array}{lllllllll}
\phantom{-}a_{3j-2}'&1&0&-c_{3j-2}''& & & &\;&d_{3j-2}' \\
-a_{3j-1}''&0&1&\phantom{-}c_{3j-1}'& & & &&d_{3j-1}' \\
 & &a_{3j}&\phantom{-}b_{3j}&c_{3j}& & &&d_{3j} \\
 & & &\phantom{-}a_{3j+1}'&1&0&-c_{3j+1}''&&d_{3j+1}' \\
 & & &-a_{3j+2}''&0&1&\phantom{-}c_{3j+2}'&&d_{3j+2}'
\end{array}
\end{displaymath}
where
\begin{displaymath}
\begin{array}{lcl}
d_{3j+1}'&=&b_{3j+2}'d_{3j+1}-c_{3j+1}'d_{3j+2}\\
d_{3j+2}'&=&b_{3j+1}'d_{3j+2}-a_{3j+2}'d_{3j+1}\\
a_{3j+2}''&=&a_{3j+2}'a_{3j+1} \\
a_{3j+1}'&=&a_{3j+1}b_{3j+2}' \\
c_{3j+1}''&=&c_{3j+1}'c_{3j+2} \\
c_{3j+2}'&=&c_{3j+2}b_{3j+1}'
\end{array}
\end{displaymath}
replace the previous values of the indicated variables.

To complete the reduction step
we eliminate $a_{3j}$ and $c_{3j}$
to produce
\begin{displaymath}
\begin{array}{lllllllll}
\phantom{-}a_{3j-2}'&1&0&-c_{3j-2}''& & & &\;&d_{3j-2}' \\
-a_{3j-1}''&0&1&\phantom{-}c_{3j-1}'& & & &&d_{3j-1}' \\
\phantom{-}a_{3j}' &0&0&\phantom{-}b_{3j}'&0&0&\phantom{-}c_{3j}'&&d_{3j}' \\
 & & &\phantom{-}a_{3j+1}'&1&0&-c_{3j+1}''&&d_{3j+1}' \\
 & & &-a_{3j+2}''&0&1&\phantom{-}c_{3j+2}'&&d_{3j+2}'
\end{array}
\end{displaymath}
with
\begin{displaymath}
\begin{array}{lcl}
b_{3j}'&=&b_{3j}-c_{3j}a_{3j+1}'-a_{3j}c_{3j-1}'\\
d_{3j}'&=&d_{3j}-c_{3j}d_{3j+1}'-a_{3j}d_{3j-1}'\\
a_{3j}'&=&a_{3j}a_{3j-1}'' \\
c_{3j}'&=&c_{3j}c_{3j+1}''
\end{array}
\end{displaymath}
again replacing the corresponding unprimed quantities.

We recognize now that the equations in $a_{3j}'$, $b_{3j}'$,
$c_{3j}'$ and $d_{3j}'$ now form a tridiagonal system one third
the size of the orignal system.  Thus we reapply the tricyclic
reduction to this reduced system.

Having solved the reduced system for the components $d_{3j}''$
of the solution vector, we back substitute according the the formulas
\begin{displaymath}
\begin{array}{lcl}
d_{3j+1}''&=&d_{3j+1}'-a_{3j+1}'d_{3j}''+c_{3j+1}''d_{3j+3}'' \\
d_{3j+2}''&=&d_{3j+2}'+a_{3j+2}''d_{3j}''-c_{3j+2}'d_{3j+3}''
\end{array}
\end{displaymath}
to obtain the remaining components of the solution.

\subsection{Operation Counts}

\hspace*{\parindent}
To analyze this algorithm, we tally separately the operations
{\sl load, store, add, multiply, divide,} and {\sl shift}.  Here {\sl shift\/}
is the operation of moving the components of a vector one element
left or right.  On the Convex, this may be accomplished by (assembly
language) register-to-register vector
compression or by storing to memory and then 
reloading with an offset.
Forming the submatrix inverse requires 4 {\sl load}'s,
6 {\sl mult}'s, 1 {\sl add}, and 1 {\sl div}.  Applying the submatrix
inverse takes 4 {\sl load}'s, 8 {\sl mult}'s, 2 {\sl add}'s, and 6
{\sl store}'s.
Finishing the reduction uses 4 {\sl load}'s, 3 {\sl shift}'s,
6 {\sl mult}'s, 4 {\sl add}'s, and 4 {\sl store}'s.  Finally, the
backsubstitution uses 7 {\sl load}'s, 1 {\sl shift}, 4 {\sl mult}'s,
4 {\sl add}'s, and 2 {\sl store}'s.
\par
The number of reductions and backsubstitutions is approximately
\begin{displaymath}
n \sum_{j=1}^{\infty} \frac{1}{3^j} = \frac{n}{2} \hbox{\bf \ \ \ ,}
\end{displaymath}
leading to an average operation count per equation of
\begin{displaymath}
\hbox{$1\over2$}
(12 \hbox{\sl \ load} + 3 \hbox{\sl \ shift} + 20 \hbox{\sl \ mult}
+ 7 \hbox{\sl \ add} + 1 \hbox{\sl \ div} + 10 \hbox{\sl \ store} )
\end{displaymath}
for reduction and
\begin{displaymath}
\hbox{$1\over2$}
(7 \hbox{\sl \ load} + 1 \hbox{\sl \ shift} + 4 \hbox{\sl \ mult}
+ 4 \hbox{\sl \ add} + 0 \hbox{\sl \ div} + 2 \hbox{\sl \ store} )
\end{displaymath}
for backsubstitution.

The comparable operation counts for power-of-two cyclic reduction are
\begin{displaymath}
8 \hbox{\sl \ load} + 3 \hbox{\sl \ shift} + 9 \hbox{\sl \ mult}
+ 4 \hbox{\sl \ add} + 1 \hbox{\sl \ div} + 7 \hbox{\sl \ store}
\end{displaymath}
for reduction and
\begin{displaymath}
4 \hbox{\sl \ load} + 1 \hbox{\sl \ shift} + 2 \hbox{\sl \ mult}
+ 2 \hbox{\sl \ add} + 0 \hbox{\sl \ div} + 1 \hbox{\sl \ store}
\end{displaymath}
for backsubstitution.

Figure~\ref{Convex:OPS} summarizes these operation counts graphically
\begin{figure}
\vspace{3.5in}
\caption{ \label{Convex:OPS}Average per-element operation counts
for power-of-two versus
power-of-three tridiagonal cyclic reduction algorithms on the Convex.  
Tricyclic reduction is always as good or better than the power-of-two
method.}
\end{figure}
for the Convex.  Note that we have conservatively
equated each division to 2.25 multiplications \cite{Convex:ARCH}.
It is apparent that on a Convex the power-of-three algorithm is
never inferior to the power-of-two method in this measure of
performance.


\subsection{Analysis}

\hspace*{\parindent}
To better appreciate the relative merits of power-of-three
versus power-of-two cyclic
reduction, we start by looking at some limiting cases.

\subsubsection{Fast memory}

\hspace*{\parindent}
In this setting we neglect memory accesses and find that the power-of-three
method uses $ {1\over2} (24 \hbox{\sl \ mult} + 11 \hbox{\sl \ add}
+ 1 \hbox{\sl \ div} )$ operations per element and the power-of-two
method uses $ 11 \hbox{\sl \ mult} + 6 \hbox{\sl \ add} + 1 \hbox{\sl \ div}$
operations per element.  If, as typical on vector computers, a
division is 2 to 3 times slower than a multiplication,
then the power-of-three has slightly less arithmetic than
power-of-two.  An analogous finding for FFT's
was made in~\cite{CT:FFT} where it was observed that, arithmetically,
the power-of-three FFT is about 6 percent more efficient than the power-of-two
version.

\subsubsection{Fast arithmetic}

\hspace*{\parindent}
When arithmetic is much faster than memory accesses times,
we neglect computation time and find that
power-of-three reduction uses $ {1\over 2} ( 19 \hbox{\sl \ load} +
4 \hbox{\sl \ shift} + 12 \hbox{\sl \ store} )$\/
operations and power-of-two uses
$ 12 \hbox{\sl \ load} +3 \hbox{\sl \ shift} + 8 \hbox{\sl \ store}$\/
operations per equation.
In this setting the power-of-three approach saves 25 percent
over power-of-two cyclic reduction, regardless of whether the shift is
implemented by memory operations or by register-to-register vector
compression.  In our fastest Fortran versions of these cyclic reduction
algorithms all but one of the {\sl shift}'s translate into two
memory operations.

For vector computers there is an additional factor that affects memory
access speeds.  Invariably, the memory is divided into interleaved banks
in order to permit faster reference to contiguous vectors.  Accessing
vectors with an element-to-element stride differing from unity can be
significantly slower if the spacing causes particular memory banks
to be addressed more often.  In practice this usually penalizes
access to vectors with strides divisible by a power of two, exactly
the case in power-of-two cyclic reduction.  Thus the true cost of
memory accesses in power-of-two reduction may be significantly higher,
tipping the scales more heavily in favor of our tricyclic reduction.

For parallel, distributed-memory computers, the memory accesses involve
processor-to-processor communication.
This is generally much slower than arithmetic.
In practice, though, communication is usually fastest for elements
spaced a power-of-two apart and then the scales tip towards power-of-two
cyclic reduction.  An example of this is the hypercube multiprocessor
architecture wherein elements spaced by a power-of-two are at most
two hops apart in the physical processor grid.  On the other hand, with
a processor ring the amount of communication per step in the reduction
is approximately constant whence the power-of-three is better by the
ratio $\log 3 / \log 2 $ because there are fewer reduction steps.

\subsubsection{Early termination}

\hspace*{\parindent}
One of the more useful, though infrequently mentioned, advantages of cyclic
reduction is the rapid decrease in the off-diagonal terms of the reduced
systems when the original system is diagonally dominant.  By predicting
or monitoring this decrease one can terminate the reduction stage early when
the off-diagonal components become negligible.

To attain the same off-diagonal reduction, the power-of-three method will
use two thirds of the number of reduction steps needed by the power-of-two
version.  This slightly increases the ratio of the number of elements
processed by the two methods from its asymptotic value of 1:2 improving
the competitiveness of power-of-two reduction a little.  More importantly,
the power-of-two has more places it can stop and thus might shave one or two
iterations from the asymptotic 3/2's of the iteration count of power-of-three.
This is potentially a significant savings for large tridiagonal systems.
In some cases this savings may be partially offset by changing the last
reduction by three into a reduction by two.

\subsubsection{Cache}

\hspace*{\parindent}
The first-generation Convex C1 required that vector
operations on noncontiguous elements be first staged into cache memory.
(The C200 Series bypasses cache for all vector operations.)  Since cyclic
reduction, whether by powers of two or three, involves noncontiguous
vector access, cache usage becomes a concern.  One's aim is to minimize
the number of passes over the data so that cache is ``missed'' as few
times as possible.  Because it uses fewer passes over the data, tricyclic
reduction has an advantage.  This advantage is nullified, however, if the
partial reduction in~\cite{R:Cache} is employed.  In this variation,
a partial reduction plus a shuffle is used to reduce the tridiagonal system
by one (or two) vector register lengths.  The result is that the data is
traversed only twice: once forward during reduction and once backwards
during resubstitution.  Cache misses are minimized.  Additionally,
vector register usage is maximized.  We are not aware of any comparative
round-off study for this variation.

\subsubsection{Parallelism}

\hspace*{\parindent}
Cyclic reduction can be parallelized to a good degree on the Convex C200 Series
by the conventional method of distributing the iterations of the strip-mine
sections of the relevant DO-loops among the available processors.
The {\sl shift\/} operations
cause some difficulty for the calculations at each end of a vector segment
in a reduction step because some of the required data are assigned
to a cooperating processor.  This can be handled by synchronized memory
references or by duplicating the end calculations in both relevant
processors.  The latter is advantageous when the {\sl shift\/} operation
is handled by vector-to-vector compression.

Memory bank conflicts are a significant concern for parallel operations.
The more simultaneous power-of-two memory accesses occurring in
parallel, the slower each becomes.  For this reason, tricyclic reduction
is to be preferred because memory bank conflicts are resolved by a systematic
increase in startup time for memory operations rather than a systematic
increase in per-element memory transfer time.

\subsection{Convex timings}

\hspace*{\parindent}
The analysis above has indicated that tricyclic reduction is a choice method
for solving tridiagonal equations on Convex computers.
We have programmed
the algorithm in both Fortran (compiler version 6.0) and Convex assembly
language
for the C200 Series.
Both the Fortran and the assembler versions vectorize and parallelize.
We have also compiled the Fortran version with parallelization disabled.
Lastly, for comparative timings, we have programmed power-of-two
cyclic reduction in FOrtran, using auxiliary
space to compress intermediate results.
This we label
{\it sericr\/} following standard terminology \cite{HJ:Tri}.
Figure~\ref{fig:CPU} displays
CPU times for these four algorithms.
\begin{figure}
\vspace{3.5in}
\caption{ \label{fig:CPU}CPU times for solving tridiagonal equations using
power-of-three and power-of-two cyclic reductions.   Because
it uses fewer memory operations, tricyclic reduction
is around 25 percent faster than power-of-two cyclic reduction
({\it sericr\/}) on the Convex.}
\end{figure}
As predicted, Fortran savings are nearly 25 percent for the
power-of-three algorithm over the power-of-two version.  Assembly
programming buys an additional ten percent efficiency over Fortran.
The small gain achieved with assembler language demonstrates that the
algorithm can be expressed in Fortran without hiding the available
paralellism, that the method is matched well to the Convex architecture,
and that the Convex Fortran compiler's vectorization, parallelization,
and instruction scheduling are quite mature.
Figure~\ref{fig:WALL} displays parallel performance of the tricyclic algorithm.
\begin{figure}
\vspace{3.5in}
\caption{\label{fig:WALL}Parallel speedups of wall-clock
times over CPU times for
power-of-three and power-of-two cyclic reduction on a four
processor Convex C240.
Tricyclic reduction is better because it
has fewer loops with more work per loop than power-of-two
reduction.}
\end{figure}
For large
systems, parallel speedup levels off at around 3.5 out of a possible 4.
In comparison, the power-of-two speedup peaks at a lower value of
about 3.1.  This reflects the fact that the power-of-two algorithm
uses more DO loops with less work per loop and so synchronization
occurs more frequently with concomittant loss of parallelism.

\section{History and Related Work}

\hspace*{\parindent}
To our knowledge, tricyclic reduction first appeared in 1980 when it was
employed by one of the authors to speed up an
implicit finite-difference program running on an IBM 3838 array processor.
It was later resurrected for a Cray-1
assembly language exercise in 1986.
About the same time,
Phuong Vu of Cray Research independently devised
the algorithm for use on the Cray-2 (pers. comm.)

Tricyclic reduction is closely related to algorithm GECR
(\cite{J:Ens}, \cite{W:GECR})
when specialized to partitions of size 3.
Indeed, neglecting round-off variations, a single reduction step
in both procedures will produce identical results.
There are two differences, however.
First, Cramer's Rule is used in tricyclic reduction as
compared to Gaussian elimination in GECR.
This cuts the number of divide operations in half compared to GECR.
Second, GECR does not recursively reapply the reduction to the
smaller resulting system as it is assumed that a proper choice of partition size
renders further parallel reduction unnecessary.  For the Convex,
GECR would subdivide the equations so that the reduced size would
be as close as possible to the vector register length of 128.  At that
point another algorithm, cyclic reduction, would be applied
to this reduced system.  We do not yet  know whether, in practice, GECR is
faster or slower than tricyclic reduction for the Convex.

\section{Conclusions}

\hspace*{\parindent}
Cyclic reduction by powers of three is significantly faster than cyclic
reduction
by powers of two for solving tridiagonal equations on the Convex.
The algorithm achieves its better performance primarily 
by reducing the number
of memory accesses needed to solve the system.
Additionally, it uses memory more efficiently because vector strides
are always odd.  Lastly, its use of half the division operations in fewer
loops with more arithmetic per loop produces a somewhat better
arithmetic balance and compiler utilization.    These performance gains
have been also proven on the Cray-2 by independent invention of the
algorithm.

\begin{thebibliography}{9}

\bibitem{Convex:ARCH} {\em CONVEX Architecture Reference (C100 Series,
C200 Series)}, 4th ed., 1989, CONVEX Computer Corp. Doc. 081-000050-202,
Richardson, TX.

\bibitem{CT:FFT} Cooley, J.W. and Tukey, J.W., 1965, An algorithm for the
machine calculation of complex Fourier series, Math. Comput. 19, 297-301.
Reprinted in: {\em Digital Signal Processing}, Rabiner, L.R. and Rader, C.M,
editors, IEEE Press, New York, 1972.

\bibitem{JJD:Lin} Dongarra, J.J., Moler, C.B., Bunc, J.R., and
Stewart , G.W., 1979,
{\em LINPACK: Users' Guide}, SIAM, Philadelphia, 320 p.

\bibitem{HJ:Tri} Hockney, R.W. and Jesshope, C.R., 1981,
{\em Parallel Computers}, Adam Hilger, Bristol, 423 p.

\bibitem{J:Ens} Johnsson, S.L., 1987, Solving tridiagonal systems on
ensemble architectures, SIAM J. Sci. Stat. Comput. 8, 354-392.

\bibitem{R:Cache} Reuter, R., 1988, Solving tridiagonal systems of linear
equations on the IBM 3090 VF, Parallel Computing 8, 371-376.

\bibitem{W:GECR} Wang, H.H., 1981,
A parallel method for tridiagonal equations,
ACM TOMS 7, 170-183.

\end{thebibliography} 
\end{document}
