

\documentstyle[psfig]{article}

\begin{document}

\sloppy

\pagestyle{plain}

\title{\vspace{2in}Introduction to Java-Ad}

\author{Bryan Carpenter, Guansong Zhang, Geoffrey Fox \\
        Xinying Li, Yuhong Wen\vspace{5mm} \\
        \em NPAC at Syracuse University \\
        \em Syracuse, NY 13244 \\
        \em \{dbc,zgs,gcf,xli,wen\}\@npac.syr.edu}

\date{November 14, 1997}

\maketitle

\begin{abstract}
We outline an extension to Java for programming
with distributed arrays.  The basic programming style is Single Program
Multiple Data (SPMD), but parallel arrays are provided as new language
primitives.  Further extensions include three {\em distributed control}
constructs, the most important being a data-parallel loop construct.
Communications involving distributed arrays are handled through a
standard library of collective operations.  Because the underlying
programming model is SPMD programming, direct calls to MPI or other
communication packages are also allowed in a Java-Ad program.
\end{abstract}

\newpage

\section{Introduction}

The Java-Ad programming model is a flexible hybrid of the
SPMD\footnote{Single Program, Multiple Data} and data parallel
approaches to parallel programming.  It provides HPF-like distributed
arrays as language primitives, and new {\em distributed control}
constructs to facilitate access to the local elements of these arrays.
In the SPMD mold, the model allows processors the freedom to
independently execute complex procedures on local elements: it is not
restricted by SIMD-style array syntax.

In Java-Ad, all access to {\em non-local} array elements must go through
library functions---typically collective communication operations.
This puts an extra onus on the programmer; but making communication
explicit encourages the programmer to write algorithms that exploit
locality, and simplifies the task of the compiler writer.
On the other hand, by providing distributed arrays as language
primitives we are able to simplify error-prone tasks
such as converting between local and global array subscripts and
determining which processor holds a particular element.  In the Java-Ad
model the programmer never has to deal explicitly with local array
subscripts.  As in HPF, it is possible to write programs at a natural level of
abstraction where the meaning is insensitive to the detailed mapping of
elements.  Lower-level styles of programming are also possible.

Our Java-Ad compiler will be implemented as a translator
to ordinary Java with calls to a suitable run-time library.
At the time of writing the underlying library is already available, and
the interface needed by the Java-Ad translator is under development.
The translator itself is being implemented in a compiler construction
framework developed in the PCRC project.

\section{Multidimensional arrays\label{multiDim}}

First we describe a modest extension of Java that would add a class of true
multi-dimensional arrays to the standard Java language.  The new arrays
allow regular section subscripting, similar to Fortran 90 arrays.  The
syntax described in this section is a subset of the Java-Ad syntax for
parallel arrays and algorithms: the motivation for introducing the sequential
subset first is just to simplify the overall presentation.

No attempt is made to integrate the new multidimensional arrays
with the standard Java arrays: they are a new kind of entity that
coexists in the language with ordinary Java arrays.  There
are good technical reasons for keeping the two kinds of array
separate\footnote{The run-time representation of our multi-dimensional arrays
includes extra descriptor information that would simply
encumber the large class ``non-scientific'' Java applications.}.

The type-signatures and constructors of the multidimensional array use
double brackets to distinguish them from ordinary arrays:
\small
\begin{verbatim}
  int [[,]] a = new int [[5, 5]] ;

  float [[,,]] b = new float [[10, n, 20]] ;

  int [[]] c = new int [[100]] ;
\end{verbatim}
\normalsize
{\tt a}, {\tt b} and {\tt c} are respectively 2-, 3- and one- dimensional
arrays.  Of course {\tt c} is very similar in structure to the standard
array {\tt d}, created by
\small
\begin{verbatim}
  int [] d = new int [100] ;
\end{verbatim}
\normalsize
{\tt c} and {\tt d} are not identical, though\footnote{For example,
{\tt c} allows section subscripting, whereas {\tt d} does not.}.

Accessing individual elements of a multidimensional array goes through
a subscripting operation involving single brackets
\small
\begin{verbatim}
  for(int i = 0 ; i < 4 ; i++)
    a [i, i + 1] = i + c [i] ;
\end{verbatim}
\normalsize
For reasons that will become clearer in later sections, this style of
subscripting is called {\em local subscripting}.  In the current
sequential context, apart from the fact that a single pair of brackest
may include several comma-separated subscripts, this kind of
subscripting behaves just like ordinary Java array subscripting.
Subscripts always start at zero, in the ordinary Java or C style (there
is no Fortran-like lower bound).

In general Java-Ad has no idea of Fortran-like array assignments.  In
\small
\begin{verbatim}
  int [[,]] e = new int [[n, m]] ;
  ...
  a = e ;
\end{verbatim}
\normalsize
the assignment simply copies a handle to object referenced by {\tt e}
into {\tt a}.  There is no element-by-element copy involved.  Similarly Java-Ad
has no idea of elemental arithmetic or elemental function application.  If
{\tt e} and {\tt a} are arrays, the expressions
\small
\begin{verbatim}
  e + a
  Math.cos(e)
\end{verbatim}
\normalsize
are type errors.

Java-Ad {\em does} import a Fortran-90-like idea of array {\em regular
sections}.  The syntax for {\em section subscripting} is different to
the syntax for local subscripting.  Double brackets are used.  These
brackets can include scalar subscripts or {\em subscript triplets}.

A section is an object in its own right---its type is that of a suitable
multi-dimensional array.  It describes some subset of the elements of
the parent array.  This is slightly different to the situation in
Fortran, where sections cannot usually be captured as named
entities\footnote{Unless a section appears as an actual argument to
a procedure, in
which case the dummy argument names that section, or it is the target
of a pointer assignment.}.
\small
\begin{verbatim}
  int [[]] e = a [[2,  2 :]] ;

  foo(b [[ : , 0, 1 : 10 : 2]]) ;
\end{verbatim}
\normalsize
{\tt e} becomes an alias for the 3rd row of elements of {\tt a}.
The procedure {\tt foo} accepts a two-dimensional array as argument.
It can read or write to the set of elements of {\tt b} selected
by the section.
As in Fortran, upper or lower bounds can be omitted in triplets,
defaulting to the actual bound of the parent array, and the
stride entry of the triplet is optional.
Note that the subscripts of {\tt e}, like any other array, start at 0,
although the first element is identified with {\tt a [2, 2]}.

In Java-Ad, unlike Fortran, it is not allowed to use vectors of
integers as subscripts.  The only sections
recognized are regular sections defined through scalar and triplet
subscripts.

Java-Ad provides a library of functions for manipulating its
arrays, closely analogous to the array transformational intrinsic
functions of Fortran 90:
\small
\begin{verbatim}
  int [[,]] f = new int [[5, 5]] ;
  Adlib.shift(f, a, -1, 0, CYCL) ;

  float g = Adlib.sum(b) ;

  int [[]] h = new int [[100]] ;
  Adlib.copy(h, c) ;
\end{verbatim}
\normalsize
The {\tt shift} operation with shift-mode {\tt CYCL} executes a cyclic
shift on the data in its second argument, copying the result to its
first argument---an array of the same shape.  In the example
the shift amount is -1, and the shift is performed in dimension 0 of
the array---ie, the first of its two dimensions.  The {\tt sum} operation
simply adds all elements of its argument array.  The {\tt copy} operation
copies the elements of its second argument to its first---it is something
like an array assignment.  These functions may have to be overloaded
to apply to some finite set of array types, eg they may be defined for
arrays with elements of any suitable Java primitive type, up to some maximum
rank of array.  Alternatively the type-hierarchy for arrays may be defined
in a way that allows these functions to be more polymorphic.

\section{Process arrays\label{procArrays}}

Java-Ad adds class libraries and some additional syntax for dealing
with {\em distributed arrays}.  These arrays are viewed
coherent global entities, but their elements are divided
across a set of cooperating processes.  As a pre-requisite to
introducing distributed arrays we discuss the {\em process arrays} over
which their elements are scattered.

An abstract base class {\tt Procs}
has subclasses {\tt Procs1}, {\tt Procs2}, \ldots, representing
one-dimensional process arrays, two-dimensional process arrays, 
and so on.
\small
\begin{verbatim}
  Procs2 p = new Procs2(2, 2) ;
  Procs1 q = new Procs1(4) ;
\end{verbatim}
\normalsize
These declarations set {\tt p} to represent a 2 by 2 process array and {\tt q}
to represent a 4-element, one-dimensional process array.  In either case
the object created describes a group of 4 processes.
At the time the {\tt Procs} constructors are executed the program should be
executing on four or more processes.  Either constructor selects four
processes from this set and identifies them as members of
the constructed group\footnote{There is no cooperation between
the two constructor calls for {\tt p} and {\tt q}, so an individual
physical process might occur in both groups or in
neither.  As an option not illustrated here, vectors of
ids can be passed to the {\tt Procs} constructors to specify
exactly which processes are included in a particular group.}.

{\tt Procs} has a member function called {\tt member}, returning
a boolean value.  This is {\tt true} if the local process is a member
of the group, {\tt false} otherwise.
\small
\begin{verbatim}
  if(p.member()) {
    ...
  }
\end{verbatim}
\normalsize
The code inside the {\tt if} is executed only if the
local process is a member {\tt p}.  We will say that inside this
construct the {\em active process group} is restricted to {\tt p}.

The multi-dimensional structure of a process array is reflected in
its set of {\em process dimensions}.  An object is associated
with each dimension.  These objects are accessed through the
inquiry member {\tt dim}:
\small
\begin{verbatim}
  Dimension x = p.dim(0) ;
  Dimension y = p.dim(1) ;

  Dimension z = q.dim(0) ;
\end{verbatim}
\normalsize
The object returned by the {\tt dim} inquiry has class {\tt Dimension}.
The members of this class include the inquiry {\tt crd}.  This
returns the coordinate of the local process with respect to the
process dimension.  The result is only well-defined if the local
process is a member of the parent process array.  The inner
body code in
\small
\begin{verbatim}
  if(p.member())
    if(x.crd() == 0)
      if(y.crd() == 0) {
        ...
      }
\end{verbatim}
\normalsize
will only execute on the first process from {\tt p}, with coordinates
$(0, 0)$.


\section{Distributed arrays\label{distArrays}}

Some or all of the dimensions of a multi-dimensional array can be
declared to be {\em distributed ranges}.  In general a distributed
range is represented by an object of class {\tt Range}.  A {\tt Range}
object defines a range of integer subscripts, and defines how they are mapped
into a process array dimension.  In fact the {\tt Dimension}
class introduced in the previous section is a subclass of {\tt
Range}.  In this case the integer range is just the range of coordinate
values associated with the dimension.  Each value in the range
is mapped, of course, to the process (or slice of processes) with
that coordinate.  This kind of range is also called a {\em primitive range}.
More complex subclasses of {\tt Range} implement more
elaborate maps from integer ranges to process dimensions.
Some of these will be introduced in later sections.  For now we concentrate
on arrays
constructed with {\tt Dimension} objects as their distributed ranges.

The syntax of section \ref{multiDim} is extended in the following way
to support distributed arrays
\begin{itemize}
\item
A distributed range object may appear in place of an integer extent in the
``constructor'' of the array (the expression following the {\tt new}
keyword).
\item
If a particular dimension of the array has a distributed range, the
corresponding slot in the type signature of the array should include
a {\tt \#} symbol.
\item
In general the constructor of the distributed
array must be followed by an {\tt on} clause, specifying the process
group over which the array is distributed.  Distributed ranges
of the array must be distributed over distinct dimensions of this
group\footnote{The {\tt on} clause can be omitted in some circumstances---see
section \ref{on}.}.
\end{itemize}
Assume {\tt p},
{\tt x} and {\tt y} are declared as in the previous section, then
\small
\begin{verbatim}
  float [[#,#,]] a = new float [[x, y, 100]] on p ;
\end{verbatim}
\normalsize
defines {\tt a} as a 2 by 2 by 100 array of floating point
numbers.  Because the first two dimensions of the array are distributed
ranges---dimensions of {\tt p}---{\tt a} is actually realized as four
segments of 100 elements, one in each of the processes of {\tt p}.
The process in {\tt p} with coordinates {\tt i}, {\tt j}
holds the section {\tt a [[i, j, :]]}.

The distribute array {\tt a} is equivalent in terms of storage
to four local arrays defined by
\small
\begin{verbatim}
  float [] b = new float [100] ;
\end{verbatim}
\normalsize
But because {\tt a} is declared as a collective object we can apply
collective operations to it.  The {\tt Adlib} functions introduced in
section \ref{multiDim} apply equally well to distributed arrays, but
now they imply inter-processor communication.
\small
\begin{verbatim}
  float [[#,#,]] a = new float [[x, y, 100]] on p,
                 b = new float [[x, y, 100]] on p ;

  Adlib.shift(a, b, -1, 0, CYCL) ;
\end{verbatim}
\normalsize
The {\tt shift} operation causes the local values of {\tt a} is overwritten
with the values of {\tt b} from a processor adjacent in the {\tt x} dimension.

There is a catch in this.
When subscripting the distributed dimensions of an array it
is {\em simply disallowed} to use subscripts that refer to
off-processor elements.  While this:
\small
\begin{verbatim}
  int i = x.crd(), j = y.crd() ;

  a [i, j, 20] = a [i, j, 21] ;
\end{verbatim}
\normalsize
is allowed, this:
\small
\begin{verbatim}
  int i = x.crd(), j = y.crd() ;

  a [i, j, 20] = b [(i + 1) % 2, j, 20] ;
\end{verbatim}
\normalsize
is forbidden.  The second example could apparently be implemented using
a nearest neighbour communication, quite similar to the {\tt shift}
example above.  But Java-Ad imposes an strict policy
distinguishing it from many data parallel languages: while
library functions may introduce communications, language
primitives such as array subscripting {\em never} imply
communication.

If subscripting distributed dimensions is so restricted, why are the
{\tt i}, {\tt j} subscripts on the arrays needed at all?  In the
examples of this section these subscripts are only allowed one value on each
processor.
Well, the inconvience of specifying the subscripts will be reduced by
language constructs introduced later, and the fact that only one subscript
value is local is a special feature of the {\em primitive ranges}
used here.  The higher level distributed ranges introduced later
map multiple elements to individual processes.  Subscripting will no
longer look so redundant.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs1 p = new Procs1(P) ;
  if(p.member()) {
    Range x = p.dim(0) ;

    float a [[#,,]] = new float [[x, B, N]] on p ;
    float b [[#,,]] = new float [[x, N, B]] on p ;

    ... initialize `a', `b'

    float c [[#,,]] = new float [[x, B, N]] on p ;

    for(int s = 0 ; s < P ; s++) {

      const int ip = x.crd() ;
      const int base = B * ((ip + s) % P) ;

      // c [[ip, :, base : ...]] =
      //         a [[ip, :, :]] * b' [[(ip + s) % P, :, :]] ...

      for(int ib = 0 ; ib < B ; ib++)
        for(int kb = 0 ; kb < B ; kb++) {

          float sum = 0 ;
          for(int j = 0 ; j < N ; j++)
            sum += a [ip, ib, j] * b [ip, j, kb] ;

          c [ip, ib, base + kb] = sum ;
        }

      float tmp [[#,,]] = new float [[x, N, B]] on p ;

      Adlib.shift(tmp, b, 1, 0, CYCL) ;
      Adlib.copy(b, tmp) ;
    }
  }
\end{verbatim}
\normalsize
\caption{\label{matmul1d}A parallel matrix multiplication program.}
\end{figure}

We finish this section with a fairly complex example using
the notation established so far.  The algorithm of figure
\ref{matmul1d} implements multiplication of two {\tt N} $\times$ {\tt
N} matrices.  One dimension of each of the two matrices is
block-distributed over the {\tt P} processors of {\tt p}, so {\tt N} is
equal to {\tt P} $\times$ {\tt B} where {\tt B} is the (constant)
local block size.

The matrices are represented as three dimensional arrays, with their
distributed dimensions explicitly split into a distributed range of
extent {\tt P} and a local sequential range of extent {\tt B}.  In
later sections we will see how to represent this distribution format
with a single block-distributed {\tt Range} object.  Even with that
facility available, the representation used here may still be more natural
for algorithms like the current one, where the block structure is
an integral to the algorithm.  The undistributed dimensions of
the matrices are just sequential ranges of extent {\tt N}.  The
operation of the algorithm for {\tt P} $= 2$ is visualized in figure
\ref{matmul1dpic}.  There are two phases.  Between the phases the data
in {\tt b} is exchanged by the {\tt shift} operation\footnote{In fact
it is neccessary to use a {\tt shift} and {\tt copy} operation because
the source and destination arguments of {\tt shift} must be distinct
arrays.  In the comment explaining the inner block matrix multiplication,
by the way, the symbol {\tt b'} means the original unshifted value of 
the array {\tt b}.}.

\begin{figure}[btp]
\centerline{\psfig{figure=matmul1d.eps,height=1.8in}}
\caption{\label{matmul1dpic}Operation of the program in figure \ref{matmul1d}
for {\tt P} $= 2$.}
\end{figure}

\section{The {\em on} construct and the active process group\label{on}}

In the last two section the idiom 
\small
\begin{verbatim}
  if(p.member()) {
    ...
  }
\end{verbatim}
\normalsize
has appeared.  
Our language provides a short way of writing
this construct
\small
\begin{verbatim}
  on(p) {
    ...
  }
\end{verbatim}
\normalsize
In fact the {\em on} construct provides some extra value.
Informally we said in section \ref{procArrays} that the {\em active
process group} is restricted to {\tt p} inside the body of the {\tt
p.member()} conditional construct.  As part of the language, Java-Ad
includes a more formal idea of an active process group (APG).  At any
point of execution some process group is singled out as the APG.  An
{\tt on(p)} construct specifically changes the value of the APG to {\tt
p}.  On exit from the construct, the APG is restored to its value on
entry.

Elevating the active process group to a part of the language
allows some simplifications.  For example, it provides a natural
default for the {\tt on} clause in array constructors.
In the matrix multiplication program of the previous section the code
\small
\begin{verbatim}
  if(p.member()) {
    ...
    float a [[#,,]] = new float [[x, B, N]] on p ;
    float b [[#,,]] = new float [[x, N, B]] on p ;
    ...
  }
\end{verbatim}
\normalsize
can be simplified to
\small
\begin{verbatim}
  on(p) {
    ...
    float a [[#,,]] = new float [[x, B, N]] ;
    float b [[#,,]] = new float [[x, N, B]] ;
    ...
  }
\end{verbatim}
\normalsize
More importantly, formally defining the active process group will
simplify the statement of various rules about what operations are legal
{\em inside} distributed control constructs like {\em on}.

\section{Higher-level ranges and locations}

The class {\tt BlockRange} is a subclass of {\tt Range} which describes
a simple block-distributed range of subscripts.  Like {\tt BLOCK}
distribution format in HPF, it maps blocks of contiguous subscripts
to each element of its target process dimension\footnote{Other higher-level
ranges include {\tt CyclicRange}, which produces the equivalent
of {\tt CYCLIC} distribution formation in HPF.}.  The constructor
of {\tt BlockRange} usually takes two arguments: the extent of the range
and a {\tt Dimension} object defining the process dimension over
which the new range is distributed.
\small
\begin{verbatim}
  Procs2 p = new Procs2(3, 2) ;
  
  Range x = new BlockRange(100, p.dim(0)) ;
  Range y = new BlockRange(200, p.dim(1)) ;

  float [[#,#]] a = new float [[x, y]] on p ;
\end{verbatim}
\normalsize
{\tt a} is created as a 100 $\times$ 200 array, block-distributed in
over the 6 processes in {\tt p}.  The fragment is
essentially equivalent to the HPF declarations
\small
\begin{verbatim}
  !HPF$ PROCESSORS p(3, 2)
  
  REAL a(100, 200)

  !HPF$ DISTRIBUTE a(BLOCK, BLOCK) ONTO p
\end{verbatim}
\normalsize
Subscripting distributed arrays with non-primitive ranges
introduces some new problems.  An array access such as
\small
\begin{verbatim}
  a [17, 23] = 13 ;
\end{verbatim}
\normalsize
is perfectly legal {\em if} the local process holds the element in question.
But deterimining whether an element is local is not so easy any more.
When arrays had only {\em primitive} distributed ranges, it
was straightforward to check that accesses were local---the
subscript simply had to be equal to the local coordinate.  With
higher-level ranges, that simple condition no longer holds.

In practise it is unusual to use integer values directly as local
subscripts.  Instead the idea of a {\em location} is introduced.  A
location can be viewed as an abstract element, or ``slot'', of a
distributed range.
Conversely, a range can be thought of as a set of locations.
This model of a range is visualized in
figure \ref{rangeSlots}.
An individual location is described by an object of the class {\tt Location}.
Each {\tt Location} element is mapped
to a particular slice of a process grid.
In general two locations are identical only
if they come from the same position in the same range.
A subscripting syntax is used to represent
location {\tt n} in range {\tt x}:
\small
\begin{verbatim}
  Location i = x [n]
\end{verbatim}
\normalsize

\begin{figure}[btp]
\centerline{\psfig{figure=rangeSlots.eps,height=0.6in}}
\caption{\label{rangeSlots}A range regarded as a set of locations,
or slots.}
\end{figure}

This is an important idea in HPJava.  By working in terms of
abstract locations---elements of distributed ranges---on can usually
respect locality of reference without resorting explicitly to low-level
local subscripts and process ids.  In fact the location can be viewed
as an abstract data type incorporating these lower-level offsets.

Publically accessible fields of {\tt Location} include {\tt dim}
and {\tt crd}.  The first is the process dimension of the parent
range.  The second is coordinate in that dimension to which the element
is mapped.  So the access to element {\tt a [17, 23]}
could now be guarded by conditionals as follows:
\small
\begin{verbatim}
  Location i = x [17], j = y [23] ;

  if(i.crd == i.dim.crd())
    if(j.crd == j.dim.crd())
      a [17, 23] = 13 ;
\end{verbatim}
\normalsize
This is still quite clumsy and error-prone.  The language provides a
second {\em distributed control construct} (analogous to {\em on}) to deal
with this common situation.  The new construct is called {\em at}, and
takes a location as its argument.  The
fragment above can be replaced with
\small
\begin{verbatim}
  Location i = x [17], j = y [23] ;

  at(i)
    at(j)
      a [17, 23] = 13 ;
\end{verbatim}
\normalsize
This is more concise, but still involves some redundancy because the
subscripts 17 and 23 appear twice.  A natural extension is to allow
locations to be used directly as array subscripts:
\small
\begin{verbatim}
  Location i = x [17], j = y [23] ;

  at(i)
    at(j)
      a [i, j] = 13 ;
\end{verbatim}
\normalsize
Locations used as array subscripts must be elements of the
corresponding ranges of the array.

The range class has a member function
\small
\begin{verbatim}
  int Range.idx(Location i)
\end{verbatim}
\normalsize
which can be used to recover the integer subscript, given a
location in the range.

There is a restriction that an {\tt at(i)} construct should only appear
at a point of execution where {\tt i.dim} is a dimension of the active
process group.  In the examples of this section this means that an {\tt
at(i)} construct, say, should normally be nested directly or indirectly
inside an {\tt on(p)} construct.

\section{Distributed loops}

As a matter of observation, good parallel algorithms don't usually
expend many lines of code assigning to isolated elements of distributed
arrays.  Sequential access to elements of parallel arrays is best
avoided.  The {\em at} mechanism of the previous section is
sometimes useful, but a more pressing need is
an idiom for {\em parallel} access to distributed array elements.

The last and most important distributed control construct in Java-Ad
is called {\em over}.
It implements a distributed parallel loop.  Conceptually it is quite
similar to the {\tt FORALL} construct of Fortran, except that the
{\em over} construct specifies exactly where its parallel iterations
are to be performed.
The argument of {\em over} is a member of the special class
{\tt Index}.  An index is associated with a particular range, which
appears in the constructor of the object.
The class {\tt Index} is a subclass of {\tt Location}, so it is syntactically
correct to use an index as an array subscript\footnote{But the effect
of such subscripting is only well-defined inside an {\em over} construct
parametrised by the index in question.}.  Here is an example of
a pair of nested {\em over} loops:
\small
\begin{verbatim}
  float [[#,#]] a = new float [[x, y]],
                b = new float [[x, y]] ;
  ...
  Index i = new Index(x), j = new Index(y) ;
  over(i)
    over(j)
      a [i, j] = 2 * b [i, j] ;
\end{verbatim}
\normalsize
The body of an {\em over} construct executes, conceptually in parallel,
for every location in the range of its index.  An individual
``iteration'' executes on just those processors holding the location
associated with the iteration.  In a particular iteration, the location
component of the index (the base class object) is equal to that location.
The net effect of the example above
should be reasonably clear.  It assigns twice the value of each element
of {\tt b} to the corresponding element of {\tt a}.  Because of the
rules about {\em where} an individual iteration iterates, the body
of an {\em over} can usually only usually combine elements of arrays
that have some simple alignment relation relative to one another.

The {\tt idx} member of range can be used in parallel updates
to give expressions that depend on global index values, as in
\small
\begin{verbatim}
  Index i = new Index(x), j = new Index(y) ;
  over(i)
    over(j)
      a [i, j] = x.idx(i) + y.idx(j) ;
\end{verbatim}
\normalsize

With the {\em over} construct we can give some more useful
examples of parallel programs.  Figure \ref{Jacobi} is the
famous Jacobi iteration for a two dimensionsional Laplace
equation.  We have used cyclic {\tt shift} to implement nearest
neighbour communications\footnote{Laplace's equation with
cyclic boundary conditions is not particularly useful,
but it illustrates the language features.
More interesting boundary conditions can easily be incorporated later.
Incidentally, this is a suitable place
to mention that the array arguments of {\tt shift} must be
{\em aligned} arrays---they must have identical distributed ranges.}.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs2 p = new Procs2(2, 2) ;

  on(p) {
    Range x = new BlockRange(100, p.dim(0)) ;
    Range y = new BlockRange(200, p.dim(1)) ;

    int [[#,#]] u = new int [[x, y]] ;

    // ... some code to initialise `u'

    int [[#,#]] unx = new int [[x, y]], upx = new int [[x, y]],
                uny = new int [[x, y]], upy = new int [[x, y]] ;

    Adlib.shift(unx, u,  1, 0, CYCL) ;
    Adlib.shift(upx, u, -1, 0, CYCL) ;
    Adlib.shift(uny, u,  1, 1, CYCL) ;
    Adlib.shift(upy, u, -1, 1, CYCL) ;

    Index i = new Index(x), j = new Index(y) ;
    over(i)
      over(j)
        u [i, j] = 0.25 * (unx [i, j] + upx [i, j] +
                           uny [i, j] + upy [i, j]) ;
  }
\end{verbatim}
\normalsize
\caption{\label{Jacobi}Jacobi iteration using {\tt shift}.}
\end{figure}

Copying whole arrays into temporaries is not an efficient way
of accessing nearest neighbours in an array.  Because this is
such a common pattern of communication, Java-Ad supports
{\em ghost regions}.  Distributed arrays can be created in such a way
that the segment stored locally is extended with some halo.  This halo
caches values stored in the segments of adjacent processes.  The cached
values are explicitly bought up to date by the library operation {\tt
writeHalo}.

An optimized version of the Jacobi program is give in figure
\ref{JacobiGhost}.  This version only involves a singe array temporary.
A new constructor for {\tt BlockRange} is provided.  This
allows the width of the ghost extensions to be specified.
The arguments of {\tt writeHalo} itself are an array with
suitable extensions and two vectors.  The first defines in each dimension
the width of the halo that must actually be updated, and the
second defines the treatment at the ends of the range---in
this case the ghost edges are updated with cyclic wraparound.
The new constructor and new {\tt writeHalo} function are simply
standard library extensions.  One new piece of syntax is needed:
the addition and subtraction operators are overloaded so that
integer offsets can be added or subtracted to {\tt Location}
objects, yielding new, shifted, locations.  The
usual access rules apply---this kind of shifted access is illegal
if it implies access to off-processor data.  It only works if
the subscripted array has suitable ghost extensions.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs2 p(2, 2) ;

  on(p) {
    Range x = new BlockRange(100, p.dim(0), 1) ;  // ghost width 1
    Range y = new BlockRange(200, p.dim(1), 1) ;  // ghost width 1

    int [[#,#]] u = new int [[x, y]] ;

    // ... some code to initialise `u'

    int [] widths = {1, 1} ;        // Widths actually updated
    Mode [] modes = {CYCL, CYCL} ;  // Wraparound at ends.

    Adlib.writeHalo(u, widths, modes) ;

    int [[#,#]] v = new int [[x, y]] ;

    Index i = new Index(x), j = new Index(y) ;
    over(i)
      over(j)
        v [i, j] = 0.25 * (u [i + 1, j] + u [i - 1, j] +
                           u [i, j + 1] + u [i, j + 1]) ;

    Adlib.copy(u, v) ;
  }
\end{verbatim}
\normalsize
\caption{\label{JacobiGhost}Jacobi iteration using {\tt writeHalo}.}
\end{figure}

\section{Subranges\label{subranges}}

A {\em subrange} is a section of a range, parametrized by a subscript
triplet.  Logically a subrange can be viewed as a subset of the
locations of the original range.  Subranges are members of the class
{\tt Range}.  Because locations in a subrange are locations of the
parent range, subranges retain an {\em alignment relation} to their
parent range.  Note that the {\em integer} subscripts for a subrange
are in the range $0, \ldots, N - 1$ where $N$ is the extent of the {\em
subrange}.  See figure \ref{subrangeSlots}.

A triplet-subscripting syntax is used for creation of subranges:
if {\tt x} is a range, then \verb$x [0 : 49]$ is a contiguous subrange
and \verb$x [1 : 98 : 2]$ is a strided subrange.

\begin{figure}[btp]
\centerline{\psfig{figure=subrangeSlots.eps,height=0.7in}}
\caption{\label{subrangeSlots}Locations of a subrange (shaded slots).}
\end{figure}

As a first application of subranges, we can uses strided subranges to
transform the Jacobi update of the previous section to a more
efficient red-black form.  The result is shown in figure \ref{RedBlack}.
The iteration is split into two phases, the first with {\tt parity} $= 0$
and the second with {\tt parity} $= 1$.  The range of the
inner {\em over} construct is either \verb$y [0 : : 2]$ or \verb$y [1 : : 2]$,
according to whether the global {\tt x} index of the outer loop
has the same or different parity (odd/even) as the current phase.
This version eliminates {\em all} temporary
arrays\footnote{Incidentally, subranges, and particularly strided
subranges, introduce an ambiguity in the definition of the shift
operators {\tt +} and {\tt -} on locations.  Is the numeric shift
measured terms of subscript relative to the subrange or of the parent
range?  As a matter of definition, the shift is always in terms of
subscript in the {\em template} range---the ultimate parent from which a
subrange is derived (by zero or more stages of triplet subscripting).}

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs2 p(2, 2) ;

  on(p) {
    Range x = new BlockRange(100, p.dim(0), 1) ;  // ghost width 1
    Range y = new BlockRange(200, p.dim(1), 1) ;  // ghost width 1

    int [[#,#]] u = new int [[x, y]] ;

    // ... some code to initialise `u'

    int [] widths = {1, 1} ;        // Widths actually updated
    Mode [] modes = {CYCL, CYCL} ;  // Wraparound at ends.

    for(int parity = 0 ; parity < 2 ; parity++) {

      Adlib.writeHalo(u, widths, modes) ;

      Index i = new Index(x) ;
      over(i) {
        Index j = new Index(y [(x.idx(i) + parity) % 2 : : 2]) ;
        over(j)
          u [i, j] = 0.25 * (u [i + 1, j] + u [i - 1, j] +
                             u [i, j + 1] + u [i, j - 1]) ;
      }
    }
  }
\end{verbatim}
\normalsize
\caption{\label{RedBlack}Red-black iteration.}
\end{figure}

As a second application involving subranges, figure \ref{Cholesky}
is a parallel version of Cholesky decomposition.  In pseudocode the
algorithm is
\begin{displaymath}
\begin{array}{l}
l_{11}=a_{11}^{1/2} \\
For\; k=1\; to\; n-1 \\
\;\;\;\; For\; s=k+1\; to\; n \\
\;\;\;\;\;\;\;\; l_{sk}=a_{sk}/l_{kk} \\
\;\;\;\;  For\; j=k+1\; to\; n \\
\;\;\;\;\;\;\;\;    For\; i=j\; to\; n \\
\;\;\;\;\;\;\;\;\;\;\;\;      a_{ij}=a_{ij}-l_{ik}l_{jk} \\
\;\;\;\; l_{k+1, k+1}= a_{k+1, k+1}^{1/2} \\
\end{array}
\end{displaymath}
The array is distributed by columns, using cyclic distribution to
improve load balancing.  The collective communication function {\tt
remap} is used to broadcast updated columns.  The {\tt remap}
function is one of the more powerful functions in the communication
library.  Like {\tt copy}, its effect is to copy data from one distributed
array to another of the same shape and type.  But {\tt copy}
(like {\tt shift}) has a restriction that its array arguments must
be {\em aligned}---{\tt copy} never introduces communication.
With {\tt remap} there is no such restriction---the mapping of the
two arrays can be unrelated.  One common application
of {\tt remap} is to {\em broadcast} data.  If the target array has no
ranges distributed over a dimension of the {\em process group} on which
it lives, {\tt remap} assumes that the result is to be stored in a
replicated fashion.  It therefore implements a broadcast.  In the
current example {\tt remap} actually implements something more
sophisticated than a simple broadcast.  In MPI terms it executes a {\em
gather-to-all} operation.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs1 p = new Procs1(P) ;
  on(p) {
    Range x = new CyclicRange(N, p.dim(0));

    float [[,#]] a = new float [[N, x]] ;

    float [[]]   b = new float [[N]] ;  // buffer

    Location l = x [0] ;
    at(l) 
      a [0, l] = sqrt(a [0, l]) ;

    for(int k = 0 ; k < N - 1 ; k++) {

      at(l)
        for(int s = k + 1 ; s < N ; s++)
          a [s, l] /= a [k, l] ;

      Adlib.remap(b [[k + 1 : ]], a [[k, k + 1 : ]]);

      Index m = new Index(x [k + 1 : ]) ;
      over(m) 
        for(int i = x.idx(m) ; i < N ; i++)
          a [i, m] -= b [i] * b [x.idx(m)] ;

      l = x [k + 1] ;
      at(l)
        a [k + 1, l] = sqrt(a [k + 1, l]) ;
    }
  }
\end{verbatim}
\normalsize
\caption{\label{Cholesky}Choleksy decomposition.}
\end{figure}

Some final comments on subranges.
Creating a triplet subscripted section 
of a distributed array implicitly creates subranges of the ranges
in the parent array.  Also,
arrays can be created directly with subranges, as in
\small
\begin{verbatim}
  Range xs = x [0 : 50] ;
  Range ys = y [1 : 198 : 2] ;

  int [[#,#]] e = new int [[xs, ys]] ;
\end{verbatim}
\normalsize
In HPF terms, {\tt e} has a non-trivial linear alignment to the
template spanned by {\tt x} and {\tt y}.
By allowing subranges (and subgroups, see section \ref{subgroups}) to
appear in array constructors we reproduce the two-level alignment model
of HPF in full generality, at little cost in terms of syntax
extensions.

\section{Subgroups\label{subgroups}}

A {\em subgroup} is some slice of a process array,
formed by restricting the process coordinates in one or more dimensions
to single values.  Process arrays (class {\tt Procs}) and subgroups
have a common base class, {\tt Group}.  In general the argument
of an {\em on} construct and the {\tt on} clause in an array constructor
is a member of {\tt Group}.  This implies that the active process group or the
group over which an array is distributed may be just some a slice of a
complete process array.

By definition, any group has a parent process array and a dimension set.
In general the dimension set is
some subset of the dimensions of the parent array.  The {\em
restriction} operation on a group takes a slice in a particular
dimension.  It is quite natural and convenient to express this
restriction procedure is in
terms of a {\em location}.  If {\tt i} is a location
in a range distributed over a dimension of {\tt p}, then
\small
\begin{verbatim}
  p / i
\end{verbatim}
\normalsize
represents a subgroup of p---the slice of {\tt p} to which location
{\tt i} is mapped.

Using the {\tt /} operator on groups explicitly is fairly unusual
practise.  But subgroups are occur naturally in two ways:
\begin{itemize}
\item
If an array {\tt a} is distributed over {\tt p}, a section of {\tt
a} will generally be
distributed over some subgroup of {\tt p}.  For example, if the only
{\em scalar} subscript in the section subscript list was location {\tt i},
the section would be distributed over the subgroup \verb$p / i$.
Triplet subscripts don't change the group---only scalar subscripts.
\item
The distributed control constructs {\em over} and {\em at} change the
active process group in a way that has not been described so far.
Assume the current active process group is {\tt p}, and {\tt i} is
a location.  Then inside the construct
\small
\begin{verbatim}
  at(i) {
    ...
  }
\end{verbatim}
\normalsize
the active process group is equal to \verb$p / i$.
If the current active process group is {\tt p}, and {\tt i} is
an {\em index}, then inside the construct
\small
\begin{verbatim}
  over(i) {
    ...
  }
\end{verbatim}
\normalsize
the active process group is equal to \verb$p / i$.  This case is
slightly more subtle, because in different parallel ``iterations''
of the loop the location component of {\tt i} has different values.  In
other words, the {\em over} construct {\em partitions} the original
active process group into several subgroups (slices) operating
independently.
\end{itemize}

To illustrate how subgroups can be used---in particular how the effect
of {\em over} on the active process group can be exploited---we
return to the matrix multiplication example of figure \ref{matmul1d}.
As a preliminary step, figure \ref{matmul1dver2} is transcribes that program
using the distributed control constructs developed over the last
few sections.  The changes are very minor.  Because the active
process group is formally changed to {\tt p}, the \verb$on p$ clauses
can be omitted from the array constructors.  Use of the {\tt crd}
inquiry to obtain the integer subscript in the {\tt x} range is
replaced by use of an {\em over} construct.
Now we want to change from distribution over a one-dimensional
process array to a two-dimensional grid.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs1 p = new Procs1(P) ;
  on(p) {
    Range x = p.dim(0) ;

    float a [[#,,]] = new float [[x, B, N]] ;
    float b [[#,,]] = new float [[x, N, B]] ;

    ... initialize `a', `b'

    float c [[#,,]] = new float [[x, B, N]] ;

    for(int s = 0 ; s < P ; s++) {

      Index ip = new Index(x) ;
      over(ip) {
        const int base = B * ((x.idx(ip) + s) % P) ;

        // c [[ip, :, base : ...]] =
        //         a [[ip, :, :]] * b' [[(ip + s) % P, :, :]] ...

        for(int ib = 0 ; ib < B ; ib++)
          for(int kb = 0 ; kb < B ; kb++) {
            float sum = 0 ;
            for(int j = 0 ; j < N ; j++)
              sum += a [ip, ib, j] * b [ip, j, kb] ;

            c [ip, ib, base + kb] = sum ;
          }
      }

      float tmp [[#,,]] = new float [[x, N, B]] ;

      Adlib.shift(tmp, b, 1, 0, CYCL) ;
      Adlib.copy(b, tmp) ;
    }
  }
\end{verbatim}
\normalsize
\caption{\label{matmul1dver2}Matrix multiplication program using
distributed control constructs.}
\end{figure}

The result is given in figure \ref{matmul2d}.  Inside the {\tt over(ip)}
construct there are {\tt P} independent active process groups corresponding
to the rows of the original grid.  The temporary array {\tt tmp} is replicated
over these rows.  The {\tt remap} operation, working independently in
the {\tt P} separate groups, implements a broadcast of the array section
representing the block of the {\tt a} stored on process \verb$(ip + s) % P$
within each group.

\begin{figure}[btp]
\small
\begin{verbatim}
  Procs2 p = new Procs2(P, P) ;
  on(p) {
    Range x = p.dim(0), y = p.dim(1) ;

    float a [[#,#,,]] = new float [[x, y, B, B]] ;
    float b [[#,#,,]] = new float [[x, y, B, B]] ;

    ... initialize `a', `b'

    float c [[#,#,,]] = new float [[x, y, B, B]] ;

    for(int s = 0 ; s < P ; s++) {

      Index ip = new Index(x), jp = new Index(y) ;
      over(ip) {

        // Broadcast a [[ip, (ip + s) % P, :, :]]...

        float [[,]] tmp = new float [[,]] ;
        Adlib.remap(tmp, a [[ip, (x.idx(ip) + s) % P, :, :]]) ;

        over(jp) {

          // c [[ip, jp, :, :]] +=
          //             a [[ip, (ip + s) % P, :, :]] *
          //                 b' [[(ip + s) % P, jp, :, :]] ...

          for(int ib = 0 ; ib < B ; ib++)
            for(int kb = 0 ; kb < B ; kb++) {
              float sum = 0 ;
              for(int jb = 0 ; jb < B ; jb++)
                sum += tmp [ib, jb] * b [ip, jp, jb, kb] ;

              c [ip, jb, ib, kb] += sum ;
            }
         }
      }

      float tmp [[#,#,,]] = new float [[x, y, B, B]] ;

      Adlib.shift(tmp, b, 1, 0, CYCL) ;
      Adlib.copy(b, tmp) ;
    }
  }
\end{verbatim}
\normalsize
\caption{\label{matmul2d}Matrix multiplication on a grid of processors.}
\end{figure}

The freedom to embed calls to collective communication functions inside
distributed control constructs is a distinctive feature of Java-Ad.
There are a few restrictions on what operations are allowed.  In
general the requirement is that all array arguments of a collective operation
should be {\em accessible} at the point of call.  An array is accesible if
it is distributed over a group contained in the active process group.

\section{Class libraries or syntax extensions?}

We have presented the Java-Ad language using a fairly liberal number
of syntax extensions.  The main extensions are:
\begin{itemize}
\item
the syntax for distributed array type signatures,
\item
the syntax for the distributed array constructors,
\item
the syntax for local subscripting of distributed arrays,
\item
the syntax for section subscripting of distributed arrays,
\item
the syntax for the three distributed control constructs,
\item
the range subscripting syntax for creating locations and subranges,
\item
the overloaded division operator for creating subgroups, and
the overloaded {\tt +} and {\tt -} operators for shifting locations.
\end{itemize}
Besides distributed arrays, we introduced four base classes
that have a special significance in the context of the above syntax
extensions:
\begin{itemize}
\item
{\tt Group}
\item
{\tt Range}
\item
{\tt Location}
\item
{\tt Index}
\end{itemize}

Some of the syntax extensions are less important than others.  The
special syntax for creating locations, subranges and subgroups could be relaced
by members of {\tt Range} and {\tt Group}:
\small
\begin{verbatim}
  Location Range.loc(int i)

  Range Range.subrng(int lb, int ub, int stride)

  Group Group.subgrp(Location i)
\end{verbatim}
\normalsize
without much inconvenience to the programmer.

It would be possible to replace the special type signatures and
constructors for
distributed arrays with a series of ordinary class types,
\begin{verbatim}
  Array1int
  Array2int
  ...
  Array1float
  Array2float
  ...
  ...
\end{verbatim}
and associated constructors.
{\tt Array1int} represents a the class of one-dimensional arrays
of {\tt int}, {\tt Array2int} represents the class of two-dimensional
arrays of {\tt int}, and so on.  The set of classes has to be quite large,
and we must fix a finite limit on the rank of an array, and the set of
array elements supported.  Unless
we allow the number of distinct array classes to grow {\em
exponentially}, it becomes impractical to distinguish between
distributed and sequential dimensions in the static type signature of
an object.  If an array is passed to a function (probably most arrays are)
it is difficult for the compiler to deduce that a particular dimension
is a sequential dimension.  The compiler may end up generating
code for global to local address conversion, even in sequential dimensions.
This is probably a serious problem for efficiency, because
sequential dimensions tend to be subscripted in an irregular way,
not in uniform {\em over} loops.  This problem could probably be
overcome by requiring the subscripting operations for sequential and
distributed dimensions explicitly differ (somehow) in the user's
program.  But that seems to complicate the language.

The syntax for the distributed control constructs and local
subscripting seems more difficult to eliminate.  In particular, the
special syntax for {\em over} is critical.  A Java-Ad translator will
do the simple but tedious code transformations needed to replace
distributed loops involving distributed arrays by local sequential
loops involving local arrays, subscripted by expressions linear in
the local loop indices.  Unless the translator does this work, data-parallel
programming becomes much less appealing.  From our point of view
this functionality is central to the SPMD data parallel model.

\end{document}

