
\documentclass{article}

\usepackage{psfig}

\begin{document}

\title{Issues in Translation of High Performance Fortran}

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

\maketitle

\begin{abstract}
This lecture reviews some issues in the translation of an HPF
program to a SPMD program.  The need for a run-time array descriptor
is discussed.  An abstract model of a suitable Distributed Array
Descriptor is developed.  Illustrative translations of HPF fragments
to Fortran and C++ are presented.
\end{abstract}

\newpage

\tableofcontents

\newpage

\section{Introduction}

\label{sec:introduction}

Here is an exceedingly simple HPF program:
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

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

      FORALL (I = 1:50)  A(I) = 1.0 * I
\end{verbatim}
\normalsize
We want to convert this data-parallel program into a program that can be
executed directly on a typical MIMD computer or cluster of workstations.
We assume that the target machine provides a Fortran compiler and an MPI
library.  A reasonable translation is given in Figure \ref{fig:simpleMPI}.
This looks a little complicated.  We will break it down into its parts.

\begin{figure}[tbp]
\small
\begin{verbatim}
      INTEGER W_RANK, W_SIZE, ERRCODE

      INTEGER BLK_SIZE
      PARAMETER (BLK_SIZE = (50 + 3) / 4)

      REAL A(BLK_SIZE)

      INTEGER BLK_START, BLK_COUNT
      INTEGER L, I

      MPI_INIT()

      MPI_COMM_RANK(MPI_COMM_WORLD, W_RANK, ERRCODE)
      MPI_COMM_SIZE(MPI_COMM_WORLD, W_SIZE, ERRCODE)

      IF (W_RANK < 4) THEN
          BLK_START = W_RANK * BLK_SIZE
          IF(50 - BLK_START >= BLK_SIZE) THEN
              BLK_COUNT = BLK_SIZE
          ELSEIF(50 - BLK_START > 0) THEN
              BLK_COUNT = 50 - BLK_START
          ELSE
              BLK_COUNT = 0
          ENDIF

          DO L = 1, BLK_COUNT
              I = BLK_START + L
              A(L) = 1.0 * I
          ENDDO
      ENDIF

      MPI_FINALIZE()
\end{verbatim}
\normalsize
\caption{MPI translation of a simple HPF program.}
\label{fig:simpleMPI}
\end{figure}

First we have a series of statements that simply set up the environment:
\small
\begin{verbatim}
      INTEGER W_RANK, W_SIZE, ERRCODE
      ...

      MPI_INIT()

      MPI_COMM_RANK(MPI_COMM_WORLD, W_RANK, ERRCODE)
      MPI_COMM_SIZE(MPI_COMM_WORLD, W_SIZE, ERRCODE)
      ...

      MPI_FINALIZE()
\end{verbatim}
\normalsize
So far as the present example is concerned, their only important role is
to define the variables {\tt W\_RANK}, {\tt W\_SIZE}.

Secondly we have statements associated with the declaration of the
distributed array:
\small
\begin{verbatim}
      INTEGER BLK_SIZE
      PARAMETER (BLK_SIZE = (50 + 3) / 4)

      REAL A(BLK_SIZE)
\end{verbatim}
\normalsize
These compute a bound on the number of elements of the
distributed array held on any processor, and allocate a local array
large enough the hold those elements.  In the present case the size
of the distributed array (50) and the number of processors over
which it is distributed are both known at compile-time.  So the block
size, $\lceil 50 / 4 \rceil$, can be computed at compile-time, and
the array segment can be declared statically.

The remainder of the code is associated with transation of the
FORALL statement assigning {\tt A}.  There is a test to see if the
local processor is a member of the processor arrangement {\tt P}
over which the elements of the assigment variable {\tt A} are
distributed:
\small
\begin{verbatim}
      IF (W_RANK < 4) THEN
          ...
      ENDIF
\end{verbatim}
\normalsize
We will always adopt the policy that an HPF processor arrangement
may have any number of abstract processors up to the number of
physical processors available at run-time (or conversely number of
physical processors provided at run-time must be at least the size
size of the largest processor arrangement declared in the program).
This is not the only possible policy, but it is a reasonable one,
and it is consistent with what the HPF standard says\footnote{To be
thorough we ought to associate a run-time check that {\tt W\_SIZE >= 4}
with the declaration of {\tt P}.}.

If the local processor {\em is} in {\tt P}, we compute some parameters of
the locally held subset of distributed array elements---the values {\tt
BLK\_START} and {\tt BLK\_COUNT}.  These characterize the position of
the local segment in the global index space, and the number of elements
in that segment:
\small
\begin{verbatim}
      INTEGER BLK_START, BLK_COUNT
      ...

          BLK_START = W_RANK * BLK_SIZE
          IF(50 - BLK_START >= BLK_SIZE) THEN
              BLK_COUNT = BLK_SIZE
          ELSEIF(50 - BLK_START > 0) THEN
              BLK_COUNT = 50 - BLK_START
          ELSE
              BLK_COUNT = 0
          ENDIF
\end{verbatim}
\normalsize
Finally we have the actual loop over the local elements:
\small
\begin{verbatim}
      INTEGER L, I
      ...
          DO L = 1, BLK_COUNT
              I = BLK_START + L
              A(L) = 1.0 * I
          ENDDO
\end{verbatim}
\normalsize
The global index is computed as a linear function of the local loop
index {\tt L}, and the assignment to the array element is performed.

Of course there was quite a lot of book-keeping in the MPI translation,
but nothing conceptually difficult.

Now consider this superficially similar fragment of HPF:
\small
\begin{verbatim}
      SUBROUTINE INIT(A)

      REAL A(50)
!HPF$ INHERIT A

      FORALL (I = 1:50)  A(I) = 1.0 * I
      END
\end{verbatim}
\normalsize
The INHERIT directive specifies that the mapping of the dummy
argument {\tt A} should be the same as the actual argument,
whatever that is.  At compile-time we have no information about
that mapping, other than it is a legal HPF decomposition.  If
the main program was
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

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

      CALL INIT(B)
\end{verbatim}
\normalsize
then the translated subroutine ought to behave exactly as in the
previous example.  But at the time {\tt INIT} is translated this
fact is generally not known.  The subroutine may be part of a library
compiled long before the calling program is written.

\begin{figure}[tbp]
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)

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

      CALL INIT(B)
\end{verbatim}
\normalsize
\psfig{figure=callblock.eps,width=4.8in}
\small
\begin{verbatim}

!HPF$ PROCESSORS P(4)

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

      CALL INIT(B)
\end{verbatim}
\normalsize
\psfig{figure=callcyclic.eps,width=4.9in}
\small
\begin{verbatim}

!HPF$ PROCESSORS P(4)

      REAL B(100)
!HPF$ DISTRIBUTE B(BLOCK) ONTO P

      CALL INIT(B(1:100:2))
\end{verbatim}
\normalsize
\psfig{figure=callstride.eps,width=4.9in}
\small
\begin{verbatim}

!HPF$ PROCESSORS Q(2, 2)

      REAL B(6, 50)
!HPF$ DISTRIBUTE B(BLOCK, BLOCK) ONTO Q

      CALL INIT(B(2, :))
\end{verbatim}
\normalsize
\centerline{\psfig{figure=callrow.eps,width=3.5in}}
\caption{Mapping of {\tt A} for some possible calls to {\tt INIT}.}
\label{fig:callexamples}
\end{figure}

A few of the possibilities for mapping of the actual argument are
illustrated in Figure \ref{fig:callexamples}.  The first case has
already been discussed.  The array {\tt A} will be
divided into 4 contiguous blocks of sizes (13,13,13,11).  In the
second case the blocking is different---(13,13,12,12)---and the
formula for computing the index value {\tt I} inside the local loop
will be quite different.  The third case illustrates that {\tt A}
might actually be some section of an array.  In this example
the blocking is (13,12,13,12) and the translation must somehow
take account that the subscripting into the local segment of
the array is strided, and there is an offset which is sometimes
zero and sometimes one.  In the final case {\tt A} is again a section,
but now the parent array is two-dimensional, and we have to deal with
the fact that {\tt A} as a whole is mapped only to a portion of the
processor arrangement.  The physical processors associated with
the bottom row of {\tt Q} should do nothing.

Somehow the procedure {\tt INIT} needs to be translated in such
a way that it can deal with all these cases, and many others.

We have deliberately started with a very simple case.  The parallel part
of the program only deals with a single one-dimensional array, and it
is clear that the translation requires absolutely no inter-processor
communication.  Nevertheless we have already stumbled into a whole
wad of complexity associated with the translation of procedures.

This is not an artificial example.  General purpose libraries for dealing
with distributed arrays are likely to encounter just this kind of problem.
As a specific example, the array communication libraries discussed
in the next lecture have exactly this requirement: they must deal with
distributed arrays that have unrestricted HPF mappings, unknown at the
time the library code is compiled\footnote{The functions in such a library
need not actually be {\em written} in HPF.  Regardless, the author of the
library must deal with the same problems as an HPF compiler translating
procedures with INHERIT mappings.}.

This lecture mainly discusses a run-time parameterization of HPF arrays
that is designed to deal with this kind of problem.

\section{Requirements for an array descriptor}

\label{sec:requirements}

It seems that, in the translation of procedure calls like the one
discussed in the last section, some non-trivial data structure
must be passed to the translated procedure to describe the layout
of the actual argument.  This data structure is called a
Distributed Array Descriptor or DAD.  We wish to understand the requirements
on the DAD, and how best to organize it.  In view of the times, that
organization should probably be informed by object-oriented principles.

One of the most obvious structural features of an HPF array is that
it is a multi-dimensional entity.  Each dimension can be mapped essentially
independently by many different strategies including:
\begin{itemize}
\item
collapsed (serial),
\item
simple block distribution,
\item
simple cyclic distribution,
\item
block cyclic distribution,
\item
general block distribution\footnote{In HPF 2.0 this is a block distribution
with an arbitrary pattern of local block-sizes, specified by some vector.},
\item
through a linear alignment relation to some template dimension
with {\em any} of the above distribution formats.
\end{itemize}

In the translation of Figure \ref{fig:simpleMPI} there were various
formulae and code fragments that were special to the choice of simple
block distribution format.  These included the formula:
\small
\begin{verbatim}
      BLK_SIZE = (N + P - 1) / P
\end{verbatim}
\normalsize
for the size of locally allocated memory in terms of the array extent,
$N$, and the extent, $P$, of the processor arrangement; the
formula
\small
\begin{verbatim}
      BLK_START = R * BLK_SIZE
\end{verbatim}
\normalsize
for the least global index associated with the local block in
terms of the process coordinate $R$; the recipe
\small
\begin{verbatim}
      IF(N - BLK_START >= BLK_SIZE) THEN
          BLK_COUNT = BLK_SIZE
      ELSEIF(N - BLK_START > 0) THEN
          BLK_COUNT = N - BLK_START
      ELSE
          BLK_COUNT = 0
      ENDIF
\end{verbatim}
\normalsize
for the actual number of elements held locally\footnote{The value assigned
to {\tt BLK\_COUNT} could be expressed as:
\begin{displaymath}
\max(\min(\mbox{\tt BLK\_SIZE}, N - \mbox{\tt BLK\_START}), 0)
\end{displaymath}
}, and the formula
\small
\begin{verbatim}
      I = BLK_START + L
\end{verbatim}
\normalsize
which computes the global index in terms of the local loop index $L$.

All of these prescriptions can differ for other mapping strategies.
For example, if the array had been distributed in CYCLIC fashion the
formula for {\tt BLK\_START} would be simply
\small
\begin{verbatim}
      BLK_START = R
\end{verbatim}
\normalsize
the recipe for {\tt BLK\_COUNT} would be just
\small
\begin{verbatim}
      BLK_COUNT = (N - R + P - 1) / P
\end{verbatim}
\normalsize
(ie, $\lceil (N - R) / P \rceil$) and the global index would be given by
\small
\begin{verbatim}
      I = BLK_START + (P * L)
\end{verbatim}
\normalsize
These formulae may be checked against the second picture in Figure
\ref{fig:callexamples}.  If we are dealing with, say, block-cyclic
distribution, or a dimension with strided alignment, the corresponding
prescriptions get considerably more complicated.

The important thing is that we see evidence of a common set of operations
that must be performed on different kinds of distributed array dimensions.
The actual computations differ, depending on the distribution format and
alignment of the individual dimensions.  To the object-oriented programmer
it looks as if a class hierarchy is emerging, with the common operations
implemented as virtual functions.  In the organization of the DAD described
in this lecture this will be the {\tt Range} hierarchy---a series of
classes describing distributed index ranges.  The details will be
left to Section \ref{sec:ranges}.

Apart from ranges, there is another important kind of subcomponent in
our DAD.  We can explain the need for this component by considering the
last example in Figure \ref{fig:callexamples}, wich involved a rank-1
section of a rank-2 array.

In most respects this section should behave exactly like an ordinary
rank-1 array---this is a requirement of the Fortran language, and a
very convenient one.  It means, for example, that sections can be freely
passed to all the transformational intrinsic functions of Fortran.
A reasonable expectation is that the DAD for the rank-1 section should
contain just one {\tt Range} object describing its single dimension.
In the current example it would be most natural for that {\tt Range}
object to simply be a copy of the corresponding object in the parent
DAD.

But the implementation of the section is not {\em exactly} like an ordinary
one-dimensional array.  It retains some vestige of the orthogonal
dimension of its parent array.  In particular, in the example, the
section is only mapped to the top row of the processor arrangement {\tt Q}.
The translation of {\tt INIT} must somehow take account of the
possible existence of hidden ``higher dimensions'' of some parent array.
If they exist, the translated code should only execute the assignments on
processors in the appropriate row, or column---or more general slice---of
the processor arrangement.

In fact a related problem was addressed in our first example, when
we insisted on a test
\small
\begin{verbatim}
      IF (W_RANK < 4) THEN
          ...
      ENDIF
\end{verbatim}
\normalsize
to ensure that the local processor was a member of the processor
arrangement over which elements of the assignment variable were
distributed.  The solution of the current problem is to generalize
the idea of distributing an array over a processor arrangement,
and explicitly allow that an array can be distributed over some
{\em restricted slice} of a processor arrangement.  A class will be
introduced to represent this more general idea of a process group,
and an object from this class will be added to the DAD.  Then
the translated code just needs to check that the local process
is a member of the group, by a suitable method.  As we will see
in the next section, the required class has a compact and efficient
implementation.

\section{Groups}

Before giving a detailed specification of the interesting
{\tt Range} hierarchy we need a more concrete language for discussing
the processor arrangements over which index ranges will be distributed.
We will also take the opportunity to discuss representation of
DAD process groups in general.

We are familiar with the MPI idea of a process group, which is represented
by a data structure of type {\tt MPI\_GROUP}.  The groups referenced
by a DAD are conceptually similar, but likely to be implemented in a
different way.  If we are going to create and manipulate DADs with
the same freedom that typical data parallel programs create references
to array sections, for example, the group object in the DAD better
be very lightweight---highly streamlined.

\subsection{The process grid}

Any DAD group will have a parent processor arrangement.  The preceding
discussion about the DAD was supposed to be at a logical, design
level.  But from now on we will need to become more concrete.  So we
give specific interfaces in C++\footnote{Of course one can do equivalent
things in Fortran, but they look relatively clumsy and difficult to
read.}.  A C++ run-time representation of an HPF processor arrangement
would naturally be implemented as an object (just as a group or
communicator is logically an object in MPI).  We can easily see that
the HPF directive
\small
\begin{verbatim}
!HPF$ PROCESSORS P(4)
\end{verbatim}
\normalsize
might map to a C++ declaration something like
\small
\begin{verbatim}
    Procs1 p(4);
\end{verbatim}
\normalsize
and
\small
\begin{verbatim}
!HPF$ PROCESSORS Q(2, 2)
\end{verbatim}
\normalsize
might map to a C++ declaration something like
\small
\begin{verbatim}
    Procs2 q(2, 2);
\end{verbatim}
\normalsize
This supposes that {\tt Procs1} is a class of one-dimensional process
grids and {\tt Procs2} is a class of two-dimensional grids\footnote{In
MPI groups themselves do not have any grid structure.  On the other
hand a {\em Cartesian communicator} does have a grid structure.  So
one possible implementation of a {\tt Procs} object might be in terms
of this kind of communicator.}.  We will generalize to a family
(class hierarchy) of process grids illustrated in Figure \ref{fig:procs}.

\begin{figure}[tbp]
\centerline{\psfig{figure=procs.eps,height=1.2in}}
\caption{A hierarchy of process grids.}
\label{fig:procs}
\end{figure}

What methods are needed on the {\tt Procs} class?  Looking back
at Figure \ref{fig:simpleMPI} we see that we need the ability
to check if the local process is a member of the process grid
in question, and we need to find the coordinates of the local
process relative to the grid.  Later we will find that it is
useful to have a separate concept of {\em process dimension}---one
dimension of a process grid.  We introduce another
class---{\tt Dimension}---to stand for this concept.  Finding the
coordinate becomes a method on the dimension class.  Concrete
interfaces for {\tt Procs} and {\tt Dimension} are given
in Figure \ref{fig:procsAPI}.

\begin{figure}[tbp]
\small
\begin{verbatim}
    class Procs {
    public:
        int member() const ;

        Dimension dim(const int d) const ;
        ...
    } ;

    class Dimension {
    public:
        int size() const ;

        int crd() const ;
        ...
    } ;
\end{verbatim}
\normalsize
\caption{Interfaces of {\tt Procs} and {\tt Dimension}.}
\label{fig:procsAPI}
\end{figure}

If we decided we wanted to use these primitives in a C++ version
of our original translation example, the lines
\small
\begin{verbatim}
      INTEGER W_RANK, ...
      ...

      MPI_COMM_RANK(MPI_COMM_WORLD, W_RANK, ERRCODE)
      ...

      IF (W_RANK < 4) THEN
          BLK_START = W_RANK * BLK_SIZE
          ...
      ENDIF
\end{verbatim}
\normalsize
could be replaced with something like
\small
\begin{verbatim}
    Procs1 p(4);
    ...

    if (p.member()) {
        blk_start = p.dim(0).crd() * blk_size ;
        ...
    }
\end{verbatim}
\normalsize

A {\tt Procs} object is probably not particularly ``lightweight''.
It is likely to contain a vector of {\tt Dimension} objects
and probably some vector defining how logical processes are mapped
to physical processes, and so on.  Most likely a DAD should not
contain a {\em copy} of this information, but some {\em reference}
to an externally created {\tt Procs} object.

\subsection{Restricted process groups}

At the end of Section \ref{sec:requirements} we saw that in
fact the DAD should be able to reference a {\em restricted slice}
of a process grid.  Specifically what is needed is the ability to
reference a portion of a process grid selected by fixing the
coordinates in any subset of the dimensions to single values.
Figure \ref{fig:restricted} gives some examples.

\begin{figure}[tbp]
\centerline{\psfig{figure=restricted.eps,width=4in}}
\caption{Example restricted groups in a 2-dimensional process grid}
\label{fig:restricted}
\end{figure}

Luckily there is a very compact way to represent this kind of
subgroup.  If we call the horizontal dimension in the figure
{\tt dim(0)} and the vertical dimension {\tt dim(1)}, then
{\bf a)} represents the whole of the grid and has dimension
set $\{\mbox{\tt dim(0)},\mbox{\tt dim(1)}\}$,
{\bf b)} has dimension set $\{\mbox{\tt dim(0)}\}$,
{\bf c)} has dimension set $\{\mbox{\tt dim(1)}\}$, and
{\bf d)} has the empty dimension set $\{\}$.  With the
dimension set given, a group can be uniquely specified by
its lead process---the process with coordinates 0 relative
to the effective dimensions.  In our examples
{\bf a)} has lead process 0,
{\bf b)} has lead process 8,
{\bf c)} has lead process 1, and
{\bf d)} has lead process 6.

The dimension set can be represented as a subset of
the parent grid dimensions using a bitmask.  Since any reasonable
grid will have less than 32 dimensions, this bitmask will always
fit in a single word.  So we see that our 4 examples can be
parametrized by the tuples $(11_2, 0)$,
$(10_2, 8)$,
$(01_2, 1)$ and 
$(00_2, 6)$, respectively, together with a pointer to the {\tt Procs}
object representing the parent grid.

The {\tt Group} class will have an interface like
\small
\begin{verbatim}
    class Group {
    public:
        Group(const Procs& p) ;

        void restrict(Dimension d, const int coord) ;

        int member() const ;
        ...
    } ;
\end{verbatim}
\normalsize
and it can be implemented in about 3 words of memory.  So {\tt Group}
objects can be freely copied and discarded: they fit naturally in the DAD.
The conversion constructor provides a user-defined type conversion from
a process grid object to a group object representing the whole of
the grid.  The {\tt restrict()} method restricts the group in one of
its dimensions.

\section{Ranges}

\label{sec:ranges}

\begin{figure}[tbp]
\centerline{\psfig{figure=ranges.eps,height=2in}}
\caption{A hierarchy of ranges.}
\label{fig:ranges}
\end{figure}

Inside our DAD a range object will describe the extent and distribution
format of one dimension of a distributed array.  We can expect a class
hierarchy of different kinds of range object (Figure \ref{fig:ranges}).
Each subclass represents a different kind of distribution format for
an array dimension.  The simplest distribution format is {\em collapsed}
(sequential) format.  Other distribution formats include {\em regular block}
and {\em simple cyclic} decomposition.  We will see later that some
range classes allow {\em ghost extensions} to support stencil-based
computations.

\begin{figure}[tbp]
\small
\begin{verbatim}
    class Range {
    public:
        int size() const ;

        Dimension dim() const ;

        int volume() const ;

        Range subrng(const int extent, const int base,
                     const int stride = 1) const ;

        void block(Block* blk, const int crd) const ;

        void location(Location* loc, const int glb) const ;
        ...
    } ;
\end{verbatim}
\normalsize
\caption{Interface of the {\tt Range} class.}
\label{fig:rangeAPI}
\end{figure}

Figure \ref{fig:rangeAPI} gives a possible interface of a {\tt Range}
class in C++.  To understand this API, let's give an example of its
use.  Our original HPF example
and a possible translation to C++ as shown in Figure \ref{fig:simpleC++}.
We declare a process grid {\tt p} as described in the previous
section, then create a {\tt Range} object, {\tt x}, of subclass
{\tt BlockRange}.  The arguments of the constructor are the process dimension
over which the range is to be distributed, and the extent of the range.
The value previously called {\tt BLK\_SIZE} is returned by the method
call {\tt x.volume()}, and a local array segment is allocated with this
size.  In any process that is a member of {\tt p} we call the method
{\tt block}, passing it the local process coordinate.  This initializes
a {\tt Block} structure which has fields:
\small
\begin{verbatim}
    struct Block {
        int count ;

        int glb_bas ;
        int glb_stp ;

        int sub_bas ;
        int sub_stp ;
    } ;
\end{verbatim}
\normalsize
The fields define the count of the local loop and the base and step
for the local subscript and global index\footnote{In the translation of
Figure \ref{fig:simpleC++}, note that the expressions {\tt b.glb\_bas +
b.glb\_stp * l} and {\tt b.sub\_bas + b.sub\_stp * l} are ``induction
variables''---linear functions of the loop index.  The compiler
can compute these efficiently by constantly incrementing temporary
variables.}.

\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
    !HPF$ PROCESSORS P(4)

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

          FORALL (I = 1:50)  A(I) = 1.0 * I
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        ... do some global initialization

        Procs1 p(4);

        BlockRange x(p.dim(0), 50);

        float* a = new float [x.volume()];

        if (p.member()) {
            Block b;
            x.block(&b, p.dim(0).crd());

            for(int l = 0 ; l < b.count ; l++) {
                const int i = b.glb_bas + b.glb_stp * l + 1 ;
                a [b.sub_bas + b.sub_stp * l] = 1.0 * i ;
            }
        }

        ... do some global finalization
\end{verbatim}
\normalsize
\caption{C++ translation of simple HPF program.}
\label{fig:simpleC++}
\end{figure}

Figure \ref{fig:simpleC++} is shorter and (arguably) more readable
than the original translation in Figure \ref{fig:simpleMPI}.  It
also has the interesting property that if we change the distribution
directive in the HPF code to
\small
\begin{verbatim}
!HPF$ DISTRIBUTE A(CYCLIC) ONTO P
\end{verbatim}
\normalsize
the only change needed in the C++ translation is that the declaration
of {\tt x} becomes
\small
\begin{verbatim}
    CyclicRange x(p.dim(0), 50);
\end{verbatim}
\normalsize
This suggests we are making progress towards the goal of being able
to produce generic code that works for data with any distribution format.

\newpage

\section{A DAD}

Before ranges and groups can be put together to make a DAD we need to
deal with one other detail.  This is the issue of {\em memory strides}.

An array dummy argument may be matched with an array {\em section} in the
calling program.  So for a dummy the layout of elements in memory
may be different to what would be expected for a general Fortran array.
This problem is not specific to HPF---it is a feature
of Fortran 90.  In:
\small
\begin{verbatim}
      REAL D(100, 100)
      ...
      CALL FOO(D(1, :))


      SUBROUTINE FOO(C)
      REAL C(:)
      ...
      END
\end{verbatim}
\normalsize
the first dimension of {\tt D} is most-rapidly-varying, so the stride in
that dimension---the memory displacement between adjacent elements---is
1 word.  The second dimension has a stride of 100 words.  So the stride
for the dummy argument {\tt C} will also be 100.  Fortran compilers
typically take this into account by passing a {\em dope vector} for an
array argument.  For a rank-$r$ argument the dope vector contains the
$r$ extents and the $r$ strides that describe the shape of the array and
its memory layout \cite{Wolfe96}.  The procedure code uses the stride
information from the dope vector when resolving array references.

In our DAD the extent information fits naturally in the {\tt Range}
objects.  We find it convenient to keep the memory stride information
separate.  It seems useful to allow range objects to be shared by
arrays that have the same mapping in particular dimensions, although
they may have different overall shape and memory layout.

The upshot is that the DAD for a rank-$r$ array will contain:
\begin{itemize}
\item
A distribution group,
\item
$r$ range objects, and
\item
$r$ integer strides
\end{itemize}
Figure \ref{fig:DADAPI} gives an interface for the {\tt DAD} class.
Here {\tt Map} is a small structure that bundles a range with a stride:
\small
\begin{verbatim}
    struct Map {
    public :
        Map(Range _range, const int _stride) ;

        Range range ;
        int stride ;
    } ;
\end{verbatim}
\normalsize

\begin{figure}[tbp]
\small
\begin{verbatim}
    struct DAD {
        DAD(const int _rank, const Group& _group, Map _maps []) ;

        int rnk() const {return rank ; }

        const Group &grp() const ;

        Range rng(const int r) const ;

        int str(const int r) const ;
        ...
    } ;
\end{verbatim}
\normalsize
\caption{Interface of the {\tt DAD} class.}
\label{fig:DADAPI}
\end{figure}

With this DAD, we can finally give a translation of
the second example in section \ref{sec:introduction}.
Source and translation is given in Figure \ref{fig:inherittedC++}.
To fully appreciate how this translation works, we need to see how
the calling program sets up the DAD.  For this we revist the four
examples in Figure \ref{fig:callexamples}.  Sources and translations
are given in Figures \ref{fig:callblock} through \ref{fig:callrow}.

The translations in Figures \ref{fig:callblock} and \ref{fig:callcyclic}
are fairly self-explanatory.  A DAD is set up describing array {\tt b}.
This contains a group object describing {\tt p} (the conversion from
{\tt Procs} to {\tt Group} is implicitly used) and and a map object
containing the range {\tt x} and stride 1.  A pointer to the DAD is
passed to the procedure, together with the pointer to the elements.

Figure \ref{fig:callstride} is a little more complicated.  We again
create a range parameterizing the 100-element array {\tt b}.  In
mechanically translated code we might go on to set up a DAD describing 
the whole of {\tt B}, but for simplicity we skip this stage here.
We only set up a DAD describing the strided section of {\tt b}.
This involves using the {\tt subrng()} method on {\tt Range} to
create a new range object representing a subrange.  The subrange has
stride 2.  When the the procedure {\tt init()} calls the {\tt block()}
method on this object, it will return appropriately adjusted values for
subscript and global index bases and steps.

Figure \ref{fig:callrow} is more complicated again.
To deal with the scalar subscript it uses the {\tt location()}
method.  This takes the global subscript as argument and
returns a {\tt Location} structure, which has fields:
\small
\begin{verbatim}
    struct Location {
      int crd ;
      int sub ;
      ...
    } ;
\end{verbatim}
\normalsize
These define the corresponding process coordinate and local subscript.
In some sense {\tt location()} is the dual method to {\tt block()}---{\tt
block()} takes a process coordinate and returns a range of indices (and
local subscripts)---{\tt location} takes an index and returns a coordinate
(and a local subscript)\footnote{We assumed the argument of {\tt location()}
is a C-style index, starting from 0, rather than a Fortran-style index
starting from 1.}.  Note the memory stride for the second dimension,
{\tt x.volume()}, goes into the {\tt Map} object, and the base address
passed to the procedure is offset in memory by an amount determined by
the scalar subscript.  Of course the group in the DAD was suitably
restricted.

\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
          SUBROUTINE INIT(A)

          REAL A(50)
    !HPF$ INHERIT A

          FORALL (I = 1:50)  A(I) = 1.0 * I
          END
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        void init(float* a, DAD* a_dad) {

            Group p = a_dad->grp();
            if (p.member()) {
                Range x = a_dad->rng(0);
                int   s = a_dad->str(0);

                Block b;
                x.block(&b, x.dim().crd());

                for(int l = 0 ; l < b.count ; l++) {
                    const int i = b.glb_bas + b.glb_stp * l + 1 ;
                    a [s * (b.sub_bas + b.sub_stp * l)] = 1.0 * i ;
                }
            }
        }
\end{verbatim}
\normalsize
\caption{C++ translation of HPF program with inheritted mapping.}
\label{fig:inherittedC++}
\end{figure}


\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
    !HPF$ PROCESSORS P(4)

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

          CALL INIT(B)
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        ... do some global initialization

        Procs1 p(4);

        BlockRange x(p.dim(0), 50);

        float* b = new float [x.volume()];

        Map maps [1];
        maps [0] = Map(x, 1);

        DAD dad(1, p, maps);

        init(b, &dad);

        ... do some global finalization
\end{verbatim}
\normalsize
\caption{Translation of first example from Figure \ref{fig:callexamples}.}
\label{fig:callblock}
\end{figure}


\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
    !HPF$ PROCESSORS P(4)

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

          CALL INIT(B)
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        ... do some global initialization

        Procs1 p(4);

        CyclicRange x(p.dim(0), 50);

        float* b = new float [x.volume()];

        Map maps [1];
        maps [0] = Map(x, 1);

        DAD dad(1, p, maps);

        init(b, &dad);

        ... do some global finalization
\end{verbatim}
\normalsize
\caption{Translation of second example from Figure \ref{fig:callexamples}.}
\label{fig:callcyclic}
\end{figure}

\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
    !HPF$ PROCESSORS P(4)

          REAL B(100)
    !HPF$ DISTRIBUTE B(BLOCK) ONTO P

          CALL INIT(B(1:100:2))
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        ... do some global initialization

        Procs1 p(4);

        BlockRange x(p.dim(0), 100);

        float* b = new float [x.volume()];

        // Create DAD for section b(::2)

        Range x2 = x.subrng(50, 0, 2);

        Map maps [1];
        maps [0] = Map(x2, 1);

        DAD dad(1, p, maps);

        init(b, &dad);

        ... do some global finalization
\end{verbatim}
\normalsize
\caption{Translation of third example from Figure \ref{fig:callexamples}.}
\label{fig:callstride}
\end{figure}


\begin{figure}[p]
\textbf{SOURCE:}
\small
\begin{verbatim}
    !HPF$ PROCESSORS Q(2, 2)

          REAL B(6, 50)
    !HPF$ DISTRIBUTE B(BLOCK, BLOCK) ONTO Q

          CALL INIT(B(2, :))
\end{verbatim}
\normalsize
\textbf{TRANSLATION:}
\small
\begin{verbatim}
        ... do some global initialization

        Procs2 q(2, 2);

        BlockRange x(q.dim(0), 6),
                   y(q.dim(1), 50);

        float* b = new float [x.volume() * y.volume()];

        // Create DAD for section b(1, :)

        Location i;
        x.location(&i, 1);
        
        Group p = q;
        p.restrict(q.dim(0), i.crd)

        Map maps [1];
        maps [0] = Map(y, x.volume());

        DAD dad(1, p, maps);

        init(b + i.sub, &dad);

        ... do some global finalization
\end{verbatim}
\normalsize
\caption{Translation of fourth example from Figure \ref{fig:callexamples}.}
\label{fig:callrow}
\end{figure}


\bibliographystyle{plain}
\bibliography{translation}

\newpage

\appendix

\section{Other features of the Adlib DAD}

The DAD described in this lecture was a simplified version of
the Adlib descriptor developed at Syracuse during the PCRC project
\cite{Carpenter97}.  In fact we strictly followed the current version
of the Adlib spec, and the C++ examples are all executable using Adlib.
But in the interests of simplifying the presentation many fields and
methods were omitted.  In this appendix we will briefly summarize some
of the more important things that were missed out.

\subsection{Support for block cyclic distributions}

The approach to translating parallel loops given in this lecture
is too simple to deal with block-cyclic distribution.  In that case
the local loop needs an outer loop enumerating the set of locally-held
blocks.  Adlib provides some ``iterator'' classes ({\tt LocBlocksIndex},
etc) for enumerating these blocks.  It also has an extra field in the
{\tt Location} class to specify the block.  An {\tt offset()} method can
be used to compute the total memory offet implied by a {\tt Location}
object, taking account of this generalized case.

\subsection{Support for ghost extensions and other memory layouts}

Adlib supports ghost extensions on suitable ranges.  However the
local subscript field computed by {\tt block()} or {\tt location()}
does not incorporate the associated shifts in base address.  Also
Adlib supports a packed layout for arrays with strided alignment,
but again the local subscript field is not divided by the packing
factor.  Instead local subscripts or {\tt Location} objects must be
converted by methods called {\tt disp()}, {\tt offset()} to allow for
these layout options.  Steps in local subscripts must be converted
by a {\tt step()} method.  The need for the extra calls is slightly
unfortunate, but the advantages of having a universal definition of
the local subscript that is independent of these details of the memory
layout seem to warrant them.

\subsection{Support for loops over subranges}

The {\tt block()} method has variants that take triplet parameters to
efficiently enumerate local blocks of a {\em sub}range.  This is an
important optimization---more efficient than creating a subrange {\em
object} and using the simple {\tt block()} method on that.

Certain subranges only have blocks defined for a subset of coordinate
values, and there are methods to compute these ranges of coordinate
values.  As a matter of fact it is illegal to call {\tt block()} for
coordinate values outside these ranges---a trap for the unwary, but
a restriction that improves the efficiency of the implementation.

\subsection{Features to support communication libraries}

The components of the Adlib DAD have accumulated various other methods to
support the communication library.  Associated infrastructure includes
iterator classes ({\tt AllBlocksIndex}) for iterating over {\em all}
blocks of a range (not just the ones in a single process).

\subsection{Miscellaneous inquiries and predicates}

There are inquiry functions to discover many properties of ranges, groups
and DADs, and various predicate methods desighed for run-time checking
for correctness of programs.

\end{document}

