\section{Programming multithreaded GCM}
\label{section:programming-id}

The data parallel model of computation is ideally suited for an application
like GCM, that uses large, regular data structures, and where most of the
work is in element-wise operations over these data structures.  For other,
less regular applications, the data parallel model is too restrictive in
the way that parallelism can be expressed.  Multithreading is a more
general model of parallel computation that allows different processors to
work on different tasks at the same time in a less synchronous fashion.
Consequently, multithreading has some overheads that result from more
dynamic scheduling and synchronization.  In the case of GCM and other
highly structured codes, the multithreading overheads may be unnecessary 
because the applications may not require the expressiveness multithreading
offers.  One goal of this study was to quantify and characterize
multithreading overhead and contrast it with some of the overheads of data
parallel execution.

Figure \ref{figure:fully-parallel} gives a high-level view of the
multithreaded execution model.  Each procedure may be associated with an
activation frame, which is local to a processor; a procedure may fork off
child procedures that may execute in parallel with the original procedure,
either on the same processor or a different one.  Similarly, the iterations
of a loop may also execute in parallel.  These procedures and loops may
access data in a global heap of shared objects.  Multithreading can be
expressed in various languages and can be executed on different kinds of
machines.  The semantics of threads and styles of communication and
synchronization between threads vary significantly, as can the typical
length of threads, depending upon the language and architecture.

\begin{fig}{figure:fully-parallel}{A high-level view of multithreaded
execution; activation frames are associated with procedure invocations
and parallel loop iterations while data structures reside in a shared
global heap.}
\illustrationx{3.2in}{figures/fully-parallel.ps}
\end{fig}

The programming language used to exploit multithreading may be explicitly
parallel such as Cilk \cite{Cilk95} or Mul-T \cite{Kranz89} (allowing the
programmer to annotate procedures or loops that should be executed
in parallel), or implicitly parallel such as Id \cite{id} (allowing the
compiler to decide how to parallelize the code).  Furthermore, the memory
models for the languages can differ considerably. For example, most data
structures in Id are based on I-structures, which are implicitly
synchronizing; a fetch of an I-structure element will not return until
the store to that element is executed.

The degree of hardware support for multithreading may also vary, as Cilk
has been implemented on the CM-5, networks of workstations, and SMPs, none
of which have special support for multithreading.  Mul-T has been
implemented on Alewife, which has support for scheduling and synchronization
primarily at the basic block level, \cite{Alewife} and Id has been
implemented on Monsoon \cite{monsoon}, which has support for multithreading
down to the instruction level.  Monsoon also has direct support for
I-structures. Other architectures that support instruction-level
multithreading and I-structure-like memory include HEP \cite{HEP}, Sigma-1
\cite{Hiraki87b}, EM-4 \cite{EM4a}, and Tera \cite{Tera}.

We have implemented GCM in Id on Monsoon, which is one extreme of
multithreaded architectures.  In doing so, we try to exploit parallelism at
the instruction-level, procedure-level, and loop-level, incurring various
execution overheads despite extensive hardware support for multithreading
in Monsoon.  We begin our discussion of the multithreaded implementation of
GCM with a description of the Monsoon hardware and its execution model
that is based on dataflow graphs.  We then discuss how GCM was coded in
Id.  Next, we describe some practical problems in running GCM for realistic
data sets on Monsoon.  Finally, we describe optimizations that are
performed by the compiler and programmer to obtain good performance on
Monsoon and give some performance figures for the multithreaded
implementation.

\subsection{The Monsoon dataflow machine}

Figure \ref{figure:monsoon} shows a high-level view of the Monsoon
multiprocessor.  Monsoon \cite{monsoon} consists of eight 64-bit processing
elements (PE's) and eight I-structure memory modules (or IS's) connected
via a multi-stage packet switched network.  Each PE is an eight-stage
pipelined processor that handles dataflow {\it tokens}.  On each processor
cycle, a token enters the pipeline and its result is available eight cycles
later.  Each token carries with it (1) a value, (2) a pointer to an
instruction to execute, and (3) a pointer to a context, or activation
frame, in which to execute the instruction.  The execution of an
instruction may cause the creation of up to two more tokens, which may be
circulated locally or sent out on the network to other PE's.  Each PE can
process 10 million tokens per second and has 256K 64-bit words of local
memory where activation frames are allocated.

Globally shared heap objects reside in the IS's, each of which contains 4M
64-bit words of memory.  Both frame and heap memory have presence tags
associated with each word to support fine-grain synchronization.  Each IS
memory access goes over the network.  Monsoon's network has a bandwidth of
100 MB per second, and the network interface can accept or deliver a
message every cycle.  Messages are sent from PE's to IS's for heap
references, or from PE to PE to communicate data between procedures or
loops.  All heap locations are equidistant to all the PE's.  Heap objects
are interleaved by hardware across the IS's and in general, data mapping
has no effect on performance.

\begin{fig}{figure:monsoon}{Monsoon consists of processing elements (PE's)
connected to I-structure memory modules (IS's) through a multi-stage
network.  Each PE is an eight-stage pipelined dataflow processor.}
\illustrationx{3.2in}{figures/monsoon.eps}
\end{fig}

Every Monsoon instruction executes in one or two cycles.  A binary
instruction synchronizes on the arrival of its two inputs: the data from
the first input token is stored in the activation frame and a {\it bubble}
is passed in the pipeline, while the arrival of the second token causes the
instruction to execute, producing more tokens.  Unary instructions execute
in one cycle without a bubble, as do binary instructions with one constant
input.  Global memory accesses are handled in a split-phase manner; a
message is sent out to the network to the IS memory, and a result or
acknowledgement comes back from the network at an arbitrary time later.
During this time, Monsoon can process other tokens, using parallelism to
mask the latency of the memory operation.  Keeping an eight processor
Monsoon busy requires at least 64-fold parallelism because each stage of
the pipeline executes a different thread, and some additional parallelism
to mask the latency of heap memory references.

To measure performance, each Monsoon processor has 16 cycle counters that
can keep exact instruction counts (these counters were used for all
measurements in this paper).  We also have a cycle-level simulator, called
MINT (Monsoon Interpreter), which is capable of running small programs.
MINT has played a crucial role in micro-benchmarking and enhancing our
understanding of various factors that affect performance.

\subsubsection{Dataflow graphs for parallel execution}
\label{sss:dataflow-graphs}

\begin{fig}{figure:dataflow}{Monsoon dataflow graph for the body 
of the innermost loop of a three-dimensional inner product code.}
\illustrationx{2.7in}{figures/dfg-3d.ps}
\end{fig}

Unlike conventional RISC processors, the machine language of Monsoon is
{\it dataflow graphs}.  Dataflow graphs specify only a partial order on
instructions and thus implicitly represent instruction-level parallelism;
in contrast, conventional superscalar RISC processors dynamically detect
opportunities for instruction-level parallelism within a linear instruction
stream.  To illustrate several features of the Monsoon execution model,
consider the innermost loop of the 3D inner product:

\[\sum_{\mbox{\tt k} = \mbox{\tt kl}}^{\mbox{\tt ku}} \mbox{\tt a[k,j,i]*b[k,j,i]}\]

Assuming the arrays are stored in row major order, Figure
\ref{figure:dataflow} shows the Monsoon dataflow graph for the body of the
loop, including the loop predicate.

The nodes of a dataflow graph represent instructions, and dataflow tokens
flow along the arcs to pass values generated by parent instructions to
their children.  Because Monsoon instructions can send a value to at most two
destination instructions, the {\it fork} operator is required when a value
values is needed by more than two destination instructions.  Because of 
instruction encoding limitations, some instructions,
can send a value
to only one destination instruction.  The arcs leading from the two {\it
i-fetch} instructions to the {\it floattimes} instruction are drawn as
dotted arcs to emphasize that I-structure references are split-phase.

Literals and loop constants are represented in light gray in Figure
\ref{figure:dataflow}.  Literals are stored into local memory at load-time
by the loader, while loop constants must be stored into local frame memory
at run-time using the {\it constant-store} instruction.  Both types of
constants can be accessed at run-time without incurring a bubble in the
pipeline.

\begin{fig}{figure:dataflow-loop}{Sequential loop implementation of the
inner loop of 3D inner product.  The sequential loop executes on one PE,
and only exploits instruction-level parallelism.}
\illustrationx{3.2in}{figures/dataflow-loop-ip3d.eps}
\end{fig}

Control in Monsoon is implemented by using the {\it switch} operator that
steers an input value to one of two locations, depending upon the value of
an input predicate.  Figure \ref{figure:dataflow-loop} shows the
``sequential'' loop schema for the loop body shown in Figure
\ref{figure:dataflow}.  Before the first iteration of the loop begins, 
loop constants are stored into the activation frame.  The initial values
for the input variables are fed to the top of the loop body, and then the
loop body in Figure \ref{figure:dataflow} is executed.  The resulting values
are either steered back to the top of the loop, or else output as the final
values.  The values are steered using {\it switch} instructions, with one
{\it switch} for each induction variable.  If there are more than two
induction variables, forks are needed to distribute the predicate value to
the switch's.  Because the switch instruction can steer a value to only one
destination instruction, additional forks are often required after a
switch.

Sequential loops only execute on one processor, and generally use the
activation frame of the surrounding procedure.  In the next section, we
show how inter-processor parallelism is exploited in Monsoon for both
procedure-level parallelism and loop-level parallelism.

\subsubsection{Run-time system and inter-processor parallelism}
\label{sss:kbounded-loops}

In Monsoon, dataflow-style instruction-level parallelism is exploited only
within a PE, to keep the eight pipeline stages busy.  Inter-processor
parallelism is exploited by allocating activation frames on remote PE's.
Activation frames are allocated when a procedure is called or a parallel
loop is initiated, and the frame may be allocated on an arbitrary Monsoon
PE, although frames do not migrate once they are allocated.  Thus, the
mapping of work to processors is governed by the frame allocation policy,
which is a part of the Monsoon run-time system (RTS).

In the Monsoon RTS \cite{frame-manager}, frame memory on each processor is
divided equally among all processors, so that each processor manages part
of the frame space of every other processor.  Frame allocation requests are
then handled locally in a round-robin fashion across all of the processors
to provide load balancing.  Since the RTS is invoked on every procedure and
parallel loop invocation, care has been taken to minimize RTS costs, to the
point where much of the RTS is written in assembly and supported by special
``microcoded'' machine instructions.  Despite all this care, allocating and
deallocating a frame on Monsoon typically takes about 40 to 50 cycles.
This shows that frame management in this parallel environment is still more
expensive than in a sequential environment, where a simple pointer-bump can
allocate or deallocate a stack frame.

A procedure call begins with allocating an activation frame.  Once the
frame is returned by the RTS, the arguments to the procedure and the return
continuation are sent as messages to the frame, which is usually on a
remote processor.  Although the composition, sending and receiving of these
messages are well supported in Monsoon hardware, argument-passing is still
relatively more expensive than register-style argument passing in
sequential processors.  In general, the programmer should define most small
leaf procedures as inlinable to avoid procedure call overhead.

Loops can trivially exploit inter-processor parallelism if they are
implemented as recursive procedure calls.  However, such an implementation
usually uses too many frames and consequently has a very high overhead both
in terms of storage and RTS calls.  Id on Monsoon exploits an alternative
compiling scheme known as $k$-bounded loops \cite{Culler90}.  The
$k$-bounded loop schema employs a ring of $k$ activation frames that are
allocated and linked at the beginning of the loop execution and deallocated
at the end of the loop execution.  Loop constants must be stored into each
frame in a $k$-bounded loop, and induction variables are passed from
iteration to iteration like arguments are passed for procedure calls.

\begin{fig}{figure:loop-cycle-counts}
{Comparison of cycles per iteration for unrolled, sequential and
$k$-bounded loop implementations.  The number of
bubbles incurred by each operation per iteration are set off in parentheses.}
\begin{center}
\small
\begin{tabular}{|l|rrr|}
\hline
\multicolumn{4}{|c|}{\bf Cycles/iteration of the 3D inner product inner loop} \\
\hline
           & completely &            &           \\
{\bf Opcode} & unrolled & sequential & k-bounded \\
\hline
fork         &        2 &         3 &        6  \\
intplus      &        1 &         1 &        1  \\
i-fetch      &        2 &         2 &        3  \\
floattimes   &    2 (1) &     2 (1) &    2 (1)  \\
floatplus    &   1 (.5) &    1 (.5) &    2 (1)  \\
\hline
switch       &          &     4 (2) &        1  \\
intle        &          &         1 &        1  \\
\hline
message      &          &           &    4 (2)  \\
i-take       &          &           &        1  \\
i-put        &          &           &    2 (1)  \\
sync         &          &           &    3 (1)  \\
\hline
{\bf Total}  &  8 (1.5) &  14 (3.5) &   26 (6)  \\
\hline
\end{tabular}
\end{center}
\end{fig}

Exploitation of inter-processor parallelism can be expensive even on
Monsoon.  In Figure \ref{figure:loop-cycle-counts}, the cycle counts for
$k$-bounded, sequential, and completely unrolled implementations of the
innermost loop of the 3D inner product are shown.  Recall that the
sequential loop schema only exploits intra-processor parallelism and a loop
can be completely unrolled only if the loop bounds are known at compile
time.  It takes three times as many cycles per iteration for the
$k$-bounded loop as the completely unrolled loop!

We explain these cycle counts further to satisfy the reader's curiosity.
The completely unrolled implementation reflects the instructions in Figure
\ref{figure:dataflow} excluding the instructions needed to evaluate the
loop predicate.  The count for the sequential implementation can be derived
by adding the predicate and the two switches shown in Figure
\ref{figure:dataflow-loop}.  The $k$-bounded implementation uses two
message instructions to send the induction variables to the next
iteration.  The i-take and i-put instructions are used as semaphores to
synchronize between loop iterations to determine when a frame may be reused
for the next iteration.  The sync instructions gather signals from work
being performed in the current iteration to ensure that all work is
completed in the current iteration.  Additional fork instructions are used
to distribute values to the new instructions.

In addition to the per iteration costs shown in Figure
\ref{figure:loop-cycle-counts}, loops have startup and shutdown costs
that may be significant, depending upon the loop schema.  Completely
unrolled loops have no startup cost, whereas sequential loops must store
loop constants into the frame before beginning loop iterations.  The
startup and shutdown costs of $k$-bounded loops can be from about 100 to
several hundred cycles per $k$, to allocate and deallocate frames, to set
up loop constants and to wire up the frames in a cycle.

Given the widely varying costs of loop execution, even for the same source
program, it is clear that the programmer or compiler must make decisions to
expose loop parallelism while minimizing instruction counts.  We will
discuss these issues further in Section \ref{sss:loop-optimizations}.

\subsection{Programming in Id}

Id is a layered language whose core consists of a higher-order, statically
typed, polymorphic, non-strict functional language.  Id also has mutable
and synchronizing data structures called I-structures and M-structures.
The GCM code is written using only the functional subset of Id, but the Id
compiler transforms functional data structures into I-structures for
execution on Monsoon.  I-structures are implicitly synchronizing, so that a
fetch of an I-structure element will not return until the store to that
element is executed.  I-structure elements can be written only once
(although they can be read any number of times) to ensure that programs are
deterministic.

\subsubsection{Array Abstractions}
\label{sss:array-abstractions}

Id has only three built-in operations on arrays; fetch, store and fetch
array extents.  However, it is possible to define powerful array
abstractions in Id and use them repeatedly.  Before defining {\tt
make\_3D\_array}, an abstraction to create a three dimensional array, we
show how it is used:
\begin{code}
\small
\begin{verbatim}
 x = make_3D_array extents f
\end{verbatim}
\end{code}
where {\tt extents} defines the extents of the array and {\tt f} is a {\it
filling} function. The net effect is that {\tt x[index] = f index}, for
each index {\tt (k,j,i)} that lies within the {\tt extents}.  By giving a
suitable definition for {\tt f} we can define each GCM state variable.
Note that {\tt x} is {\em immutable} in the sense that once it has been
defined it cannot be modified.  Consequently, at each time step a new array
for each state variable is created and the old is implicitly discarded.

An (inlinable) function definition of 3D daxpy can be expressed as:
\begin{code}
\small
\begin{verbatim}
 Zones = 
    ((min_W_longitude, max_E_longitude),
     (min_N_parallel , max_S_parallel ),
     (surface        , max_depth      ));

 defsubst aX_plus_Y a X Y zone = 
   a * X[zone] + Y[zone];

 defsubst daxpy_3d a X Y =
   make_3D_array Zones (aX_plus_Y a X Y);
\end{verbatim}
\end{code}
where {\tt zone} represents an index and {\tt Zones} represents the
dimensions of the box containing the ocean.  One can think of {\tt
aX\_plus\_Y} as a inlinable function that takes four arguments and produces a
number.  When only three arguments are passed to it, as shown above, a new
function is returned that takes one argument {\tt zone} (an index) and
returns a number.

The {\tt make\_3D\_array} function is not a built-in operator in Id: it is
defined as a library function using {\it array comprehension} syntax as
follows:
\begin{code}
\small
\begin{verbatim}
 defsubst make_3D_array extents f =
   { ((kl,ku),(jl,ju),(il,iu)) = extents;
     in
    {3D_array extents of
       | [index] = f index
                  || i <- il to iu;
                     j <- jl to ju;
                     k <- kl to ku;
                     index = (k,j,i) }};
\end{verbatim}
\end{code}
The initial binding destructures the {\tt extents} argument into its
component parts.  The phrase between ``{\tt |}'' and ``{\tt ||}'', that is
``{\tt [index] = f index}'', specifies an index and an expression to be
filled in the indexed slot of the array.  The phrases such as ``{\tt i <-
il to iu}'' that occur after ``{\tt ||}'' are called {\em generators} and
define a way to enumerate the indices.

All of the abstraction involved in the array comprehension, destructuring
and higher-order function call could be extremely expensive, but if the
programmer is careful to declare the filling function and array constructor
as inlinable, all of the overhead is compiled down into an I-structure
array allocation and a triply-nested loop that fills the I-structure using
the inlined filling function.

This level of abstraction may seem to be overkill for a simple operation
like daxpy, but it is also used at every level of the code, simplifying the
specification of the algorithm greatly and also allowing for easier
modification to the computation by changing of the underlying abstraction.

Notice that, in the Id array constructor abstraction, we
ordered of the axes as we did in the CM Fortran code, although for a quite
different reason.
In Id, we want the innermost loop to iterate over the vertical
axis of the state arrays so that we can avoid computation on land, and also
because we want to unroll the loop completely.  Given that the innermost
loop iterates over the vertical axis, we also want that axis to be
stride-1, so that we can avoid overheads from additional support
instructions to calculate memory addresses.

\subsubsection{Avoiding unnecessary computation}

In the elliptic solver, we compute four inner products involving state
variables.  Because the land zones do not contribute to the inner product, we
could optimize this computation by only iterating over water zones as shown
below:
\begin{code}
\small
\begin{verbatim}
 defsubst ip_3d a b = 
   { s = 0.0 ;
    in
     {for i <- min_W_longitude to
               max_E_longitude do
       s_jk = 0.0;
       next s = s + 
        {for j <- northmost_parallel[i] to 
                  southmost_parallel[i] do
          s_k = 0.0;
          next s_jk = s_jk +
           {for k <- surface to
                     bottom[i,j] do
              next s_k = 
               s_k + a[k,j,i]*b[k,j,i];
              finally s_k};
          finally s_jk};
       finally s} }; 
\end{verbatim} 
\end{code} 
where {\tt bottom[i,j]} is the depth of the ocean at ({\tt i,j}).  The
variable {\tt surface} has been defined as a constant instead of a 2D array
without loss of generality.  {\tt northmost\_parallel[i]} and {\tt
southmost\_parallel[i]} represent, respectively, the northern most and
southern most parallel at longitude {\tt i} where a water zone may be
found.  A note on syntax -- the {\tt next} qualifier is used to express
recurrences in {\tt for} loops -- ``{\tt next s}'' specifies the value the
variable {\tt s} will assume in the next iteration.  The {\tt finally}
qualifier is used to describe the result of a {\tt for} loop.

We can further reduce the computation by defining some of the state
variables only on water zones, if their values are not needed on land
zones.  The {\tt make\_water\_zone\_array} abstraction defined below
enumerates only water zones:
\begin{code}
\small
\begin{verbatim}
 defsubst make_water_zone_array f =
  {3D_array Zones of
     | [zone] = f zone
       || i <- min_W_longitude to 
               max_E_longitude;
          j <- northmost_parallel[i] to 
               southmost_parallel[i];
          k <- surface to bottom[i,j];
          zone = (k,j,i) }
\end{verbatim}
\end{code}

Consequently the following daxpy would do fewer floating point operations
than the one defined in Section \ref{sss:array-abstractions}.
\begin{code}
\small 
\begin{verbatim}
 defsubst daxpy_w_3d a X Y =
   make_water_zone_array (aX_plus_Y a X Y);
\end{verbatim}
\end{code}

\subsubsection{Avoiding deadlocks due to undefined elements}

Sometimes an attempt to avoid unnecessary computations can lead to
deadlocks in the Id code.  Consider the stencil computation in the PCG (as
discussed in the data parallel section, the diagonals of matrix A are
symmetric; so only four arrays need to be stored instead of seven):
\begin{code}
\small
\begin{verbatim}
 defsubst seven_pt_stencil b =
   make_water_zone_array (multiply_by_A b);

 defsubst multiply_by_A b zone =
   { C = A_C[zone] * b[zone];
     W = A_WE[Wfz zone] * b[Wzz zone];
     E = A_WE[Efz zone] * b[Ezz zone];
     N = A_NS[Nfz zone] * b[Nzz zone];
     S = A_NS[Sfz zone] * b[Szz zone];
     U = A_UL[Ufz zone] * b[Uzz zone];
     L = A_UL[Lfz zone] * b[Lzz zone];
    in 
      C + W + E + N + S + U + L }};
\end{verbatim}
\end{code}
where the function {\tt Nzz} computes the zone to the North of a zone,
while {\tt Nfz} computes the north face of a zone:
\begin{code}
\small
\begin{verbatim}
 defsubst Wzz (k,j,i) = k,j,i-1 ;
 defsubst Ezz (k,j,i) = k,j,i+1 ;
 ...
 defsubst Wfz (k,j,i) = k,j,i-1 ;
 defsubst Efz (k,j,i) = k,j,i ;
 ...
\end{verbatim}
\end{code}

The {\tt multiply\_by\_A} filling function will attempt to read elements of
the array {\tt b} that may correspond to land zones, and thus may have
undefined values.  An attempt to read an undefined array element in Id will
cause a deadlock because the I-structure fetch will never return.  In the
CM Fortran version, the ocean is padded with a layer of land to handle this
boundary condition around the edges, and the other land zones are filled with
benign values.  The same strategy can be used in the Id version, however, the
land padding cannot simply be initialized once and re-used at every PCG
iteration. Rather, the padding must be filled every time
because a new array is defined for the variable on every iteration of the
solver.

Note that the land padding is not necessary for the non-stencil operations
in the PCG, so we only pad array {\tt b}, which is needed by the stencil
computation.  To accomplish this, we define a different array constructor
that fills zeros in all the land zones and use it to redefine {\tt
daxpy\_3d} of Section \ref{sss:array-abstractions}.  In the following code,
{\tt OG} represents the ocean geometry:
\begin{code}
\small
\begin{verbatim}
 def make_zone_array' f =
   {3D_array Zones of
      | [zone] = {case OG[zone] of
                   | water = f zone
                   | land  = 0.0}
        || i <- min_W_longitude to
                max_E_longitude;
           j <- min_N_parallel  to 
                max_S_parallel;
           k <- surface to max_depth;
           zone = (k,j,i) }; 

 defsubst daxpy_3d' a X Y =
   make_zone_array' (aX_plus_Y a X Y);
\end{verbatim} 
\end{code}

\subsubsection{Putting it all together}
\label{sss:mt-preconditioner}

An important step in the elliptic solver is the multiplication by the
preconditioner $\K$, which is actually implemented as $N_{x} \times N_{y}$
independent tridiagonal linear equation solvers of size $N_{z} \times
N_{z}$ as described in Section \ref{section:conjugate-gradient}.  As in the
CM Fortran code, each of the $N_{x} \times N_{y}$
independent linear systems is initially factored into $\Lh$ and $\Uh$ matrices 
corresponding to
the decomposition of the block diagonal elements of $\A$, and these
matrices are stored into a constant matrix, which are referenced in the call to {\tt
solve\_tridiag\_LU}.
\begin{code}
\small
\begin{verbatim}
 defsubst precondition r = 
  {3D_array Zones of
   | [zone] = X[k]
     || i <- min_W_longitude to
             max_E_longitude;
        j <- northmost_parallel[i] to
             southmost_parallel[i];
        X = solve_tridiag_LU r (i,j);
        k <- surface to bottom[i,j];
        zone = (k,j,i) };
\end{verbatim}
\end{code}

\begin{fig}{figure:pcg-id}{The PCG in Id requires the definition of new state
arrays on every iteration.  The order of execution is not necessarily in
the textual order -- all of the bindings proceed implicitly in parallel and
must only wait for data dependencies to be satisfied.}
\begin{center}
\small
\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_w_3d alpha p b;
    next r  = daxpy_w_3d (-alpha) r Ab;
    next xi = precondition (next r);
    beta    = (ip_3d (next r) (next xi)) / 
                gamma;
    next b  = daxpy_3d' beta (next xi) b;
  finally p}
\end{verbatim}
\end{center}
\end{fig}

Given all these functions, it is straightforward to write PCG in Id as
shown in Figure \ref{figure:pcg-id}.  There are no other issues that arise
in coding the whole GCM code; each filling function is either some sort
of a stencil, or involves purely local computation.

After compilation, if all the functions have been declared as
substitutable, the resulting code turns out to be nothing but a nest of
nested loops.  Most of these loops and the individual bindings in them can
proceed in parallel, as soon as the data dependencies for the binding are
satisfied.  Note that because of the element-wise synchronization of
I-structures, parts of state arrays can be filled as soon as the
parts of the arrays they depend on are filled.

\subsubsection{Programmer annotations and restructuring}
\label{sss:loop-optimizations}

Once the functionality of the code has been implemented, the programmer may
have to annotate his code to get good performance.  Standard loop
optimizations such as unrolling, peeling, strip-mining, and interchange can
have a larger impact on Monsoon performance than for conventional
architectures because of the overhead of parallel asynchronous execution of
loop iterations in the dataflow model.  The Id compiler automatically
performs some of these optimizations, and the others are left to the
programmer to be indicated as pragmas.  We discuss the role of loop pragmas
here.

Unless otherwise annotated by the programmer, all loops are assumed by the
compiler to be sequential loops.  Constructs such as array comprehensions
desugar into loops, and in general, annotations that are used for loops
may also be used for array comprehensions.  It is a good rule of thumb to
annotate {\it outermost} loops with parallelism to be $k$-bounded loops.
(It makes no sense to declare the time-step loop in GCM to be $k$ bounded
because there is not much parallelism between time-steps.)  In this way,
large chunks of work are forked to remote processors, and the overheads of
$k$-bounded loop iterations are paid less often.  Inner loops should be
annotated to be unrolled to do $n$ iterations at a time to reduce
sequential loop overhead.  If the number of iterations is a constant and
small, then the pragma to unroll the loop completely should be used.

Suppose that the programmer has annotated the outer and middle loops of a
triply-nested loop to be $k$-bounded, with the $k$-bounds being {\it ib}
and {\it jb}, respectively.  In general $k$-bounds can be arbitrary
expressions that return an integer; in practice, however, they are usually
literals.  Figure \ref{figure:loop-frames} shows the pattern of frame usage
for the two loops when {\it ib} is 3 and {\it jb} is 4.

\begin{fig}{figure:loop-frames}{Pattern of frame usage for nested
$k$-bounded loops for the 3D inner product.}
\illustrationx{3in}{figures/loop-frames.eps}
\end{fig}

The amount of interprocessor parallelism exposed for such a nested loop is
then on the order of $ib \times jb$.  Given a certain machine
configuration, we are interested in setting the product $ib \times jb$ such
that it will keep the entire machine busy.  What are the optimal values of
$ib$ and $jb$, such that the number of instructions executed by the loops
are minimized?

If there are $imax$ iterations of the outer loop executed, then the number
of frame allocations and initializations (requiring hundreds of cycles
apiece) is then: $ib + imax \times jb$.  For a set amount of parallelism
$ib \times jb$ we wish to exploit in a nested loop, in general it is better
to make $ib$ larger and $jb$ smaller to minimize loop initialization costs.
A possible optimization is to lift middle loop frame allocations so that it
occurs $ib$ times as opposed to $imax$ times.  This optimization has been
implemented in the Id compiler, but has not worked reliably.

Large values of $k$ may not be useful in exposing parallelism if there are
dependencies of any kind (data, control, producer/consumer) between loop
iterations, or if the number of iterations is not significantly larger than
$k$.  Choosing loop bounds can make a tremendous difference in performance,
as shown by Culler \cite{Culler90}, but to date, no automatic
compiler-directed policy has been implemented that achieves high
performance.

Some optimizations cannot be expressed simply with pragmas and require
restructuring of code.  For example, the precondition routine in Section
\ref{sss:mt-preconditioner} requires allocating a temporary array {\tt X}
for each column {\tt (i,j)} of the ocean for each iteration of PCG.  The
only purpose of this array is to hold the values returned by the
tridiagonal solver before it is copied into the state array. 
This array is eliminated by making the tridiagonal solver store directly into array {\tt
xi}.  Such a program can be expressed in Id by writing the preconditioner
using loops and I-structures rather than array comprehensions.  A subtle
issue is if the ocean has a known uniform depth, then it is also possible
to eliminate another array temporary in the tridiagonal solver.

\subsection{Running GCM on Monsoon hardware}

Monsoon was designed at MIT and built by Motorola under a research
collaboration agreement.  We received the first Monsoon with one PE and one I
S in November 1990.  Two machines, each with 8 PE's and 8 IS's, were
delivered to MIT and Los Alamos National Labs (LANL) in the latter half of
1991.  Several more machines with 2 PE's and 2 IS's were built and placed
at various research institutions across the country.  Detailed performance
studies of Monsoon were conducted in 1991-92 (see, for example,
\cite{id-performance,boon-loops,Sur94}).  During this period, both the Id
compiler and run-time system were being tuned continually for Monsoon.

By 1993, when the GCM study was begun, the team that had built Monsoon at
Motorola had disbanded; the machine at MIT was maintained by
cannibalizing parts from other Monsoons.  At this point, an initial
unoptimized version of a fully functioning multithreaded GCM ran on an 8 PE
8 IS Monsoon approximately 100 times slower than the data parallel version
on a 32-node (128 VU) CM-5.  Given the speed difference, experiments took
many hours to run on Monsoon, and checkpointing was implemented because the
Monsoon hardware had intermittent failures.  We started doing most of the
experiments on PCG and 2 PE 2 IS systems to avoid some of these problems.
Our experience with other applications on Monsoon \cite{id-performance}
have shown that 2 PE 2 IS performance translates directly into 8 PE 8 IS
performance. 

As is often the case in any such study, both implementations were
continually being tuned based on observed performance.  However,
algorithmic innovations were always guided by the concerns for greater
efficiency on the CM-5, as this was a production code.  It was a difficult
task to keep the two implementations consistent for this study, especially
because gathering statistics on a 4 to 24 hour run on Monsoon was tedious.
A major milestone in this project was reached in early 1994 when we were
able to produce the same numerical results on both implementations for
several large ocean geometries.

The hardware maintenance issue became more serious in 1995 because key
personnel at MIT were concentrating on building new machines, and the
condition of Monsoon had deteriorated.  Most of the experiments from 1995
onwards had to be performed 2 PE 2 IS machine.  No new experiments have
been conducted on Monsoon since mid-1996.  Monsoon was decommissioned and
donated to the Computer Museum in February 1997, and the CM-5 was
decommissioned in July 1997.

Although we also have a robust cycle-level Monsoon simulator, even a 2 PE 2
IS Monsoon is many thousand times faster.  The Monsoon simulator was used
extensively for micro-benchmarking, and can run several iterations of PCG
with small data sets in a reasonable amount of time.  However, running
realistic data sets with the entire GCM code would be impossible on the
simulator, whereas it is feasible on even a 2 PE 2 IS Monsoon.

\subsection{Multithreaded GCM performance}

To factor out machine size and clock rates, as in the data parallel case,
we measure the efficiency of the multithreaded implementation of GCM by
determining the number of cycles spent for each floating point operation.
The ``Water Only'' bar in Figure \ref{figure:gcm-comparison} shows the
breakdown of the overhead operations executed in the PCG code for a particular
geometry, a $171 \times 93 \times 5$ ocean with 60\% water.

\begin{fig}{figure:gcm-comparison}{Comparison of the three versions of 
PCG for a $171 \times 93 \times 5$ geometry, 60\% water.  ($ib=5$,
$jb=3$, 2 PE, 2 IS Monsoon)}
\illustrationx{2.9in}{figures/three-versions.ps}
\end{fig}

It is clear from this graph that the vast majority of the time is spent in
overhead cycles.  ``Move'' instructions include forks and syncs, ``Misc''
instructions include switches and message-passing instructions, ``Memory''
instructions include i-fetches and i-stores, and ``Int'' instructions
include ordinary arithmetic and logical instructions.  Actual floating
point operations consist of a tiny fraction of overall time.

\subsubsection{Decreasing overhead cycles}

To reduce overhead cycles, we attempted to simplify the computation by
making it identical to the data parallel implementation.  In doing so, we
were able to completely unroll the inner loops of each of the PCG
components.  Complete loop unrolling is beneficial to 3D daxpy and 3D inner
product because the loop overheads are high relative to the small loop
bodies for these PCG components.  It does not help the 7-point stencil
computation as much because of the larger loop body of the stencil.  The
preconditioner is significantly simplified because of the elimination of a
temporary array.

I-fetches for the north-south and bottom bounds of the
ocean were eliminated, along with the calculation of the addresses for those i-fetches.  By making
the ocean bounds constant, some loop constants can be eliminated by making
them literals.  These simplifications reduce the number of cycles per
floating point operation from about 35 in the ``Water Only'' case to about
12 in the ``Land$+$Water'' case in Figure \ref{figure:gcm-comparison}.

Note that the ``Idle'' cycle count is misleading in that the $k$-bounds for
the loops are set in the context of the entire application, not for the PCG
executing in isolation.  Those idle cycles will be overlapped by other work
when PCG executes within the entire GCM code.

The reduction in overhead cycles per floating point operation comes from
fewer switches and messages for simpler loops, and from fewer bubbles and
forks incurred by eliminated switches and messages.  Address calculation is
also simpler, as evidenced by fewer ``Int'' instructions, and some fewer
memory operations are executed because of the constant ocean bounds.

It is surprising to discover the degree to which the
``Land$+$Water'' version was more efficient than the ``Water Only''
version.  Part of the explanation is in the very shallow geometry used;
 because the innermost loop executes at most four times in the ``Water
Only'' case, it is difficult to take advantage of regular loop unrolling
(as opposed to complete loop unrolling) to reduce loop overheads.  For
geometries with more layers, the difference will not be so marked, although
the completely unrolled preconditioner would remain an advantage of the
``Water$+$Land'' version.

\subsubsection{Hybrid PCG}

Although the ``Land$+$Water'' version is much more efficient per floating
point operation executed than the ``Water Only'' version, it must execute
more floating point operations because of land zones.  Despite this, the
overheads from land are more than made up for by the efficiency of the
``Land$+$Water'' version for most geometries; however, we would still
like to take advantage of some of the geometry irregularity that result from land
zones, while also using a simplified loop structure.

To do this, a {\it hybrid} version of the PCG was developed 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, much of the computation on land is eliminated, while still
allowing complete unrolling of the innermost loop.  Figure
\ref{figure:gcm-comparison} shows that this version has more overhead per
floating point operation than the ``Land$+$Water'' version, but less than
the ``Water Only'' version.

The hybrid approach is only superior to the Land$+$Water approach in
situations where it can take advantage of not having to compute on land.
Typically, geometries that we are interested in simulating have between
30\% to 70\% of their surface covered with water, so the hybrid approach
actually does pay off most of the time.

\subsubsection{PCG component performance and GCM performance}

\begin{fig}{figure:operation-comparison}{The best case performance
for the components of the hybrid PCG, the PCG itself, and the whole GCM
ocean modeling code for a $32 \times 32 \times 32$ ocean with 100\% water.
($ib=5$, $jb=3$, 2 PE, 2 IS Monsoon).}
\illustrationx{3.1in}{figures/operation-comparison.ps}
\end{fig}

Figure~\ref{figure:operation-comparison} shows the best case performance in
terms of 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.  In this
case, the same criteria as the data parallel
case is used; it is assumed that all floating point operations executed are required
operations, without taking into account land.

The problem size we use is proportionately smaller than the CM-5 version
because the Monsoon machine configuration is smaller.  Note that the Id
program is the same for a 2 PE Monsoon as for an 8 PE Monsoon, except for
$k$-bounds for loops.

We optimized the $k$-bounds shown in Figure \ref{figure:gcm-comparison} and
Figure \ref{figure:operation-comparison} for the entire GCM computation.
By doing so, when we extract portions of the code such as the PCG
computation or components of the PCG, additional idle cycles appear because
there is less work that can be overlapped to eliminate idle cycles.  These
idle cycles are misleading, because they could be eliminated by increasing
$k$-bounds, but doing so could distort other numbers in the overhead
accounting.

It is interesting to compare the performance of the 3D inner product in
Figure \ref{figure:operation-comparison} with the detailed cycle count for
the completely unrolled inner loop in Figure
\ref{figure:loop-cycle-counts}.  In Figure
\ref{figure:loop-cycle-counts}, eight cycles are spent for two floating point
operations, yielding a ratio of four cycles/floating point operation.  In
Figure \ref{figure:operation-comparison}, the full triply-nested loop shows
a ratio of about eight cycles per floating point operation, indicating that much of
the overhead arises from the outer two loops of the inner product.  This is
similar to the CM-5, where the detailed cycle count for the innermost loop
was much better than the actual measured time of the entire 3D inner
product.

Compared to the best case numbers for the CM-5, the Monsoon numbers show
that approximately three times as many cycles are necessary for every
floating point operation.  However, our choice of performance measure
showcases the best aspects of data parallel computing.  In the
next section, we attempt to adjust these raw ``best case'' numbers for
situations that are more realistic.
