\section{Performance profiling}
\label{s:experiments}

To obtain high performance for the MIT model on a new computer platform, 
we
take an incremental approach, beginning with performance profiling of
elementary kernels, and moving towards profiling of the full application.
At each stage, we gain a better quantitative understanding of the
architecture, programming environment, and algorithmic match to the
machine.  The successive tests give us a step-by-step guide to the design
of the full application.

\subsection{{\bf DAXPY} test}

The {\bf DAXPY} test measures the performance of a system for vector-style
operations.  {\bf DAXPY} is a primitive of the BLAS1 \cite{BLAS-DESC}
linear algebra library which performs the operation
$\underline{y}=a\underline{x}+\underline{y}$, where $\underline{y}$ and
$\underline{x}$ are 64-bit floating point arrays.

\begin{figure}[t]
\illustrationx{4in}{figures/daxpy-comparison.ps}
\caption{Performance  of {\bf DAXPY} BLAS1 primitive for  different 
lengths
of vector on one, two and four processors of a DEC 8400, and on a 32 node
CM-5 computer.}
\label{figure:DAXPY}
\end{figure}

Figure \ref{figure:DAXPY} shows the floating point performance as a
function of the length of the input arrays for the 32 node CM-5, and for 
different
processor counts of the DEC 8400.  Examining this graph, we clearly see 
the
impact of the cache size on performance.  For the 8400 system, performance
tails off sharply when cache size is exceeded, which is indicated at the
``cache shoulders''.  As the number of processors increases, the total
amount of cache memory also increases, as seen by the cache shoulder 
moving
towards the right as more processors are used in the 8400.

For the CM-5, the performance increases with the size of the problem,
because the CM-5 can more fully use its vector units on larger sized
problems.  In comparing the absolute performance of the CM-5 to the 8400,
we should try to determine what region of interest is delineated by the
curves.  Once the size of the problem is too large to fit into the caches
of the 8400, the performance is uncompetitive with the CM-5, or like
machines.

In some ways, {\bf DAXPY} is the worst-case situation for the SMP, because
no data is reused, and therefore, the caches are poorly exploited.  We can
attempt to alleviate the effects of computations like {\bf DAXPY}, by
reusing data that is in the cache.  However, traditional vector
supercomputers favored applications which were structured as compositions
of DAXPY-like operations, so programming and compiling for
microprocessor-based SMP's requires a re-thinking of the mapping of an
algorithm to the computer.

\subsection{{\bf 3D DAXPY} test}

{\bf 3D DAXPY} is a three-dimensional version of {\bf DAXPY}.  Although 
the
{\bf 3D DAXPY} computation seems extremely straightforward, there are
actually some subtleties about programming even such a simple piece of
code.  Figure \ref{fig:axy-axis-layout} shows some of the results obtained
when we ordered the axes of the arrays differently, and when we used
different HPF data layout directives for the same code using a fixed
problem size and a fixed number of processors.  We ran the code with 
arrays
of size $200 \times 100 \times 10$, which is representative of runs we do
daily on the CM-5.

\begin{fig}{fig:axy-axis-layout}{The performance of {\bf 3D DAXPY} on the 
DEC
8400 under the DEC HPF compiler is  highly dependent on the axis order
declaration and data layout directives.  The  3D arrays are of extents
$X_{len}=200,Y_{len}=100,Z_{len}=10$, and  runs were  on a 4 processor
8400.} 
\begin{center}
\begin{tabular}{|c|c|c|r|}
\hline
\multicolumn{4}{|c|}{{\bf 3D DAXPY }performance on DEC 8400}\\
\hline
{\bf Version} & {\bf Axis Order} & {\bf Data Layout} & {\bf Mflops} \\ 
\hline
1 & {\tt ($z$,$x$,$y$)} & {\tt (BLOCK,*,*)} & 60 \\ 
\hline
2 & {\tt ($z$,$x$,$y$)} & {\tt (*,BLOCK,BLOCK)} & 210 \\ 
\hline
3 & {\tt ($x$,$y$,$z$)} & {\tt (BLOCK,BLOCK,*)} & 260 \\ 
\hline
4 & {\tt ($y$,$x$,$z$)} & {\tt (BLOCK,BLOCK,*)} & 230 \\ 
\hline
5 & {\tt ($x$,$y$,$z$)} & {\tt (*,BLOCK,*)} & 270 \\ 
\hline
\end{tabular}
\end{center}
\end{fig}

There is almost a factor of 5 difference between the performance of 
Version
1 and Version 5.  Version 2 actually uses the same axis order and layout 
as
the CM-5 version of the MIT model, but it is not optimal on the 8400. It 
is
clear that Version 5 gives the maximum possible performance, because its
performance is about the same as the one-dimensional {\bf DAXPY} test.

\subsection{{\bf Forward\_Step} test}
\label{forward_step_test}

The {\bf Forward\_Step} part of the code consists of updates to a three
dimensional grid of double-precision floating point numbers representing
the state of the fluid.  In our study, in keeping with our goal of a highly
intuitive and readable code, the grid is mapped directly to
three-dimensional arrays.  Figure \ref{fig:forward-scatter} shows the
performance of a representative fragment from {\bf Forward\_Step}, which
calculates the $G$ terms in Equation \ref{eq-Vt}.  Problem sizes were
chosen to be wide and thin, which is representative of geometries we use
for our research.  Times were gathered by extracting the section of code
from the model, and using the optimal axis ordering and data layout as
determined by the {\bf 3D DAXPY} test.

\begin{fig}{fig:forward-scatter}{Performance of the {\bf Forward\_Step}
test as a function of problem size.  The shoulder for the DEC 8400 has
moved to about $10^{5}$ elements, because many more arrays are used for
{\bf Forward\_Step} than for {\bf DAXPY}.}
\illustrationx{4in}{figures/forward-scatter.ps}
\end{fig}

The results in Figure \ref{fig:forward-scatter} represent the performance
after modest optimizations were made to the DEC 8400 code relative to the
CM-5 code.  The CM-5 version actually uses about 40 arrays, some of which
are used to minimize communication.  With the DEC 8400 shared-memory
architecture it is more important to reduce the overall memory usage, so
some of these buffer arrays are removed.  Memory requirements were further
reduced by storing time invariant terms as 32-bit rather than 64-bit
values.

Note that the performance of the DEC 8400 does not drop off as much as for
{\bf DAXPY} with this more representative code extract, but that the
drop-off occurs at an earlier point.  The improvement in asymptotic
performance reflects the higher ratio of floating point operations to
memory references of the {\bf Forward\_Step} test (which has a ratio of
$\approx 26:8$) relative to the {\bf DAXPY} tests (where the ratio is
$1:1$).  This higher ratio translates into increased opportunities for
reusing data loaded from main memory into cache.  The degree to which this
is exploited is highly dependent upon optimizations performed by the
compiler and the way in which the code is written by the user.

The earlier drop-off point relative to {\bf DAXPY} is the result of using
more arrays -- whereas {\bf DAXPY} uses 2 arrays, this code segment uses
the equivalent of 24 64-bit arrays, and accordingly, the drop-off point is
about a factor of 12 earlier than for {\bf DAXPY}.  The full MIT model
requires about the same number of arrays (to hold six state variables,
their time derivatives at two time levels and geometric terms), so we
expect that the rest of the {\bf Forward\_Step} code should behave as this
code fragment does.  To put this problem size in perspective, it is about
equivalent to a global simulation with a four degree horizontal grid and
twenty vertical layers.

\subsection{{\bf Inversion\_Step} test}

The other important stage in the model algorithm is the {\bf
Inversion\_Step} which solves an elliptic problem to ensure non-divergent
flow.  In discrete form this inversion can be written:
\begin{equation}
\label{eq-Axb}
\underline{\underline{A}}\underline{p}=\underline{f}
\end{equation}
in which, for hydrostatic simulations, $\underline{\underline{A}}$ 
is the horizontal Laplacian operator,
a matrix with five diagonals, and $\underline{p}$ and $\underline{f}$
represent the surface pressure and the right-hand side term of Equation
\ref{eq-Ps} respectively.

In the CM-5 implementation of the model, Equation \ref{eq-Axb} is solved
using an iterative pre-conditioned conjugate gradient (PCG) scheme.  PCG is
well suited to data parallel implementations and the same algorithm was
implemented in HPF on the DEC 8400.

\begin{fig}{fig:cg2d-scatter}{The preconditioned conjugate gradient
algorithm used to solve for the pressure is competitive to the CM-5 over
the region of interest ($10^{4}$ to $10^{5}$).}
\illustrationx{4in}{figures/cg2d-scatter.ps}
\end{fig}

Figure \ref{fig:cg2d-scatter} shows the performance of the PCG algorithm on
the CM-5 and DEC 8400.  Note that because we are solving a two-dimensional
problem instead of a three-dimensional problem (required in the
non-hydrostatic limit formulation of the model), the problem sizes we
consider are much smaller, and thus the performance of the DEC 8400 does
not show the characteristic cache shoulder of the previous tests.

In addition, we perform further optimizations to minimize memory use.  In
the CM-5 implementation of PCG, to reduce communication we duplicate data
structures that hold the diagonals of $\underline{\underline{A}}$ shifted
by one point in $x$ and $y$.  On the DEC SMP system, this duplication is
not necessary, and we exploit the symmetry of $\underline{\underline{A}}$
to eliminate another array by recalculating its value at run-time.
Finally, we also store $\underline{\underline{A}}$ in 32-bit precision
rather than 64-bit.  In total, we were able to reduce the memory required
to represent $\underline{\underline{A}}$ by a factor of 5.

The PCG iterative solution can require anywhere between 10 and 200
iterations to converge -- the exact number will depend on the physics of
the particular simulation.  For a domain with a surface of size $NX \times
NY$, if we assume that convergence requires $nits$ iterations, the number
of floating point operations executed is:
\[nits \times 32 \times NX \times NY\]
where 32 is the number of floating point operations performed per grid
point per iteration in our PCG implementation.

An alternative approach to an iterative PCG solver is a direct solver.
Because $\underline{\underline{A}}$ is symmetric and sparse, the LU
decomposition can be stored within the profile of the
$\underline{\underline{A}}$ matrix.  For any $\underline{\underline{A}}$
representing a closed physical domain, this requires about $NX \times NX
\times NY$ elements of storage.  Once the LU factorization has been
accomplished, we can find $\underline{p}$ using forward-backward
substitution.  Solving for $\underline{p}$ using this {\it direct} method
requires the following number of floating point operations:
\[3 \times NX \times NX \times NY\]
where 3 floating point operations are performed per word of the LU
factorization for the forward-backward substitution.

\begin{fig}{fig:skyline-vs-pcg}{The direct solution to $Ax=b$ is divided
into a factorization of $A$ and the solve for $x$.  Factorization only
occurs once during the run of a program, and its cost is amortized out
on long runs.  The Skyline direct solver is faster  in wall clock time
than  the PCG iterative  solve, even though it   has a lower Megaflops
rating.} 
\begin{center}
\begin{tabular}{|c|r|r|c|r|c|}
\hline
\multicolumn{6}{|c|}{\bf Skyline (Direct) vs. PCG (Iterative) solvers on DEC 8400}\\
\hline
{\bf Ocean} & \multicolumn{3}{|c|}{{\bf Skyline} (1 proc)} & 
\multicolumn{2}{|c|}{{\bf PCG} (4 proc)} \\
\cline{2-6}
{\bf   Size} & {\bf  Factor}  & {\bf Solve}   & {\bf Megaflops} & {\bf
Solve} & {\bf Megaflops} \\ 
\hline
$100 \times 100$ & 1.1 s. & 0.1 s. & 9.1 & 0.1 s.  & 350 \\ 
\hline
$200 \times 100$ & 2.4 s. & 0.2 s. & 8.7 & 2.7 s.  & 380 \\ 
\hline
$200 \times 200$ & 23.0 s. & 0.9 s. & 8.8 & 1.6 s.  & 467 \\ 
\hline
$400 \times 200$ & 49.7 s. & 1.9 s. & 8.4 & 4.6 s.  & 480 \\ 
\hline
$400 \times 400$ & 238.0 s. & 10.0 s. & 6.4 & 25.8 s.  & 230 \\ 
\hline
$800 \times 400$ & 633.3 s. & 20.9 s. & 6.1 & 83.5 s.  & 210 \\ 
\hline
\end{tabular}
\end{center}
\end{fig}

Because of the way the problem is posed, we can choose to appropriately
switch $NX$ and $NY$ to arrive at a problem with a smaller LU
decomposition.  We should note, however, the lateral aspect ratio of
typical geometries is not usually greater than $2:1$, so that the direct
method has an algorithmic complexity $O(n^{3})$, whereas the iterative
methods generally require about $O(n^{2})$.

Despite the higher order complexity of the direct method, the high constant
factors in the iterative method make the direct method competitive for a
certain range of problem sizes.  Suppose that a $100 \times 100$ sized
problem required 100 PCG iterations to converge -- the PCG algorithm would
require $100 \times 32 \times 100 \times 100 = 32,000,000$ floating point
operations, whereas the direct method would only require $3 \times 100
\times 100 \times 100 = 3,000,000$ floating point operations.  However,
even though the direct method may execute fewer floating point operations,
each floating point operation will be much slower than the direct method
because the direct method must stream through the entire LU factorization,
performing only 3 operations per element, leading to a fairly bad ratio of
floating point operations to memory operations.

To test the effectiveness of a direct approach versus the iterative
approach, we solved the same 2D problem using a sparse direct solver
included with the DEC DXML scientific library \cite{DEC-DXML}.  The solver
uses the skyline algorithm, which assumes matrices have elements which are
primarily close to the main diagonal.  This particular algorithm does not
assume any other structure (other than symmetry) within the
$\underline{\underline{A}}$ matrix, and currently only runs on 1 processor
of the 8400 system -- a more specialized or parallel direct solver may be
faster than this very general-purpose, single-processor algorithm.

Nevertheless, the results in Figure \ref{fig:skyline-vs-pcg} show that a
direct solver could be a competitive solution method for certain problems.
For example, Hill and Marshall\cite{PCFD95} make a non-hydrostatic model
competitive with a hydrostatic model, even in the presence of complex
geometry, by separating out, and solving to high-accuracy, the
two-dimensional surface pressure problem we examine here.  Because the
direct solver returns a solution to Equation \ref{eq-Axb} which is accurate
to machine precision, the solver would be ideal for accelerating a
non-hydrostatic model.

\subsection{Full application performance}
 
We evaluate the full model's DEC 8400 implementation by comparing its
performance with the CM-5 implementation for a specific physical
simulation.  The simulation we chose was an experiment modeling the
collapse and break up, through baroclinic instabilities, of an anomalously
dense column of fluid in a rotating frame of reference.\cite{Jones96} The
chosen test has a computational grid of size $100 \times 100 \times 26$,
which is large enough to exceed the cache size of a four processor DEC
8400, but not so large that there are no opportunities for optimization.

\begin{figure}[t]
\illustrationx{3.5in}{figures/fullmodel.ps}
\caption{Performance of the full model on a 32 node CM-5 and a
DEC 8400 5/300 server.}
\label{figure:FULL_TEST}
\end{figure}

Figure \ref{figure:FULL_TEST} shows the relative performance of the CM-5
and DEC 8400 for this simulation.  The DEC 8400 four processor system is
about one half as fast as the CM-5.  From our previous tests, we should be
able to predict the performance of the DEC 8400 given the problem size.
For the {\bf Forward\_Step}, referring to Figure \ref{fig:forward-scatter},
the problem size is $2.6 \times 10^{5}$, which shows that the DEC 8400
should be about 2.5 times slower than the CM-5.  For the {\bf
Inversion\_Step}, referring to Figure \ref{fig:cg2d-scatter}, the problem
size is $10^{4}$, so the DEC 8400 should be about 2.5 times faster than the
CM-5.

For the problem that we simulated in Figure \ref{figure:FULL_TEST}, the
time spent in one time step on the CM-5 is about equally divided between
the {\bf Forward\_Step} and the {\bf Inversion\_Step}.  This ratio will
depend on the physics of the simulation, in particular, the amount of work
required by the {\bf Inversion\_Step} varies greatly from one problem to
another.  Notice that our performance prediction was fairly accurate: the
time spent in the {\bf Inversion\_Step} for the DEC 8400 is indeed about
2.5 times less, and the time spent in the {\bf Forward\_Step} is roughly
2.5 times greater than the CM-5 time.

Although the particular integration considered here had a 50/50 balance
between forward and inversion steps on the CM-5, other simulations have a
different split. For problems that spend a large fraction of time in the
inverter, the 4-processor DEC 8400 should out-perform the CM-5.

We have also compared single processor performance of our HPF test codes
expressed in Fortran 77, using the DEC compilers.  In these checks, the HPF
test code performance was always within 10\% of its Fortran 77 equivalent.
This, along with the linear scaling shown in Figure \ref{figure:FULL_TEST},
leads us to conclude that an explicit message-passing version of the MIT
model (for instance, coded using PVM or MPI) would have similar performance
to our HPF implementation.  Furthermore, unlike a message-passing code, we
are able to use a global address space throughout the model, making the
code much easier to understand and modify.
