Suppose that we found the LU factorization \cite{tridiagonal-LU-decompo} of
the $A$ matrix as shown in Figure \ref{fig:a-matrix}.  Because $A$ is
symmetric and sparse, it turns out that the LU decomposition can be stored
within the profile of the $A$ matrix -- that is, the spaces between (and
including) the main diagonal and the outermost diagonal.  For the $A$
matrix, this requires about $NX \times NX \times NY$ elements of storage.
However, once we find the LU factorization, we can find the solution to $x$
given the $b$ vector using forward-backward substitution.  This algorithm
requires 3 floating point operations per element of the LU factorization.
Solving for $x$ using this {\it direct}, as opposed to {\it iterative}
method requires the following number of floating point operations:
\[3 \times NX \times NX \times NY\]

Because of the way the problem is posed, we can choose to switch $NX$ and
$NY$ so that we have a problem with a smaller LU decomposition.  We should
note, however, the aspect ratio of ocean 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})$.  Due to the computational constraints imposed by the rest of the
program, we are usually considering problem sizes within a certain range,
and the constant factors due to the number of iterations to convergence
play an important factor.

The direct method will require much more memory than the PCG solution.  To
store the LU factorization, we need $NX \times NX \times NY$
double-precision floating point numbers.  For a $200 \times 100$ sized
ocean, we require about 16 Megabytes to store the LU factorization.
Furthermore, the forward-backward substitution is similar to the daxpy code
in that once an array element is used, it will not be used again, so that
the computation must stream through the entire factorization, paying the
penalty for accessing main memory.

Note that the direct approach is only feasible with a 2D problem, because
the 3D problem would require an $A$ matrix whose factorization is many
orders of magnitude larger.  The transformation of the 3D problem to a 2D
problem enables us to consider the direct solution of the linear system.

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 $A$ matrix, and
currently only runs on 1 processor of the 8400 system.

Figure \ref{fig:skyline-vs-pcg} shows the results of a simple experiment
comparing the direct and iterative methods.  Although the factorization for
the direct method can be expensive, it is only performed once at the
beginning of the program.  The ``solve'' portion of the direct method is
faster for the sizes we have tested, even though the Megaflops ratings are
much smaller, simply because much fewer operations are needed for the
direct method.

We should note that there are direct methods superior to Skyline,
especially for matrices with more structure.  Also, there are parallel
direct solvers, although we have used the sequential one provided in the
DEC DXML library, and we have compared its performance to the
four-processor performance of the iterative solver.

These results are preliminary, but serve to show that more Megaflops do not
necessarily indicate faster execution.  The direct method performs well
because it requires fewer operations, within the problem size range we have
tested, and because the DEC 8400 provided us with enough real memory to
hold the factorization.
