\section{The ocean model}
\label{section:ocean-model}

\begin{fig}{figure:simple-ocean}{Ocean states are represented as three
dimensional arrays.  Some elements of the arrays represent land zones and
some represent water zones.}
\illustrationx{2.9in}{figures/simple-ocean.eps}
\end{fig}

\begin{fig}{figure:pacific}{A significant percentage of real ocean
geometries may be land.  In this $171 \times 93 \times 4$ Pacific basin
geometry, 40.0\% of the zones represent land.  Land zone states remain
constant over time.}
\illustrationx{3.1in}{figures/ocean.eps}
\end{fig}

The model used in this study is the MIT ocean circulation model, referred
to in this paper as the ``General Circulation Model'' or GCM. It is being
actively used in oceanographic research, and daily production runs are
performed on a CM-5 and, more recently, on high-end symmetric
multi-processor systems.

The physical basis of the model and its continuous and discrete forms are
described briefly in the Appendix, and more detailed accounts can be found
in \cite{Marshall97a}, \cite{Marshall97b} and \cite{Hill95}.  The model is
based on the incompressible Navier-Stokes equations and can be used to
study the ocean from convective ($\sim$100m in the horizontal) up to global
scales.  The numerical procedure involves alternating prognostic and
diagnostic steps: prognostic to step forward in time the ocean currents and
thermodynamic variables, diagnostic to find the pressure field.  The
``pressure correction'' method is used to ensure that the evolving velocity
field remains non-divergent. A Poisson equation with Neumann boundary
conditions is solved to diagnose for the pressure field.  This latter step
is a computationally challenging one because it involves communication of
information across the whole model grid to the boundary, in a geometry as
complicated as that of an ocean basin.

\subsection{Ocean geometry and state}
\label{section:ocean-state}

\noindent The ocean is confined to a basin that can have a highly
irregular geometry.  A {\it geometry}, which defines the shape of the ocean
basin, is represented as a three-dimensional array of finite volumes that 
we call {\it zones}. The faces of the zones must be chosen to coincide with
coordinate surfaces, except when they abut a solid boundary (the coast or
bottom). At these solid boundaries the zones may be ``sculptured'' to fit
them to irregularities in the topography \cite{Marshall97b,adc97}.  
Each element
of the geometry array specifies whether the corresponding zone is land or
water. We will assume that the ocean has a rigid lid at the upper surface.  
The rigid lid
is a device to filter from the model rapidly propagating surface
gravity waves that would severely restrict the possible length of the model
time-step.  For some integrations a free surface is admitted, but handled 
implicitly as discussed in \cite{Marshall97b}.
The geometry array is constructed from the bathymetry -- a two-dimensional 
array of depths for each column of the ocean at a chosen
horizontal resolution.  A highly idealized geometry for a $4 \times 4$
ocean with $4$ layers of water is shown in
Figure~\ref{figure:simple-ocean}. The white zones are water; the gray zones
land.  A more realistic geometry is shown in
Figure~\ref{figure:pacific}, representing the Pacific basin at the
resolution of $1$ degree of latitude by $1$ degree of longitude; there
are $171$ zones along each line of parallel, $93$ zones along each line of
longitude and $4$ zones in the vertical.  In typical ocean modeling
applications the number of zones stacked on top of one another is between 1
and 100, generally fewer than the number employed in the horizontal.

Each zone has six faces, and these faces categorized according to orientation:
\begin{itemize}
\item {\it longfaces}: vertical faces running along longitudes
\item {\it parfaces}: vertical faces running along parallels
\item {\it horfaces}: horizontal faces
\end{itemize}
The faces of the zones remain fixed during the integration over time; they
are dependent only on the geometry of the ocean basin and the resolution at
which it is discretized.

The state of the ocean is described by the velocity $\vel = (u,v,w)$, the
pressure $p$ and the density $\rho$.  The density is a function of the
(potential) temperature $T$, the salinity $S$ and pressure $p$.  Zone
quantities -- $p$, $\rho$, $T$ and $S$ -- are defined at the center of the
zone and are implemented as {\it zone arrays}. Face quantities are defined
on the faces of the zones -- velocity components $\vel$ -- in {\it longface
arrays} ($u$), {\it parface arrays} ($v$) and {\it horface arrays} ($w$).
Zone and face arrays are three-dimensional arrays, and typically there are
a dozen such arrays in use at any one time.

In certain applications of the model (see \cite{Marshall97a}) appropriate
periodicity conditions must be employed.

\subsection{The numerical algorithm}
\label{section:main-loop}

The heart of the model is a main time-stepping loop that is repeated many
times.  Depending upon the dynamical process being studied each time step may 
represent several hours of time or just a few seconds; a typical
simulation entails thousands or millions of time steps. The main loop computes
the changes in the values of $\vel = (u,v,w)$, $p$, $T$ and $S$ by
applying the laws of classical mechanics and thermodynamics
to the fluid contained within the zones.
Discrete forms of the continuous equations are deduced by integrating them
over the zones and using Gauss's theorem.  From $p$, $T$ and $S$ the
density $\rho$ is computed.

During each time step, it is also necessary to compute the forces normal to 
the faces of the zones to yield, according to Newton's Laws, the 
acceleration of fluid parcels. This is done by evaluating the pressure-gradient $G$, inertial, Coriolis,
and applied stresses and frictional forces, as described in 
\cite{Marshall97a,Marshall97b}.

\begin{fig}{figure:gcm-high-level}{High-level specification of the GCM algorithm.}
\begin{tabbing}
a \= b \= Step 4: $\Tm$ \= \kill
\procname{GCM} () \\
\> $\velinit,\Pinit,\Tinit,\Sinit,\rinit$ = \procname{compute\_initial\_ocean\_states}() \\
\> \keyword{for} each time step $n$ \\
\> \> Step 1: $\Gveli$ \> = \procname{estimate\_Gv}($\veln$,$\Pn$,$\rn$) \\ 
\> \> Step 2: $\Pm$    \> = \procname{diagnose\_pressure}($\Pn$,$\veln$,$\Gveli$)\\
\> \> Step 3: $\velm$  \> = \procname{advance\_velocity}($\Pm$,$\vn$,$\Gveli$) \\
\> \> Step 4: $\Tm$    \> = \procname{advance\_temperature}($\Tn$,$\velm$) \\
\> \> Step 5: $\Sm$    \> = \procname{advance\_salinity}($\Sn$,$\velm$) \\
\> \> Step 6: $\rho^{n+1}$ \> = \procname{calculate\_density}($\Tm$,$\Sm$) \\
\> \keyword{return}($\velfinal,\Pfinal,\Tfinal,\Sfinal,\rfinal$)
\end{tabbing}
\end{fig}

A high-level description of GCM is shown in Figure
\ref{figure:gcm-high-level}.  Steps 1, 3, 4, 5 and 6 only require local
communication and nearest-neighbor communication to update ocean state
arrays, but Step 2 necessitates solving the following linear system,
a discrete form of the Poisson equation for the pressure:
\begin{equation}
\A \Pm = f(\Gvel,\veln) 
\label{equation:elliptic-eq}
\end{equation}
where $\A$ is a sparse $N \times N$ matrix. 

Our algorithm is designed to perform efficiently across the scales of
interest in the ocean, from non-hydrostatic phenomena on the smallest
scales, to the highly balanced hydrostatic flows on large scales.  The
approach, which involves separating the pressure field in to its
hydrostatic, surface and non-hydrostatic components, is set out in
\cite{Marshall97a}.  In the hydrostatic limit, Equation
(\ref{equation:elliptic-eq}) is a two-dimensional inversion for the surface
pressure, with $N = N_x \times N_y$.  The pressure in the interior of the
ocean is then obtained on intergation vertically using the hydrostatic
approximation. But in non-hydrostatic applications of the model $N = N_x
\times N_y \times N_z$ is the total number of zones, $p$ is a vector of
length $N$ carrying the pressure of each zone, and $f$ is the source
function.  In typical ocean modeling applications, $N_x$ and $N_y$ range
from 100 to 1,000; and $N_z$ from 1 to 100. Thus, $\A$ is potentially huge,
and special techniques must be used to invert Equation
(\ref{equation:elliptic-eq}) for $\Pm$.  Because non-hydrostatic
simulations are the most demanding of computation, that limit is the focus
of attention here.

\subsubsection{The elliptic problem}
\label{section:elliptic-problem}

Solving Equation (\ref{equation:elliptic-eq}) efficiently requires an
understanding of the form taken by the matrix $\A$; its structure is
sketched schematically in Figure~\ref{figure:matrix-a}.  The column vector
$p$ in Equation (\ref{equation:elliptic-eq}) comprises the pressure in
every zone in the ocean.  Its singly subscripted elements are $p_l$, where
first points in each vertical column are enumerated first, then points in
the horizontal directions, as follows:
\begin{displaymath}
l= k + N_z (j-1) + N_z N_y (i-1)
\end{displaymath}
where $i$ is an index increasing eastwards, $j$ is an index increasing
southwards and $k$ increases downwards.  There are $N = N_x \times N_y
\times N_z$ elements in total.

The elements of $\A$ comprise the discrete representation of the Laplacian
operator $\nabla^2$ (in the three space dimensions), suitably modified to
take account of boundary conditions on the upper, lower and lateral
boundaries; it is constant in time and depends only on the basin geometry.
$\A$ is symmetric and, although huge, has only 7 diagonals representing the
coupling in the 3 space dimensions.  The three central diagonals are
flanked, on both sides, by diagonals displaced a distance $N_z$ and $N_z
\times N_{y}$ from the central diagonal, as indicated in Figure
\ref{figure:matrix-a}.

\begin{fig}{figure:matrix-a}{Matrix $\A$, of size $N \times N$ where $N =
N_{x} \times N_{y} \times N_{z}$.}
\illustrationx{2.6in}{figures/matrix-a-new.ps}
\end{fig}

To write down the complete form of $\A$ consider the representation
of $\nabla^2 p$ for one particular zone (labelled $C$) in terms of the
six surrounding zones (labelled $U$, $L$, $N$, $E$, $S$ and $W$).  Suppose 
that the
dimensions of the zones are $\deltax$ by $\deltay$ by $\deltaz$, where
$\deltax$ is its length in the eastward direction, $\deltay$ is its length
in the southward direction, and $\deltaz$ is the length in the vertical
extent of the zone.  Then $\nabla^2 p$, to an accuracy which is second
order in the $\deltax$, can be written as follows (keeping $\deltax
=\deltay$ and $\deltaz$ constant for simplicity).\footnote{In applications,
$\deltax$, $\deltay$ and $\deltaz$ may vary and often do.}
\begin{small}
\begin{eqnarray}
\nabla^2 p & = & \frac{\pu + \pl - 2 \pc}{(\deltaz)^2} + 
                 \frac{\pn + \ps + \pe + \pw - 4 \pc}{(\deltax)^2} \nonumber \\
           & = & \frac{2}{(\deltaz)^2} 
                 [ \underbrace{\{ \half \pu - (1 + 2 \epsilon) \pc + \half \pl \}}
                              _{{\textstyle same \: vertical \: column}}
                   + \nonumber \\
           &   & \hspace{1cm}  \underbrace{\epsilon \{ \pn + \ps + \pe + \pw \}}
                              _{{\textstyle same \: horizontal \: plane}}
                 ]
\label{equation:seven-components}
\end{eqnarray}
\end{small}
where $\epsilon = (\frac{\deltaz}{\deltax})^2$ measures the aspect ratio of
the zones.

The detail of the structure underlying $\A$ is now readily understandable;
the three leading diagonals are the coefficients multiplying $p$'s in
the same {\it vertical column} of ocean ($\pu$, $\pc$ and $\pl$); the four
diagonals in the wings are the coefficients multiplying $p$'s in zones
surrounding the zone of interest in the same {\it horizontal plane}
($\pn$, $\ps$, $\pe$ and $\pw$).

\begin{fig}{figure:blocked-a}{Blocked representation of the Matrix $\A$.
$M = N_{x} \times N_{y}$, and each block is of size $N_{z} \times N_{z}$.
The $\D$ and $\E$ blocks are themselves tridiagonal and diagonal matrices,
respectively.}
\illustrationx{2.6in}{figures/blocked-a-new.ps}
\end{fig}

The inner three diagonals of $\A$ can thus be blocked and arranged as shown
in Figure~\ref{figure:blocked-a} where each block $\D$ is a tri-diagonal
matrix representing the communication between zones in the same vertical
column of ocean; $\D$ is a matrix of size $N_z \times N_z$ where there are
$N_z$ zones in each column of ocean (there are $N_x N_y$ such blocks, one
for each column of ocean).  The $\E$'s are $N_z \times N_z$ matrices made up
of $\epsilon$ along the diagonal and zero elsewhere:
\begin{displaymath}
\E = 
\left(
\begin{array}{ccccc}

\epsilon \\
         & \epsilon \\
         &          & \ddots   \\
         &          &          & \epsilon 

\end{array}
\right)
\end{displaymath}

For example, if the ocean were only four zones deep at a particular
horizontal position then by inspection of Equation
(\ref{equation:seven-components}) the $\D$ would have the form:
\begin{small}
\begin{displaymath}
\D =
\left(
\begin{array}{cccc}

\bterm & \half  & 0      & 0     \\
\half  & \dterm & \half  & 0     \\
0      & \half  & \dterm & \half \\
0      & 0      & \half  & \bterm

\end{array}
\right)
\end{displaymath}
\end{small}

Note that the elements in the top left and bottom right of the matrix have
been appropriately modified to take account of the boundary conditions at
the top and bottom of the ocean.  If at another horizontal position, the
ocean were only three zones deep, then $\D$ would have zeros in the last
row and column.

The entries of $\A$ must also be modified for zones that are laterally
adjacent to a rigid boundary, whether it be the coast, a submerged island
or an island that cuts the ocean's surface.  For example, suppose the north
face of zone C is the coast.  Then homogeneous Neumann conditions on $p$
must be imposed here (see \cite{Marshall97a,Marshall97b}) which, in
discrete form is:
\begin{displaymath}
\frac{p_N - p_C}{\deltax} = 0
\end{displaymath}
and the entries of $\E$ are appropriately modified; zeros will appear at
appropriate positions along the diagonal.  $\A$, however, will remain
symmetric.

One further property of $\A$ must be noted.  $\epsilon$ is typically very
much less than unity.  This is because the ocean is a thin film of fluid
filling the basin that is several thousand kilometers wide but only a few
kilometers deep -- the $\nabla^2 p$ operator is dominated by $\dzs$.  So
$\A$ is dominated by the elements of the blocks $\D$.  This is
exploited by the solution method.

The size of $\A$, and the distance of the outer diagonals from the main
diagonal make direct methods of solving the system difficult for large
geometries because the inverse of $\A$ is dense.  Operating with the
inverse would involve $N^2$ multiplications, and this is an unrealistic
task given that typically $N > 10^{5}$.  An iterative procedure is adopted
that exploits the sparseness of $\A$ and its diagonal dominance.  The
procedure involves repeated multiplication of the iterative solution by
$\A$ and by another sparse matrix $\K$, an approximate inverse of $\A$;
$\K$ is called the preconditioner.  The iterative algorithm employed to
solve the linear system is called the {\it preconditioned conjugate
gradient algorithm} and it is described in more detail in the next section.

\subsubsection{Preconditioned conjugate gradient algorithm}
\label{section:conjugate-gradient}

The preconditioned conjugate gradient (PCG) algorithm has been extensively
studied and is a popular algorithm for implementation on highly parallel
platforms because it is easily parallelized.  Barrett,
et.al.~\cite{templates} give a good overview of various iterative methods
for solving linear systems including conjugate gradient; Golub and Van Loan
\cite{GolL89} is a good general reference for matrix computations,
including conjugate gradient and Shewchuk \cite{Conjugate-Gradient} gives a
clear and light introduction to conjugate gradient.  Here, a brief
derivation of the preconditioned conjugate gradient method is provided.

Pre-multiply Equation (\ref{equation:elliptic-eq}) by a (carefully chosen)
preconditioning matrix $\K$, which is an approximate inverse of $\A$.
Equation (\ref{equation:elliptic-eq}) can then be written (dropping
superscripts) as
\begin{displaymath}
(\Idm - \Cm) p = \K f
\end{displaymath}
where $\Cm = \Idm - \K \A$.
If $\Cm$ is {\it close} to zero then the above suggests the following iterative
scheme, where $i$ is the iteration step: 
\begin{eqnarray*}
\pj & = & \Cm \pi + \K f \\
    & = & \pi + \bi
\end{eqnarray*}
where
\begin{displaymath}
b^{i} = \K \ri
\end{displaymath}
is called the {\em search direction} and
\begin{displaymath}
\ri = f - \A \pi
\end{displaymath}
is the {\em residual vector}. The $\ri$'s and $\bi$'s can be deduced
from the previous iteration using the relations: 
\begin{eqnarray*}
\rj & = & \ri - \A\bi \\
\bj & = & \K\rj
\end{eqnarray*}

The above iterative procedure, referred to as the Richardson iteration, can
be accelerated by choosing search directions in a more optimal way.  In
GCM, the {\em conjugate gradient method} is adopted, which selects search
directions as linear combinations of the previous search direction $b^{i}$
and current gradient $r^{i+1}$ modified by the preconditioner $\K$:

\begin{eqnarray*}
                                  b^{i+1}   & = & \xi^{i+1} + \beta b^{i} \\
\mbox{where} \mbox{\hspace{.3cm}} \xi^{i+1} & = & \K r^{i+1}
\end{eqnarray*}
Conjugate gradient also selects a parameter $\alpha$ to minimize the
magnitude of the residual vector as measured by $\ip{e}{\A e}$, where (if
there were no preconditioning) $e=p^{i}-p^{answer}$ is the error vector
along the direction $b^{i}$.  Choosing the optimal values of $\alpha$ and
$\beta$ results in the algorithm shown in Figure
\ref{figure:pcg}.  

\begin{fig}{figure:pcg}{The preconditioned conjugate gradient
algorithm, which solves the linear system $\A p=f$ in an iterative fashion.
The variable $\tau$ is the tolerance and $\pz$ is an initial estimate of
$\pfinal$.} 
\begin{tabbing}
a \= b \= \kill
\procname{PCG}($\A$, $f$) \\
\> $\rz = f - \A \pz;$ \\
\> $\bz = \K \rz$ \\
\> $\xi^{0} = \bz$ \\
\> \keyword{while} $\sqrt{\ip{\ri}{\ri}} > \tau$ \keyword{do} \\
\> \> $\gamma = (\ip{\ri}{\xi^{i}})$ \\
\> \> $\Ai = \gamma \; / \; (\ip{\bi}{\A \bi})$ \\
\> \> $\pj = \pi + \Ai \bi$ \\
\> \> $\rj = \ri - \Ai \A \bi$ \\
\> \> $\xi^{i+1} = \K \rj$ \\
\> \> $\Bi = (\ip{\rj}{\xi^{i+1}}) \; / \; \gamma$ \\
\> \> $\bj = \xi^{i+1} + \Bi \bi $ \\
\> \> $i   = i+1 $\\
\> \keyword{return}($\pfinal$)
\end{tabbing}
\end{fig}

We have designed a preconditioner $\K$ so that (i) it can be efficiently
stored, (ii) the number of operations one has to perform when multiplying
by it is as small as possible, (iii) it is a good approximation to
$\A^{-1}$, allowing the iterative procedure converge more rapidly and (iv)
it requires as little communication as possible in the data parallel
implementation.  After considerable experimentation, we have chosen a {\it
block-diagonal} preconditioner (a matrix $\K$ whose diagonal is made up of
the inverse of the tri-diagonal matrices $\D$ defined above):
\begin{displaymath}
\K =
\left(
\begin{array}{ccccc}

\D_{1,1}^{-1} \\
          & \D_{1,2}^{-1} \\
          &           & \D_{1,3}^{-1} \\
          &           &           & \ddots \\
          &           &           &        & \D_{N_x,N_y}^{-1}

\end{array}
\right)
\end{displaymath}
Evaluation of the inner products to compute $\Ai$ and $\Bi$ in the
above involves forming the vector $\xi$, where 
\begin{displaymath}
\xi = \K r
\end{displaymath}

Multiplication by the preconditioner $K$ is therefore $N_{x} \times N_{y}$
independent ($N_{z} \times N_{z}$) matrix multiplications, each one
corresponding to a column of water at $i,j$:
\begin{displaymath}
\xi_{ij} = \D_{i,j}^{-1} r_{i,j}
\end{displaymath}
Dispensing with the subscripts,
\begin{displaymath}
\xi = \D^{-1} r
\end{displaymath}
If the ocean model has many levels, then $\D^{-1}$ will also be dense;
consequently,
storing and multiplying by $\D^{-1}$ is also demanding of resources.
So, GCM exploits the fact that $\D$ is tri-diagonal, and use LU
decomposition to solve the preconditioning equations for $\xi$ in the form:
\begin{displaymath}
\D \xi = r
\end{displaymath}
We write $\D = \Lh \Uh$ where $\Lh$ is a lower-triangular matrix and $\Uh$
an upper triangular matrix.  Solving for $\xi$ is then equivalent to
solving the two sets of equations:
\begin{eqnarray*}
\Lh \xi' & = & r \\
\Uh \xi & = & \xi'
\end{eqnarray*}
first for $\xi'$ and then for $\xi$.  This is straightforward because of the
triangular structure of $\Lh$ and $\Uh$.  (See
\cite{tridiagonal-LU-decompo} for more details of the tri-diagonal LU
decomposition.)

Ignoring the operations for LU decomposition, these forward and backward
substitution steps involve $O(N_z)$ operations as compared to $O(N_z^2)$ if
we had used the inverse directly.  However, because the geometry (and
consequently $\A$) are fixed, LU decomposition of $\D$, the block diagonal
elements of $\A$, is performed only once.  Notice that the
multiplication by the preconditioner is actually implemented as $N_{x}
\times N_{y}$ independent tridiagonal linear equation solvers of size
$N_{z} \times N_{z}$.

The resulting preconditioner was found to be a good compromise: it
significantly reduces the number of iterations required to find a solution
to Equation (\ref{equation:elliptic-eq}), it is an acceptable approximation
to the inverse of $\A$, it is sparse and its application requires no
communication across the network in the data parallel implementation of the
algorithm.

\subsection{Computational characteristics}
\label{ss:computational-characteristics}

To effectively program GCM for a parallel platform, we first determine how
much computation must be performed in each step of the algorithm and then
consider what sources of parallelism are available.

\subsubsection{Distribution of overall computation time}

\begin{fig}{figure:computational-requirements}{Step 2
(\procname{diagnose\_pressure}) of the time stepping loop contains the
preconditioned conjugate gradient algorithm. In the situations we focus
on in this study, this step may iterate hundreds of times until convergence, 
making it the dominant step in the algorithm.}
\begin{center}
\begin{tabular}{|l|r|}
\hline
\multicolumn{2}{|c|}{\bf Floating Point Ops per Time Step per Zone}\\
\hline
Step 1: \procname{estimate\_Gv}         & 467 \\
Step 2: \procname{diagnose\_pressure}   & 32 $\times$ PCG iter. $+$ 32\\
Step 3: \procname{advance\_velocity}    & 20  \\
Step 4: \procname{advance\_temperature} & 122 \\
Step 5: \procname{advance\_salinity}    & 118 \\
Step 6: \procname{calculate\_density}   & 6   \\
\hline
{\bf Total} & 32 $\times$ PCG iter. $+$ 765 \\
\hline
\end{tabular}
\end{center}
\end{fig}

Figure~\ref{figure:computational-requirements} shows how many floating point 
operations are required in the GCM's main time-stepping loop.  There are also 
initialization and termination phases that are not considered in this table. 
These only occur once per model run, and will not be significant factor in long 
simulations.

From Figure~\ref{figure:computational-requirements} we can see that a fully
general implementation of GCM requires that the PCG algorithm be
implemented efficiently. For example if the number of PCG iterations
exceeds 300 per time step, the \procname{diagnose\_pressure} stage will
dominate the total count of floating point operations, containing more than
90\% of the operations.  This scenario is typical of extremely
high-resolution ($\sim$ 100m resolution laterally and vertically)
non-hydrostatic ocean simulations, the focus of this study.

Time is also spent in communications, non-floating point computation, and idling
due to load imbalances; however, for GCM, the style of computation in all of the 
steps is similar enough that floating point operation counts give a fair view 
of the distribution of time. To optimize GCM, we must concentrate on the time 
spent in non-essential floating point and other overhead operations, because 
given a particular model run, the number of essential floating point operations 
is fixed.

Although ocean modeling applications do perform substantial amounts of I/O
we do not analyze I/O costs here.
Rather, we focus on the comparison of
two different computational models. The difference 
in performance of I/O susbsystems is not predominantly a function of the 
computational model. Thus, to expose the impact of the computational models, 
we have removed I/O costs from both implementations.

\subsubsection{Distribution of computation in PCG}

\begin{fig}{figure:pcg-computational-requirements}{PCG is composed of four
simple components.  Each component performs a set number of floating point
operations per zone, and some components occur more than once in PCG.}
\begin{center}
\begin{tabular}{|l|r|r|r|}
\hline
\multicolumn{4}{|c|}{\bf Floating Point Ops per PCG iteration per Zone}\\
\hline
{\bf PCG component}  & {\bf occurrences} & {\bf FLOPs} & {\bf Total} \\
\hline
3D inner product     &                 4 &           2 &           8 \\
3D daxpy             &                 3 &           2 &           6 \\
7-point stencil      &                 1 &          13 &          13 \\
preconditioner       &                 1 &           5 &           5 \\
\hline
{\bf PCG Total}      &                   &             &          32 \\
\hline
\end{tabular}
\end{center}
\end{fig}

The PCG algorithm itself can be further decomposed into four simple components: 
3D inner product, 3D daxpy, a 7-point stencil and the preconditioner.  These 
components can be seen in the algorithmic outline in Figure \ref{figure:pcg}; 
the actual number of floating point operations contributed by each 
component is shown in Figure \ref{figure:pcg-computational-requirements}.

Each of the components operates on three-dimensional state arrays.  The
7-point stencil component implements $\A b$ (the multiplication of the
Laplacian operator) and the preconditioner component performs the
multiplication of the preconditioner $\K$, which is actually implemented as
a tridiagonal linear equation solvers (see Section
\ref{section:conjugate-gradient}).

To efficiently implement the PCG algorithm, each of the components must
also be implemented efficiently.  In Sections
\ref{section:programming-fortran} and \ref{section:programming-id}, the
implementations of these components in both the data parallel and
multithreaded models are described.

\subsubsection{Available parallelism}

Most of the potential parallelism in GCM comes from the {\it data
parallelism} in updating each of the large ocean state arrays.  This
parallelism is typically sufficient to keep a large parallel computer busy
and, because this computation is so regular, it lends itself well to vector
computation.

There is also {\it procedure-level parallelism} available: for example,
much of Step 4 can be done in parallel with Steps 1, 2 and 3, and within
individual steps some blocks can be executed in parallel.  In Step 2, for
instance, at least three separate state arrays can be computed in parallel.
This source of parallelism can be exploited by the multithreaded
implementation, but not by the data parallel implementation.

{\it Producer-consumer parallelism} is also plentiful in GCM, because the
computation of new state arrays always depends upon old state arrays.
Given appropriate linguistic and architectural support, the computation of
parts of new state arrays can begin as soon as the parts of the old state
arrays they depend upon are computed.
