
\section{Data-parallel computation in Fortran 90}

\subsection{Some array terminology}

     In Fortran, any variable or expression is either {\em scalar},
or has a non-zero {\em rank}.  In the second case we are dealing with
an array; the ``rank'' is the number of dimensions to that array.
\begin{verbatim}
  REAL A (5, 10), B (11:15, 10), C (11:15, 11:20)
\end{verbatim}
declares \verb$A$, \verb$B$ and \verb$C$ as rank-2 arrays (of real
numbers).  In Fortran the index of an array starts at 1 by default, but
index ranges with other starting points can be specified.  This has
been done for the first dimension of \verb$B$ and both dimensions of
\verb$C$.  Although the indices associated with the first dimension of
\verb$A$ and \verb$B$ take different values (1 to 5 for \verb$A$; 11 to
15 for \verb$B$), the {\em extent} they range over is 5 in both cases.
Similarly, the extent of the second dimension is 10 for all three
arrays.

     The vector of extents associated with an array is called the {\em
shape} of the array.  \verb$A$, \verb$B$ and \verb$C$ all have shape
$(5, 10)$.  Two arrays are {\em conformable} if they have the same
shape.  So \verb$A$, \verb$B$ and \verb$C$ are certainly conformable.

\subsection{Fortran 90 array expressions}

     The most primitive kind of array expression is the
``array constructor''.  This is a list of items in special
brackets.  For example
\begin{verbatim}
  (/ 2, 3, 5, 7 /)
\end{verbatim}
is a constant array of shape $(4)$.  We won't have a great deal to
say about array constructors, but two things are worth mentioning:
lists inside array constructor brackets can include implied do loops;
and array constructors only yield one-dimensional arrays, but
the intrinsic function \verb$RESHAPE$ can be used to obtain higher rank array
values.

     Turning to arithmetic, any of the built-in operations of Fortran 90 can
be applied directly to conforming arrays.  If \verb$A$, \verb$B$ and
\verb$C$ are declared as in the previous section,
\begin{verbatim}
  A + B
  B * C
  A + (B * C)
\end{verbatim}
are legal expressions.  Their results are arrays of the same shape
as the operands.  \verb$B * C$ stands for the table of values
\begin{center}
\makebox{\psfig{figure=b_times_c.ps,height=1in}}
\end{center}
Values held in corresponding positions in the arrays
are multiplied together.  In general, as here, corresponding
``positions'' may be associated with different indices, if the index
ranges are declared differently.

     As a special case, scalars are ``conformable'' with any array.
If \verb$D$ is scalar
\begin{verbatim}
  REAL D
\end{verbatim}
then \verb$C + D$ represents
\begin{center}
\makebox{\psfig{figure=c_plus_d.ps,height=1in}}
\end{center}

     In addition to arithmetic operations, many of the Fortran intrinsic
functions can be applied ``elementally''.  Hence \verb$COS (A)$ means
the array of values
\begin{center}
\makebox{\psfig{figure=cos_a.ps,height=1in}}
\end{center}
This possibility doesn't extend to user-defined functions in Fortran 90,
but see section \ref{pure}.

     Array expressions (and variables) can also be formed by taking a
{\em sections} of a named array object.  An array section is built from
some subset of the elements of an array object---those associated with
a selected subset of the index range attached to the object.
(Incidentally, the restriction to ``named object'' is significant.  In
Fortran index ranges are associated with declared, named arrays.  It
isn't possible to take an element or section of an unnamed expression,
even if it is array-valued.)

     The simplest array section is formed by subscripting an array by
colon-separated range of index values.  Suppose an array \verb$X$ is
declared by
\begin{verbatim}
  INTEGER X (0:12)
\end{verbatim}
Then \verb$X (1:3)$ is a section containing the second, third
and fourth elements of \verb$X$ (starting at the second element of
\verb$X$, because indexing of \verb$X$ starts at 0). 
It is an array-valued expression of shape $(3)$, and can be used
in the same context as any other such expression.  Either the
lower bound or the upper bound can be omitted in the index
range of a section, in which case they default to the lowest or
highest values taken by the array's index.  So \verb$X (:)$
is a section containing the whole of \verb$X$.

     This simple form of array section can be generalised with a
``stride''.  \verb$X (1:10:3)$ is an array section of size 4 containing
every third element of \verb$X$ with indices between 1 and 10 (ie,
indices 1, 4, 7, 10).  Even more generally, if \verb$Y$ is a one
dimensional integer array, \verb$X (Y)$ is an array section of the same
length as \verb$Y$, built from elements of \verb$X$ by taking each
element of \verb$Y$ in turn as an index into \verb$X$.

     For multi-dimensional arrays, some dimensions could be subscripted
with a normal scalar expression, and some could be ``sectioned''.  The
rank of the resulting array variable is the number of
dimensions that are sectioned.

\subsection{\label{assignment}Array assignments}

Naturally, we can use array-valued expressions
in assignments.  The simplest version has exactly the same syntax as
the scalar Fortran assignment.  The expression on the right hand side
of the assignment must be conformable with the array variable on 
the left hand side.

     As a first, elementary example of a ``data parallel'' program,
let's take
\begin{verbatim}
  REAL HORIZ (10), VERT (10), HYPOT (10)

  HYPOT = SQRT (HORIZ**2 + VERT**2)
\end{verbatim}
Assume HORIZ and VERT were initialised to the lengths of the
horizontal and vertical sides of ten right-angled triangles.  The result,
returned in the array HYPOT is the length of the ten
hypotenuses\footnote{We need to check that the expression on the right
hand side of the assignment is a legal array expression, according to
the rules given above.  HORIZ is an array of shape $(10)$.  The
constant $2$ is a scalar, conformable with any array, and \verb$**$ is
a built-in arithmetic operator which can be used to combine any
conformable expressions.  So the subexpression \verb$HORIZ**2$ is a
legal array expression of shape $(10)$.  The same goes for the
sub-expression \verb$VERT**2$, and the sum of these two expressions.
Finally \verb$SQRT$ is an ``elemental'' intrinsic function.  The full
expression is conformable with the array variable on the left hand
side, so we have a legal program.}.

     The effect of an array assignment is that first all elements
of the expression on the right-hand side are evaluated, {\em then}
they are assigned to the elements of the variable on the left-hand-side.
This statement is than less vacuous than it sounds.  For
example the statement
\begin{verbatim}
  X (0:9) = X (-1:8)
\end{verbatim}
replaces all elements of \verb$X$ with indices between 0 and 9
by the values of the preceding elements.  It has a different
effect from the loop
\begin{verbatim}
  DO i = 0, 9
    X (I) = X (I - 1)
  ENDDO
\end{verbatim}
which replaces all the elements of \verb$X$ with indices between
0 and 9 with the value of \verb$X (-1)$.

     Returning to parallel programming with a less trivial problem,
let's try to sum the elements of an array, using a parallel algorithm.
\begin{verbatim}
  INTEGER INPUT (16)
  INTEGER SUM, I

  SUM = 0
  DO I = 1, 16
    SUM = SUM + INPUT (I)
  ENDDO
\end{verbatim}
is a sequential program which achieves the required result.  It takes 16
sequential steps to get there.  A more parallel approach is to sum the
16 numbers pairwise in parallel, to give 8 results.  These 8 are
then summed pairwise in parallel to give 4 numbers, and the process is
repeated until a single result is obtained.
This can be written as follows
\begin{verbatim}
  INTEGER INPUT (16), TEMP1 (8), TEMP2 (4), TEMP3 (2)
  INTEGER SUM
  
  TEMP1 = INPUT (1:8) + INPUT (9:16)
  TEMP2 = TEMP1 (1:4) + TEMP1 (5:8)
  TEMP3 = TEMP2 (1:2) + TEMP2 (3:4)
  SUM   = TEMP3 (1) + TEMP3 (2)
\end{verbatim}
Four parallel steps reduce the set to a single result.  In principle
this could save a factor of $16/4 = 4$ in execution time.
The program can be generalised with a loop
\begin{verbatim}
  INTEGER INPUT (16), TEMP (16)
  INTEGER SUM, I

  TEMP = INPUT
  I    = SIZE (INPUT)
  DO WHILE I > 1
    J = I / 2
    TEMP (1 : J) = TEMP (1 : J) + TEMP (J + 1 : I)
    I = J
  ENDDO
  SUM  = TEMP (1)
\end{verbatim}
(\verb$SIZE$ is a new intrinsic function in Fortran 90, which can
be used to find the extent of an array).
This will work for any \verb$INPUT$ array whose size is a power of two.
It gives a result in $\log_2 (\mbox{\verb$SIZE (INPUT)$})$
iterations.

     A second form of array assignment allows one to restrict the
assignment to some subset of the array variable on the left-hand-side
by specifying a ``mask''.  A mask is \verb$LOGICAL$ array expression,
conformable with the variable.  Assignment of the
corresponding element of the array expression on the right-hand is only
performed where the mask takes value \verb$.TRUE.$.  Where the mask
is \verb$.FALSE.$, the variable is left unchanged.  For example,
suppose \verb$EVEN$ is a logical array with indices from 0 to 9, whose
even-indexed elements are \verb$.TRUE.$ and odd-indexed elements are
\verb$.FALSE.$.  Then
\begin{verbatim}
  WHERE EVEN
    X (0:9) = X (1:10)
  ENDWHERE
\end{verbatim}
replaces the values in the section \verb$X(0:9)$ having even indices
with copies of the values of the next element in the array.
\verb$WHERE$ constructs can have optional \verb$ELSEWHERE$ clauses, which look
rather similar to \verb$ELSE$ clauses in \verb$IF$ constructs.  This is
often useful.  It can have slightly surprising effects if variables
assigned in the main clause are referenced in the \verb$ELSEWHERE$
clause.  Take the code
\begin{verbatim}
  WHERE EVEN
    X (0:9) = X (1:10)
  ELSEWHERE
    X (0:9) = X (-1:8)
  ENDWHERE
\end{verbatim}
We might hope that the effect is to replace each
even-indexed element of \verb$X$ with the next element of the array, and
replace each odd-indexed element with the previous element, ie to swap
adjacent odd and even elements, $X (0) \leftrightarrow X (1)$, $X (2)
\leftrightarrow X (3)$, etc.  {\em In fact}, the \verb$WHERE$/\verb$ELSWHERE$
construct is defined to make the above program equivalent to
\begin{verbatim}
  WHERE EVEN
    X (0:9) = X (1:10)
  ENDWHERE
  WHERE .NOT. EVEN
    X (0:9) = X (-1:8)
  ENDWHERE
\end{verbatim}
which has a different effect---the first assignment copies the following
odd element down into each even element, but the second assignment
copies the {\em new} value back to where it originally came from.
So the second assignment effectively does nothing.

     In general both the \verb$WHERE$ and \verb$ELSEWHERE$ clauses
can contain any number of assignments, so long as all of the
array variables are conformable.  The meaning is defined by splitting
the construct into a series of consecutive single-assignment
\verb$WHERE$ constructs, as above.

\subsection{\label{array-args}Arrays as arguments and results of subprograms}

     Consider the following subprogram definition
\begin{verbatim}
  FUNCTION MATPROD (A, B)
  REAL A (:, :), B (:, :)
  REAL MATPROD (SIZE (A, 1), SIZE (B, 2))
    DO I = 1, SIZE (A, 1)
      DO J = 1, SIZE (B, 2)
        MATPROD (I, J) = DOT_PRODUCT (A (I, :), B (:, J))
      ENDDO
    ENDDO
  ENDFUNCTION
\end{verbatim}
This example introduces several new features of Fortran 90.
Firstly, an {\em assumed
shape} dummy argument is an array whose index (upper) bounds
are undefined, as for \verb$A$ and \verb$B$ (a lower bound can be
optionally specified, but it defaults to 1, as here).  The upper
bounds are set
according to the shape of the actual argument.  Secondly \verb$SIZE$
has an optional argument that defines which dimension of a
multi-dimensional argument is required, so \verb$SIZE (B, 2)$ is the
extent of the second dimension of \verb$B$.  Thirdly \verb$DOT_PRODUCT$
is a new intrinsic function which returns the scalar product of two
vectors.

     Last but not least, the result of a function can be an array.
Furthermore, the size of this array can depend on the arguments, as in the
example.

     In various circumstances, including the case when a dummy
argument is assumed size, and when a function returns an array result,
the called procedure must have an {\em explicit interface} in
the calling program.  There are various ways to make a subprogram's
interface explicit.  The simplest is to include an {\em interface
block} amongst the declarations in the calling program.  A
program calling \verb$MATPROD$ might contain the declaration
\begin{verbatim}
  INTERFACE
    FUNCTION MATPROD (A, B)
    REAL A (:, :), B (:, :)
    REAL MATPROD (SIZE (A, 1), SIZE (B, 2))
    ENDFUNCTION
  ENDINTERFACE
\end{verbatim}
This contains just the information (the attributes of the arguments
and result) needed for the compiler to generate the correct calling
code.  Typically, one might collect groups of interfaces in
text files, and \verb$INCLUDE$ such files in programs calling
procedures from a group.  In full Fortran 90, the module system
incorporates a more controlled way of making interfaces explicit.

\subsection{Fortran 90 array intrinsics}

     The set of Fortran intrinisic functions has been extended to
include many array operations.  One of these functions
(\verb$DOT_PRODUCT$) was introduced in the previous section.

     Some other notable inclusions are the ``array reduction''
functions including \verb$SUM$ and \verb$PRODUCT$ which
add or multiply together elements of an array; \verb$MAXVAL$ and
\verb$MINVAL$ which return largest or smallest elements of
an array; and \verb$ANY$ and \verb$ALL$, which take the
disjunction and conjunction of the elements of a \verb$LOGICAL$
array.

     Other functions include the \verb$SPREAD$ function which ``adds an
extra dimension'' to an array, returning an array in which values from
the original array are replicated over the range of the new index;
\verb$TRANSPOSE$, which returns an array containing in the elements of
its two dimensional argument, transposed; and shift operations which
return arrays with values shifted by a constant index offset in one of
the dimensions of the argument array.

     We can assume (or hope) that these functions have been implemented
using the best available---where appropriate, parallel---algorithms.
So, for example \verb$SUM$ may have been implemented along the lines
described in section \ref{assignment}.

     Creative use of the array manipulation functions greatly
extends the class of parallel algorithms we can express in Fortran 90.
But in many cases calls to functions like \verb$SPREAD$, \verb$TRANSPOSE$
and the shifts can be eliminated, and readability of programs
improved, by using the HPF construct \verb$FORALL$.

