
\documentclass{article}

\usepackage{psfig}

\begin{document}

\title{Communication in data parallel languages}

\author{Bryan Carpenter \\
NPAC at Syracuse University \\
Syracuse, NY 13244 \\
dbc@npac.syr.edu}

\maketitle

\begin{abstract}
We start by discussing the patterns of communication needed to implement
various constructs in High Performance Fortran.  Then several communication
libraries that have been developed to support these kinds of
communications will be reviewed.  Finally, as a case study a, scheme
for implementing the array remapping operation used in general array
assignments is described in detail.
\end{abstract}

\newpage

\tableofcontents

\newpage

\section{Patterns of communication}

\label{sec:patterns}

In this section we will discuss some patterns of communication that
arise in when a language like HPF is translated to a SPMD program.
First we will cover the regular communication patterns that
arise from simple array assignments, nearest-neighbour problems,
and array intrinsics.  Then the irregular communication patterns
implied by data-parallel statements with non-linear subscript expressions
will be introduced.  Finally we cover some essentially task-parallel
codes that can be translated most effectively in terms of ``one-sided
communication''.

\subsection{Array assignments}

\label{sec:arrayassignments}

In HPF, variants of the innocent-looking assignment,
\small
\begin{verbatim}
      A = B
\end{verbatim}
\normalsize
hide a variety of interesting communication patterns.  The terms {\tt
A} and {\tt B} may be any conforming arrays.  Let's run through some
examples.

As a first example consider:
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

      REAL A(50), B(50)
!HPF$ DISTRIBUTE A(BLOCK) ONTO P
!HPF$ DISTRIBUTE B(BLOCK) ONTO P

      A (1:49) = B (2:50)
\end{verbatim}
\normalsize
The assignment variable and expression are sections of arrays with the
same distribution.  But the alignments of their elements are shifted by
one place relative to each another.  In translating this assignment,
some communication between neighbouring processes will be needed.  
The situation is illustrated in Figure \ref{fig:assignshift}.
In a translation to a SPMD program, the communications might be
implemented using the point-to-point primitives {\tt MPI\_SEND} and
{\tt MPI\_RECV}, or in a single call using {\tt MPI\_SENDRECV}.

\begin{figure}[btp]
\centerline{\psfig{figure=assignshift.eps,width=4.8in}}
\caption{Assignment involving a shift in alignment.  Heavy arrows represent
interprocessor communication.}
\label{fig:assignshift}
\end{figure}

The ``shift'' pattern of communication is actually quite common in
practical situations, and the next subsection will discuss it in
more detail.  But it is slightly contrived as an example of an array
assignment.  Usually the right- and left-hand-side arrays will not
have such a simple alignment relationship.
Here is another example involving array sections\footnote{As a matter
of fact, any communication pattern in HPF involving array sections can be
reproduced using whole arrays.  One just has to give those arrays suitable
alignments to templates.  We choose to work with sections because they
seem more concrete.}:
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

      REAL A(50,50), B(50)
!HPF$ DISTRIBUTE A(*, BLOCK) ONTO P
!HPF$ DISTRIBUTE B(BLOCK) ONTO P

      A (:, 1) = B
\end{verbatim}
\normalsize
In this case the first dimension of {\tt A} is collapsed, so the target
of the assignment is a section that lives entirely on one processor.
The communication pattern is visualized in Figure \ref{fig:assigngather}
You may recognize this as being essentially the communication pattern
implemented by the MPI collective operation {\tt MPI\_GATHER} (or {\tt
MPI\_GATHERV}).  If the direction of the assignment had been reversed:
\small
\begin{verbatim}
      B = A (:, 1)
\end{verbatim}
\normalsize
then the communication pattern would correspond to the MPI
operation {\tt MPI\_SCATTER} or {\tt MPI\_SCATTERV}.
If the assignment target had a {\em replicated}, collapsed alignment:
\small
\begin{verbatim}
      REAL A(50,50), B(50), C(50)
!HPF$ DISTRIBUTE A(*, BLOCK) ONTO P
!HPF$ DISTRIBUTE B(BLOCK) ONTO P
!HPF$ ALIGN C(I) WITH A(I, *)

      C = B
\end{verbatim}
\normalsize
the communication pattern would be to broadcast all elements of the
distributed array {\tt B} to all processors---we recognize this as being
like the MPI collective {\tt MPI\_ALLGATHER} or {\tt MPI\_ALLGATHERV}.

\begin{figure}[tbp]
\centerline{\psfig{figure=assigngather.eps,width=4.8in}}
\caption{Assignment involving gathering array data to a single processor.}
\label{fig:assigngather}
\end{figure}

Finally we can easily write assignments that behave like {\tt MPI\_ALLTOALL}
or {\tt MPI\_ALLTOALLV}.  As an exercise the student may wish to verify
that the following example qualifies:
\small
\begin{verbatim}
      REAL A(50,50), B(50,50)
!HPF$ DISTRIBUTE A(*, BLOCK) ONTO P
!HPF$ DISTRIBUTE B(BLOCK, *) ONTO P

      A = B
\end{verbatim}
\normalsize
as does:
\small
\begin{verbatim}
      REAL A(50), B(50)
!HPF$ DISTRIBUTE A(BLOCK) ONTO P
!HPF$ DISTRIBUTE B(CYCLIC) ONTO P

      A = B
\end{verbatim}
\normalsize

On the one hand these examples illustrate the power of the HPF language.
The effect of a complicated collective call in MPI can often be
expressed very concisely in HPF as a simple array assignment.  On the
other hand, we get an idea of how complicated the communication patterns
implied by simple-looking HPF statements can be.

\subsection{Stencil problems and ghost areas}

\label{sec:stencil}

Stencil updates, in which each element of an array is updated in terms
of some fixed footprint or ``stencil'' of neighbouring elements, have
been an important target for parallel computing.
They are quite common in practise, arising for example in solution of
partial differential equations, simulation of cellular automata, and
image processing applications.  This situation is certainly amenable to
massive parallelism, but it does require some communication.

These problems {\em can} be reduced to the kinds of array assignment
discussed in Section \ref{sec:arrayassignments}.
The well-known example:
\small
\begin{verbatim}
      FORALL (I = 2:N-1, J = 2:N-1)
     &    U(I,J) = 0.25 * (U(I,J-1) + U(I,J+1) + U(I-1,J) + U(I+1,J))
\end{verbatim}
\normalsize
can be recast as:
\small
\begin{verbatim}
     U(2:N-1, 2:N-1) = 0.25 * (U(2:N-1, 1:N-2) + U(2:N-1, 3:N) +
    &                          U(1:N-2, 2:N-1) + U(3:N, 2:N-1))
\end{verbatim}
\normalsize
The compiler could translate the array assignment code by
introducing a series of temporary arrays {\tt T1}, \ldots, {\tt T1}
all aligned with the section \verb$U(2:N-1, 2:N-1)$:
\small
\begin{verbatim}
     T1 = U(2:N-1, 1:N-2) 
     T2 = U(2:N-1, 3:N)
     T3 = U(1:N-2, 2:N-1) 
     T4 = U(3:N, 2:N-1))
     U(2:N-1, 2:N-1) = 0.25 * (T1 + T2 + T3 + T4)
\end{verbatim}
\normalsize
Due to the alignment of the temporary arrays, the last statement is purely
local computation; the first four assignments absorb all communication.
If the the arrays here have been distributed in the BLOCK style---usually
optimum for this kind of problem---these communications follow a pattern
similar to Figure \ref{fig:assignshift}.

In terms of the volume of interprocessor communication this is
a reasonable approach.  But it has two big drawbacks.
First it uses a lot of temporary storage---four whole arrays.  Secondly,
although the first four assignments involve no computation and only
a modest amount of communication at the edges, they introduce a lot
of extra memory-to-memory copying.

A sequential version of the original loop has good locality properties
in terms of its use of memory and cache.  Because the elements on the
left-hand-side and the operands in the right-hand-side of each assignment
are typically clustered together in memory (spatial locality), and
because elements are often reused in consecutive iterations, or iterations
not far apart in iteration space (temporal locality), a
straightforward transcription of the algorithm is likely to make good use
of the cache of a modern microprocessor \cite{Wolfe96}.
But if the program is split into consecutive array assignments as indicated
above, much of this locality will be lost.  Generally whole arrays
will not fit in cache memory.  Cache will be loaded from main memory and
flushed back for each of the five loops associated with the five array
assignments.

In fact this is a generic problem for the ``array syntax'' style of
programming encouraged by Fortran 90.  As we have seen in preceding
lectures, this style evolved in large part to support a
generation of SIMD and vector processors.  On contemporary microprocessors,
memory access costs usually dominate over the costs of arithmetic.
Ironically a compiler may end up working quite hard to fuse the array
assignments of a Fortran 90 program into sequential loops that have
good memory locality.

A better approach than copying whole arrays of neighbours is to allocate
the local segments of the distributed array {\tt U} with a small amount
of extra space
around the edges---so-called {\em ghost regions}.
The translation of the FORALL statement above could be something like:
\small
\begin{verbatim}
      REAL U(0 : BLK_SIZE1 + 1, 0 : BLK_SIZE2 + 1)
      REAL T(BLK_SIZE1, BLK_SIZE2)

      ... Update ghost regions of U with values from neighbours

      DO L1 = 1, BLK_COUNT1
          DO L2 = 1, BLK_COUNT2
              T(L1, L2) = 0.25 * (U(L1, L2 - 1) + U(L1, L2 + 1) 
     &                          + U(L1 - 1, L2) + U(L1 + 1, L2))
          ENDDO
      ENDDO

      ... Copy T to interior part of U (local copying)
\end{verbatim}
\normalsize
The communications required to update the ghost regions are
illustrated in Figure \ref{fig:updateghost}.
For the sake of definiteness we will refer to the interior part
of the local segment---the non-ghost part---as the {\em physical} segment.
The communication pattern illustrated here updates the cells in the ghost
regions with the values from the corresponding cells in the physical
segments.

\begin{figure}[tbp]
\centerline{\psfig{figure=updateghost.eps,width=4in}}
\caption{Communications needed to update ghost regions.  Only
communications involving the top left-hand processor are shown.}
\label{fig:updateghost}
\end{figure}

\subsection{Reductions and other transformational intrinsics}

Fortran implementations come with a set of built-in functions and
subroutines called the {\em intrinsics}.
The {\em transformational intrinsics} are a subset of the intrinsics that
peform functions on whole arrays, sometimes returning whole array results.
In HPF these functions can operate on distributed arrays, implying
new patterns of communication.

Some of the transformational intrinsics just reshuffle elements of
arrays in regular ways.  For example if the {\tt CSHIFT}, {\tt EOSHIFT},
{\tt TRANSPOSE} and {\tt SPREAD} intrinsics are applied to distributed
arrays, they imply patterns of data remapping that are quite similar
to patterns in the array assignments of section \ref{sec:arrayassignments}.
The reduction intrinsics, including {\tt SUM}, {\tt PRODUCT}, {\tt
MAXVAL}, {\tt MINVAL}, {\tt ALL}, {\tt ANY} and several others, introduce
qualitatively different communication patterns.  They reduce an array
either to a scalar or a lower-rank array by arithmetically combining
elements along some or all dimensions.  Typically the associated
communications occur on some spanning tree.  In MPI terms they are
related to, and might be implemented in terms of, {\tt MPI\_REDUCE}.

{\em Parallel prefix} operations can also be expressed.  The HPF standard
includes them as standard library functions (rather than intrinsics).
They can also be expressed in terms of FORALL statements:
\small
\begin{verbatim}
      FORALL (I = 1 : N)  RES (I) = SUM(A (1 : I))
\end{verbatim}
\normalsize
The $i$th element of the result contains the sum of all elements in
{\tt A} up to and including its $i$th element.
The pattern of
communication needed to implement the prefix algorithm {\em efficiently} is
different to the one used in global reduction operations\footnote{The
CM Fortran compiler, for example, could recognize the kind of FORALL
statement above as a parallel prefix and translate it appropriately.}.
In MPI terms parallel prefixes are related to, and might be implemented
in terms of, {\tt MPI\_SCAN}.

\subsection{General subscripting in array parallel statements}

\label{sec:indirectsubscripts}

Quite often a parallel program needs to execute some code like
\small
\begin{verbatim}
      FORALL (I = 1 : 50)  RES (I) = A (IND (I))
\end{verbatim}
\normalsize
where {\tt A} is one distributed array, and {\tt IND} is some
other distributed array.  The special feature is that some subscript
expressions (the subscripts for {\tt A} in this case) are not a simple
linear functions of the index variables.
This operation cannot be reduced to the kind of a simple array assignment
described in Section \ref{sec:arrayassignments}.  Even in the most
favourable case, when the arrays {\tt RES} and {\tt IND} are identically
aligned, the data movement cannot be reduced to a single MPI collective
operation.

Assuming {\tt RES} and {\tt IND} {\em are} aligned, the owner of the element
{\tt RES(I)} also owns the subscript {\tt IND(I)}.  So the processor
that holds the target variable of the assignment can compute where the
required element of {\tt A} lives.  This is generally on some other
processor.  Unfortunately that other processor---the one that owns the
{\tt A} element---does {\em not} have the information to tell the
value is required by the first processor.
It seems that there must be at least two phases of communication---one
in which the target processors send out requests to the owners of the
required data, and one in which the owners return the data.
The situation is illustrated in Figure \ref{fig:assignind} for the
case
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

      REAL A(50)
!HPF$ DISTRIBUTE A(CYCLIC) ONTO P

      REAL RES(50)
      INTEGER IND(50)
!HPF$ DISTRIBUTE RES(BLOCK) ONTO P
!HPF$ DISTRIBUTE IND(BLOCK) ONTO P

      IND = (/ 5, 41, 7, ... /)
\end{verbatim}
\normalsize
The first processor can deal with the assignment of the values
\verb$A(5)$ and \verb$A(41)$ to the first two elements of {\tt RES} locally.
But the element \verb$A(7)$ lives on the third processor, and initially
the third processor has no way to know that this element is needed by
the first.  Some two-way communication is needed.

\begin{figure}[tbp]
\centerline{\psfig{figure=assignind.eps,width=4.8in}}
\caption{Assignment involving an indirect reference.}
\label{fig:assignind}
\end{figure}

There are many generalizations of this problem.  The non-linear subscript
expressions may appear on the left-hand-side rather than the right-hand-side:
\small
\begin{verbatim}
      FORALL (I = 1 : 50)  RES (IND (I)) = A (I)
\end{verbatim}
\normalsize
(More generally, non-linear subscripts may appear on {\em both} sides.)
The non-linear subscript
may be a function or expression rather than an indirection array.
The irregularly subscripted arrays may be multidimensional:
\small
\begin{verbatim}
      FORALL (I = 1 : 50)  RES (I) = A (IND1 (I), IND2(I))
\end{verbatim}
\normalsize
As a special case it may be that the subscripts are actually be linear, but
because the subscripts in orthogonal dimensions are not independent,
the assignment can still not be reduced to a simple array assignment:
\small
\begin{verbatim}
      FORALL (I = 1 : 50)  RES (I) = A (I, I)
\end{verbatim}
\normalsize
The loop itself may be multidimensional:
\small
\begin{verbatim}
      FORALL (I = 1 : 50, J = 1 : 50)  RES (I, J) = A (IND (I, J))
\end{verbatim}
\normalsize
These examples are special cases of two (rather complicated) general families:
\small
\begin{tabbing}
\verb$      FORALL ($$i_1$\verb$ = 1 : $$n_1$\verb$, ..., $$i_R$\verb$ = 1 : $$n_R$\verb$)$ \\
\verb$     &     RES ($$i_1, \ldots, i_R$\verb$) = SRC ($$\mbox{\tt IND}_1$\verb$ ($$i_1, \ldots, i_{R}$\verb$), ...,$
$\mbox{\tt IND}_{S}$\verb$ ($$i_1, \ldots, i_R${\tt ))}
\end{tabbing}
\normalsize
where the {\tt IND} arrays are aligned with the {\tt RES} array,
and
\small
\begin{tabbing}
\verb$      FORALL ($$i_1$\verb$ = 1 : $$n_1$\verb$, ..., $$i_S$\verb$ = 1 : $$n_S$\verb$)$ \\
\verb$     &     RES ($$\mbox{\tt IND}_1$\verb$ ($$i_1, \ldots, i_{S}$\verb$), ..., $$\mbox{\tt IND}_{R}$\verb$ ($$i_1, \ldots, i_S$
\verb$)) = SRC ($$i_1, \ldots, i_S$\verb$)$
\end{tabbing}
\normalsize
where the {\tt IND} arrays are aligned with the {\tt SRC} array.
We will call these {\em generalized gather} and {\em generalized scatter}
operations respectively\footnote{Note that these operations are {\em not}
related to the must simple {\tt MPI\_GATHER} and {\tt MPI\_SCATTER}
operations---MPI's use of the terms ``gather'' and ``scatter'' is
slightly unconventional.}.
A large class of FORALL statements can be reduced to a series of
generalized gather and scatter operations, interspersed with some
ordinary array assignments (although there is no guarantee that this is
always the most {\em efficient} way to translate a FORALL statement).

\subsection{Accessing remote data in task-parallel code}

All the previous examples considered patterns of communication
occurring in array parallel statements---array assignments
or FORALL statements.  These communication patterns are quite naturally
treated as generalized collective operations.  But there are situations
in HPF---and in general SPMD programming---where this approach is
not readily applicable.

One example is the INDEPENDENT DO loop of HPF, which takes the form:
\small
\begin{verbatim}
!HPF$ INDEPENDENT
      DO i = 1, 10
          ...
      END DO
\end{verbatim}
\normalsize
The INDEPENDENT directive asserts that there are no data dependences
between individual iterations of the following loop, and the iterations
may therefore be executed in parallel\footnote{Bear in mind that this
independence at the logical variable level implies nothing about the
home processors of accessed variables, and thus nothing about whether
interprocessor communications are needed to translate statements.}.
Unlike the FORALL statement, which explicitly limits the code executed in
parallel to simple assignments, the body of an INDEPENDENT DO can involve
any Fortran construct, including conditionals, loops and procedure calls.
So the patterns of access to remote data inside parallel ``iterations''
may vary in unpredictable ways from one iteration to the next.
It may become difficult to do any advance orchestration of data exchanges.
An HPF the compiler is free to ignore the INDEPENDENT directive
if it decides the loop is too complex to parallelize.  But this may deprive
the programmer of one of the few options in HPF for expressing the
task-farming style of parallelism.

Actually there is at least one other way to express task parallelism in HPF.
A user-defined procedure with the PURE attribute can be called from
within a FORALL statement:
\small
\begin{verbatim}
      PURE REAL FUNCTION FOO(INTEGER I)
      ...
      END

      ...
      FORALL (I = 1 : N)  RES (I) = FOO(I)
\end{verbatim}
\normalsize
There are quite strict restrictions on PURE procedures, but nothing to
prevent a procedure from reading elements of global distributed
data---a distributed array in a COMMON block, for example.
Unfortunately this makes it difficult or impossible for the compiler
to determine at the point of call of FOO exactly what
remote variables it will access.  For example, the actual behaviour of the
program might be similar to the first example of Section
\ref{sec:indirectsubscripts}:
\small
\begin{verbatim}
      PURE REAL FUNCTION FOO(INTEGER I)

      REAL RES(50)
      INTEGER IND(50)
!HPF$ DISTRIBUTE RES(BLOCK) ONTO P
!HPF$ DISTRIBUTE IND(BLOCK) ONTO P
      COMMON /GLOBALS/ RES, IND

      RETURN RES (IND (I))
      END
\end{verbatim}
\normalsize
But by the time \verb$RES(IND(I))$ is accessed, instances of the function
{\tt FOO} have already been dispatched to execute independently
across the available set of processors.  In a real sense, once inside FOO
processors are no longer sharing a single ``loosely synchronous'' thread
of control.  It is difficult to see how the parallel invocations
of FOO can behave collectively.  In particular if the underlying
model is MPI point-to-point communication it is difficult
to see how the owner of a particular array element can always be ready
to send an element when its value is accessed by a peer processor.

INDEPENDENT DO loops have similar problems, compounded because they do
not have the restrictions on PURE procedures that prevent them from
{\em writing} to global variables.  If this sort of code is to be compiled
to run in parallel the most practical approach is probably to assume
the availability of {\em one-sided communication}.  The MPI 2 standard
added this functionality to MPI, but it is still not widely implemented.

\section{Libraries for distributed array communication}

As we have seen, the communication patterns implied by languages
like HPF can be complex.  Often it is impractical for a compiler to
generate all the low-level message passing instructions needed to
execute these communications.  Instead the compiler may choose to emit
code with calls to higher-level libraries for manipulating distributed
array data.  By analogy with the run-time support libraries used for
memory management (and so on) in sequential languages, the data parallel
compiler's library for handling distributed arrays is often called
its {\em run-time library}.

In this section we will discuss some libraries that have been used,
or could be used, in this role.

\subsection{The PARTI primitives}

The CHAOS/PARTI series of libraries was developed at the University
of Maryland.

The original PARTI library was designed to deal with irregular scientific
computations.  
A classic example is a physical problem discretized on an unstructured
mesh.  A characteristic inner loop for this kind of problem might look
something like:
\small
\begin{verbatim}
      DO I = 1, NEDGE
          Y(EDGE1(I)) = Y(EDGE1(I)) + F(X(EDGE1(I)), X(EDGE2(I)))
          Y(EDGE2(I)) = Y(EDGE2(I)) + G(X(EDGE1(I)), X(EDGE2(I)))
      ENDDO
\end{verbatim}
\normalsize
We assume the problem is defined on some graph (or mesh).  The arrays
\verb$X$ and \verb$Y$ are defined over the nodes of the graph.
The {\tt I}th edge in the graph connects the nodes \verb$EDGE1(I)$ and
\verb$EDGE2(I)$.  The value of \verb$Y$ at a node is a sum of terms, each
depending on the value of \verb$X$ at the ends of an edge connected to
the node.

PARTI is particularly optimized for problems that have some ``locality
of reference''.  It assumes that a mapping of nodes to processors
can be chosen such that {\em most} edges connect pairs of nodes
on the same processor.  The situation is illustrated in Figure
\ref{fig:graphpartition}.  
The local loops will be over the
edges held on each processor.  Here only one the edge (2, 3) held on
processor P0 causes non-local references.

\begin{figure}[tbp]
\centerline{\psfig{figure=graphpartition.eps,width=4.8in}}
\caption{An irregular problem graph, exhibiting locality of reference.}
\label{fig:graphpartition}
\end{figure}

Figure \ref{fig:graphpartition} is copied from one of the
illustrations in \cite{Ponnusamy95}.  It assumes an {\em irregular}
distribution of {\tt X} and {\tt Y}: elements 1, 2, 5 of the data arrays
are stored on the first processor and elements 3, 4, 6 are stored on the
second processor.  This property is not particularly important to
the discussion here.  We could assume that the node numbering is
permuted so that arrays like {\tt X} and {\tt Y} have a block distribution.

In any case the important point is that there is a class of problems
with the following property: edges and nodes of the problem graph can
be partitioned in such a way that most references become local.  To be
specific, they can be partitioned so that the {\em majority} of locally
held elements of indirection vectors like {\tt EDGE1}, {\tt EDGE2}
reference locally held elements of data arrays, like {\tt X} and {\tt Y}.

Based on this observation, PARTI assumes irregular loops are parallelized
by a technique similar to the method of ghost regions discussed for
regular stencil problems in section \ref{sec:stencil}.  First the
indirection vectors are preprocessed to convert global index values to
{\em local subscripts}.  The locality property of the partition implies
that the majority of these local subscripts refer to locally held data
elements.  If the global index actually referenced a data element held
on another processor, the local subscript references an element in a
``ghost extension'' of the local data segment.  These ghost regions
are filled or flushed by suitable PARTI primitives: collective
communication routines, called outside the local processing loop.

Here is a simpler sequential loop with irregular accesses:
\small
\begin{verbatim}
      DO I = 1, N
          X(IA(I)) = X(IA(I)) + Y(IB(I))
      ENDDO
\end{verbatim}
\normalsize
A parallel version (from \cite{Das94}) is given in Figure
\ref{fig:simpleirreg}.

\begin{figure}[tbp]
\small
\begin{verbatim}
C Create required schedules (Inspector):

      CALL LOCALIZE(DAD_X, SCHEDULE_IA, IA, LOCAL_IA, I_BLK_COUNT,
                    OFF_PROC_X)

      CALL LOCALIZE(DAD_Y, SCHEDULE_IB, IB, LOCAL_IB, I_BLK_COUNT,
                    OFF_PROC_Y)

C Actual computation (Executor):

      CALL GATHER(Y(Y_BLK_SIZE + 1), Y, SCHEDULE_IB)

      CALL ZERO_OUT_BUFFER(X(X_BLK_SIZE + 1), OFF_PROC_X)

      DO L = 1, I_BLK_COUNT
          X(LOCAL_IA(I)) = X(LOCAL_IA(I)) + Y(LOCAL_IA(I))
      ENDDO

      CALL SCATTER_ADD(X(X_BLK_SIZE + 1), X, SCHEDULE_IA)
\end{verbatim}
\normalsize
\caption{PARTI code for simple irregular loop.}
\label{fig:simpleirreg}
\end{figure}

The first call to the subroutine {\tt LOCALIZE} deals with the
\verb$X(IA(I))$ terms.  It does two things.  It translates the {\tt
I\_BLK\_COUNT} global subscripts in {\tt IA} to local subscripts, returned
in the array {\tt LOCAL\_IA}, and it sets up a {\em communication
schedule}.  A handle to this data structure is returned in {\tt
SCHEDULE\_IA}.

A communication schedule is created by analysing the requested set of
accesses, sending lists of accessed elements to the processors
that own them where necessary, detecting appropriate aggregations and
redundancy eliminations, and so on.  The end result is some digested list
of messages that must be sent and received, including the local sources and
destinations of the data in those messages.

Another input parameter to {\tt LOCALIZE} is the distribution of the data
array---{\tt DAD\_X} in the first call.  Another output parameter is the
number of elements in the ghost region that will actually be needed---{\tt
OFF\_PROC\_X} here.

The second call to {\tt LOCALIZE} performs a similar analysis for the
term \verb$Y(IB(I))$.  

Together these calls comprise what is called the {\em inspector phase}
for the loop.  It is followed by the {\em executor phase}, in which
results are actually computed and data is actually communicated.

The collective call to {\tt GATHER} communicates necessary element
values from physical segments of {\tt Y}---the second argument---into
the target ghost regions for {\tt Y}, which starts at
\verb$Y(Y_BLK_SIZE + 1)$---the first argument.  The third argument
is the communication schedule for this operation.  The call to {\tt
ZERO\_OUT\_BUFFER} just sets all elements in the ghost region of {\tt
X} to zero.

The main loop is self-explanatory.  Contributions to locally owned
\verb$X(IA)$ elements are accumulated directly into the local physical
segment of {\tt X}.  Contributions to non-locally owned elements are
accumulated into the ghost region of {\tt X}.

Finally the call {\tt SCATTER\_ADD} sends the values in the ghost
region of {\tt X} to the relevant owners, where they are added to
the appropriate elements in the physical region of the array segment.
This is an example of a {\em combining scatter} operation.

Besides the PARTI software, a major contribution here is the elaboration
of the {\em inspector-executor model}, and the insight that construction
of communication schedules should be separated from execution of those
schedules.  One immediate benefit of this separation arises in the common
situation where the form of the inner loop (the pattern of subscripting)
is constant over many iterations of some outer loop.  The same communication
schedule can be reused many times---in other words the inspector phase
can be lifted out of the main loop.

The importance of the inspector-executor model and the idea of
communication schedules are not tied to the details of the
PARTI primitives.  For example they are not dependent on the particular
assumptions about locality, or the special use of ghost regions, or even
particularly specific to {\em irregular} computations.

\subsection{Multiblock PARTI}

\subsection{The Global Array Toolkit}

\section{Adlib}

\subsection{Background and functionality}

\subsection{Case study: implementation of the {\tt Remap} schedule}

\bibliographystyle{plain}
\bibliography{runtime}

\end{document}

