\Section{Preconditioned Conjugate Gradient}
\label{Section:pcg}

\begin{fig}{figure:pcg}{Inner loop of preconditioned conjugate gradient in Id.}
\begin{center}
\begin{minipage}{3in}
\footnotesize
\begin{verbatim}
{while (ip_3d r r) > epsilon do
   gamma   = ip_3d r xi;
   Ab      = seven_pt_stencil b;
   alpha   = gamma / (ip_3d b Ab);
   next p  = daxpy_3d p alpha b;
   next r  = daxpy_3d r (-alpha) Ab;
   next xi = fb_preconditioner (next r);
   beta    = (ip_3d (next r) (next xi)) / 
                gamma;
   next b  = daxpy_3d (next xi) beta b;
  finally p}
\end{verbatim}
\end{minipage}
\end{center}
\end{fig}

The preconditioned conjugate gradient algorithm \cite{templates} is the
core of many computational fluid dynamics applications and its regular
structure is similar to many other scientific applications.  It solves
linear systems using an indirect (i.e.\ iterative) method.  The
preconditioning step of PCG is used to speed convergence, and varies
depending upon the characteristics of the linear systems being solved.  The
algorithm itself is composed of only four basic operations; in the case of
the GCM ocean modelling code, the operations are implemented as follows:

\begin{itemize}
\item {\bf 3D inner product} -- this is the computation described in the
previous Section.
\item {\bf 3D daxpy} --  $aX+Y$, where $X$ and $Y$ are 3D vectors and $a$ is a
scalar constant.  3D daxpy requires a floating point add and multiply per
element.
\item {\bf 7-point stencil} -- the stencil takes a weighted average of the
local element and the six surrounding elements.  7-point stencil requires
six floating point additions and seven floating point multiplications per
element.
\item {\bf forward-backward substitution} -- the preconditioner we chose
for our PCG is a forward-backward substitution of a tridiagonal linear
system performed in each column of the ocean.  The forward-backward
substitution requires about five floating point operations per element.
\end{itemize}

Each of these operations can be implemented as a triply-nested loop.  The
loops are simple, structured, and have no dependencies between iterations.
Results from some of these operations feed into others, and because of
Monsoon's fine-grained word-level synchronization on memory, operations may
proceed before their input arrays are completely defined.

The structure of the state arrays directly reflects the geometry of the
ocean being modeled.  Certain elements of the arrays represent land state,
and certain elements represent water state.  We are not concerned with
elements on land because they do not change with time.  We can structure
most of the computation so that all computation is uniform over all
elements of the array, or just attempt to do computation on water.  In some
cases, land may take up to 60-70\% of the geometry.

We consider only computation performed on water to be required, and for a
given geometry, the number of {\it required floating point} operations is
constant for this algorithm.  Computations performed on land are extraneous
and overhead.  In the next Section, we consider optimization by computing
only on water elements.

\SubSection{Optimizing out computation on land}

\begin{fig}{figure:gcm-comparison}{Comparison of the three versions of 
PCG for the $168 \times 90 \times 4$ Pacific ocean, 60\% water.  ($ib=5$,
$jb=3$, 2 processor Monsoon)}
\illustrationx{2.9in}{figures/three-versions.ps}
\end{fig}

\begin{fig}{figure:operation-comparison}{The components of the hybrid PCG,
the PCG itself, and the whole GCM ocean modeling code for the $168 \times
90 \times 4$ Pacific ocean.  ($ib=5$, $jb=3$, 2 processor Monsoon)}
\illustrationx{3.1in}{figures/operation-comparison.ps}
\end{fig}

By eliminating computation on land, the maximum speedup we could expect
over an implementation that computes on both land and water is the ratio of
land and water elements over water elements.  For most real geometries,
this is a maximum speedup of three or four, and for the Pacific ocean
geometry, the maximum speedup is less than a factor of two.

We also incur overhead when we attempt to compute only on water, because we
must use sequential loops in the innermost loop (which computes over
columns in the ocean) as compared to completely unrolled loops in the case
where we compute on both land and water.  Figure
\ref{figure:gcm-comparison} shows that the version which computes on land
and water actually is faster than the version which computes only on water,
even though fewer floating point operations are performed in the water-only
version.  For this type of application on this type of architecture,
reducing floating point operations is less important than simplifying the
structure of the code to reduce loop overheads.

In an attempt to take some advantage of both simplified loop structure and
eliminating computation on land elements, we also developed a {\it hybrid}
version of the PCG which only computes on the bounding hull of the water on
the surface of the geometry, while computing on entire columns of water
within the bounding hull.  In this way, we eliminate much of the
computation on land, while still allowing complete unrolling of the
innermost loop.  Figure \ref{figure:gcm-comparison} shows that this version
is slightly faster than the version which computes on both land and water.

The number of idle cycles in the graph of Figure
\ref{figure:gcm-comparison} is misleading because PCG is only one part
(albeit the most computationally intense part) of the entire GCM code.
When the PCG is incorporated into the larger code, parallelism from the
rest of the code can mask idle cycles in PCG.

\SubSection{Breakdown by operation}

Figure~\ref{figure:operation-comparison} shows the cycles per floating
point operations for the four operations of the hybrid version of PCG, the
PCG itself and the overall GCM code.  The constituent parts were timed by
extracting the code fragment for the component from the elliptic solver and
running it in isolation.  Note that the idle cycles for some of the
components are off the scale -- however, when the four components are
integrated together into the PCG, the idle cycles are masked by work from
other components.  Also, when the PCG itself is integrated into the the
overall GCM code, the instruction counts heavily reflect the ratios of the
PCG, but idle cycles decrease because of other work being performed which
masks those idle cycles when the PCG is run in isolation.  This masking of
idle cycles is one of the primary benefits of the multithreaded execution
model.

When we set $k$-bounds for loops, we must do it in the context of an entire
application, because as Figure ~\ref{figure:operation-comparison} shows,
processor utilization for codes with loop bounds set in isolation do not
reflect utilization for integrated codes.
