A Multigrid Application

The multigrid method [5] is a fast algorithm for solution of linear and nonlinear problems. It uses a hierarchy or stack of grids of different granularity (typically with a geometric progression of grid-spacings, increasing by a factor of two up from finest to coarsest grid). Applied to a basic relaxation method, for example, multigrid hugely accelerates elimination of the residual by restricting a smoothed version of the error term to a coarser grid, computing a correction term on the coarse grid, then interpolating this term back to the original fine grid. Because computation of the correction term on the fine grid can itself be handled as a relaxation problem, the strategy can be applied recursively all the way up the stack of grids.

In our example, we apply the multigrid scheme to solution of the
two-dimensional Poisson equation. For the basic, unaccelerated, solution
scheme we use red-black relaxation. An HPJava method for red-black
relaxation is given in Figure 1.
This looks something like an HPF program with different syntax. One
obvious difference is that the base language is Java instead of Fortran.
The HPJava type signature `double [[-,-]]`
means a two dimensional distributed array of `double` numbers^{}. So the
arguments passed to the method `relax()` will be distributed arrays

The inquiry `rng()` on the distributed array `f` returns the
`Range` objects `x`, `y`. These describe the distribution
format of the array index (for the two dimensions).

The HPJava `overall` construct operates like a `forall` construct,
with one important difference. In the HPJava construct one *must*
specify how the iteration space of the parallel loop is distributed over
processors. This is done by specifying a `Range` object in
the header of the construct.

The variables `i`, `j` in the figure are called *distributed
index* symbols. Distributed indexes are scoped by the overall constructs
that use them. They are *not* integer variables, and there
is no syntax to declare a distributed index *except* through an overall
construct (or an *at* construct--see later). The usual Java
scoping rules for local variables apply: one can't for example use
`i` as the index of an overall if there is already a local variable
`i` in scope--the compiler doesn't allow it.

An unusual feature of the HPJava programming model is
that the subscripts in a distributed array element reference usually
*must* be distributed index symbols.
And these symbols must be distributed
with the essentially *same* format as the arrays they subscript.
As a special case, *shifted* index expressions like `i+1`
are allowed as subscripts, but only if the distributed array was
created with *ghost regions*. Information on ghost regions, along
with other information about distribution format, is captured in the
`Range` object associated with the array dimension or index.

These requirements ensure that a subscripting operation in an
overall construct only accesses locally held elements. They
place very stringent limitations on what kind of expression can
appear as a subscript of a distributed array.
We justify this by noting that this restricted kind of data parallel
loop is a frequently recurring pattern in SPMD programs in general, and it is
convenient to have it captured in syntax. A glance at the full source of
the applications described in this paper should make this claim more
plausible^{}.

The method `Adlib.writeHalo()` is a communication method (from the
library called Adlib). It performs the edge-exchange to fill in the
ghost regions. As emphasized earlier, the compiler is not responsible
for inserting communications--this is the programmer's responsibility.
We assume this should be acceptable to programmers currently accustomed
to using MPI and similar libraries for communication.

Because of the special role of distributed index symbols in array subscripts,
it is best not to think of the expressions `i`, `j`, `i+1`,
etc, as having a numeric value: instead they are treated as a special kind
of thing in the language. We use the notation `i``
to extract the numeric global index associated with `i`, say^{}.
In particular, use of this expression in the modulo 2 expression in
the inner overall construct in Figure 1
implements the red-black pattern of accesses.

This completes the description of most ``non-obvious'' features of HPJava syntax. Remaining examples in the paper either recycle these basic ideas, or just introduce new library routines; or they import relatively uncontroversial syntax, like a syntax for array sections.

Figure 2 visualizes
the ``restrict'' operation
that is used to transfer the error term from a
fine grid to a coarse grid.
The HPJava code is given in Figure 3.
The restrict operation here computes
the residual term at sites of the fine grid with *even* coordinate
values, then sends these values to the coarse grid.
In multigrid the restricted residual from the fine grid becomes
the RHS of a new equation on the coarse grid.
The implementation uses a temporary
array `tf` which should be aligned with the fine grid (only a subset
of elements of this array are actually used). The last line introduces
two new features: distributed array *sections*, and the library
function `Adlib.remap()`.
Sections work in HPJava in much the
same way as in Fortran--one small syntactic difference is that they use double
brackets. The bounds in the `fc` section ensure that edge values,
corresponding to boundary conditions, are not modified. The stride in the
`tf` section ensures only values with even subscripts are selected.
The `Adlib.remap()` operation is needed because in general there
is no simple relation between the distribution format of the fine and coarse
grid--this function introduces the communications necessary to perform
an assignment between any two distributed arrays with unrelated distribution
format.
As another example, the interpolation code of Figure 4
performs the complementary transformation from the coarse grid to the
fine grid.

The basic pattern here depends only on the geometry of the problem. More complex (perhaps non-linear) equations with similar geometry could be tackled by similar code. Problems with more dimensions can also be programmed in a similar way.