From: Rob Schreiber <[email-masked]>
Date: Mon, 15 Jan 1996 16:08:03 -0800

---------------------------------------------------------------------------
[email-masked] is a mailing list for discussion of control-parallel
features in HPF.  Instructions for adding or deleting yourself from this
list appear at the bottom of this message.
---------------------------------------------------------------------------

One more time.   Here is the reduction proposal WITH the external syntax,
and WITH examples of derived type reductions.   Please comment

\chapter{Proposed HPF Extension for Reduction Operations in Independent Loops}

    Rob Schreiber and Chuck Koelbel, and Joel Saltz

    January 15, 1996

\section{Changes since the last version}

Changes reflect response at the November, 1995 HPFF meeting.

There is a proposal for intrinsic reduction, followed by a proposed extension for
user defined reduction.   

The ``external" form in which the reduction variable is given as a
clause of the INDEPENDENT directive is described.  Now, the loop in
which a reduction variable is protected is labeled explicitly.  Each and
every occurrence of a reduction variable 
encountered while the loop is active must be
as a reduction variable in a reduction statement of the form specified
in Section~\ref{sec-redux}.
This alternative syntax has the advantages of economy (only one
directive regardless of the number of reduction statements) and clarity
(the protected loop is explicitly indicated).

\section{Rationale}

It is often the case that a data parallel computation cannot be
expressed in HPF Version 1.1 as an {\tt INDEPENDENT} loop because several of
the loop iterations update one or more variables.  In such cases
parallelism may be possible and desirable, because the order of updates
is immaterial to the final result.  This is most often the case with
accumulations, such as the following loop:
\CODE
        do i = 1, 1000000000
            x = x + complicated_function(i)
        end do
\EDOC
This loop can run in parallel as long as its iterations make their
modifications to the shared variable ({\tt X)} in an atomic manner.
Alternately, the loop can be run in parallel by making updates to
temporary copies of {\tt X}, with a (short) final phase to merge these
changes.  In either case, the computation is conceptually parallel but
not allowed by the original definition of {\tt INDEPENDENT}.

When programming this kind of computation to run in parallel in HPF
Version 1, the programmer must store update information in a temporary
array whose size is equal to the number of loop iterations, and then
use an intrinsic reduction function or {\tt XXX_SCATTER} library function
after the loop.  The problems with this approach are twofold:  the
temporary array may be excessively large; and the reduction operation
has to be intrinsic.

For these reasons, we have added a reduction feature to HPF that allows
loops that update shared variables (shared means that they are updated
by more than one loop iteration) to be labeled as independent.

\section{Technical proposal}

In this discussion, we use the term {\em range} (of a {\tt DO} construct) and
{\em active} (an active {\tt DO} construct) in their technical sense as defined
by Fortran 95.  
The range of a {\tt DO} construct consists of the statements
lexically bounded by the {\tt DO} statement and the terminal statement of the
construct.  Thus, a subprogram invoked from inside a {\tt DO} loop is not in
its range, although the {\tt CALL} statement or function reference that
invokes it is.    On the other hand, the {\tt DO} construct is active from
the time the {\tt DO} statement is encountered until the time that the
matching execution of the terminal statement is executed, or the
program stops, or control is transferred out of the construct in one of
the other ways allowed for {\tt INDEPENDENT} loops (see page 302 of the
Handbook and Section 4.4 of the HPF Standard.)

\subsection{Enhanced {\tt INDEPENDENT} directive}
The syntax of the INDEPENDENT directive is modified as follows
\BNF
independent-directive      \IS   !HPF$ INDEPENDENT [,new-clause] [,reduction-clause]
\FNB
The reduction clause lists the reduction variables that are active in the following
loop:
\BNF
reduction-clause           \IS  REDUCTION(reduction-var-list)

reduction-var              \IS  array-variable-name
                           \OR  scalar-variable-name
                           \OR  structure-component
\FNB

\begin{constraints}
\item
A reduction variable must be of intrinsic type.
\end{constraints}

Note that the syntax does not allow an array element, array section, or substring
to be designated as a reduction variable in the reduction clause of the INDEPENDENT
directive.    Even though such a subobject may occur in a reduction statement, it is the
entire array or character variable that is treated as a reduction variable.

Any variable whose name occurs in a reduction clause is said to be
{\em protected} while the immediately following {\tt INDEPENDENT DO} loop is active.
A protected variable may not be referred to at all, {\em except}  in a reduction statement
(see below).   In particular, it may not occur in any HPF directive, including the
variable list in a {\tt NEW} clause.   This includes any {\tt NEW} clause in the same
{\tt INDEPENDENT} directive.

\subsection{Reduction statement}
\label{sec-redux}

A reduction statement is an assignment statement of the following special form that  
occurs in the range (i.e. the lexical body) of an
independent {\tt DO} loop for which the name of its reduction variable occurs in a
reduction clause.
\BNF
reduction-stmt          \IS  

     reduction-var-ref = expr reduction-op reduction-var-ref
			\OR
     reduction-var-ref = reduction-var-ref reduction-op expr
			\OR
     reduction-var-ref = &
       reduction-function(reduction-var-ref, expr)
			\OR
     reduction-var-ref = &
       reduction-function(expr, reduction-var-ref)

reduction-var-ref   \IS  variable

reduction-op         \IS  intrinsic-op 

reduction-function   \IS        MAX
                           \OR  MIN
                           \OR  IAND
                           \OR  IOR
                           \OR  IEOR
                           \OR  MATMUL
\FNB
The operator or function that precedes or follows the reduction variable
reference after the equals sign is the reduction operator.

\begin{constraints}
\item
An assignment statement of the requisite form is a reduction statement if and only if
it occurs in the range of an {\tt INDEPENDENT DO} loop in which the name of its reduction-var
occurs in a
reduction clause on the {\tt INDEPENDENT} directive that precedes the {\tt DO}
statement.

\item  The two reduction-var-refs that occur in a
reduction statement must be lexically identical.

\item  A reduction-op must be an intrinsic, associative
operator; therefore, reduction-op may not be a rel-op; and it may not be
power-op (i.e. **).

\item  The following two forms of the reduction-stmt are not allowed:
\CODE
    reduction-var-ref = expr - reduction-var-ref

    reduction-var-ref = expr / reduction-var-ref
\EDOC

\item In strict HPF 2.0, the noncommutative intrinsics {\tt //} and {\tt MATMUL}
are forbidden in reduction statements.
\end{constraints}

The intrinsic operators and functions that are allowed are all
associative (at least in their mathematical definitions, even though the
usual implementations of the arithmetic operators by Fortran language
processors and the underlying hardware are not); they are listed below
(in Section~\ref{sec-implem}).  All are commutative except for character
concatenation and {\tt MATMUL}.  (We apply the terms associative and
commutative to the subtraction and division operators, since their
meaning is addition of the additive or multiplication by the
multiplicative inverse, respectively:  {\tt A - B} means {\tt A + (-B)}.)

{\bf  Note:
In HPF version 2.0, the associative but noncommutative intrinsic functions and operators
(i.e. MATMUL and concatenate) are disallowed in reductions.
The description of their use in reduction that occurs here is of an endorsed extension.
}

\subsection{Relaxation of the HPF standard's requirement for {\tt INDEPENDENT} loops}

In Section 4.4 of the HPF standard, 
the first proviso on what may not occur in the loop is modified to:

*  Any two operations that assign to the same object interfere with
each other, (begin added text) except that the allowed 
occurrences of a reduction variable in a
reduction statement do not cause interference.  (end added text)

Thus, the user's assertion implied by {\tt INDEPENDENT} about the behavior of
the program is weakened, and no longer applies to reduction variables
that are protected in the {\tt DO} construct.  

On exit from the outermost {\tt DO} loop in which it is protected
(the one whose INDEPENDENT directive contains a reduction clause naming it),
the value of a reduction variable is
processor dependent, although it is constrained by the model
implementation mechanism given in Section~\ref{sec-implem}.

\subsection{Example}

Here is an {\tt ADD_SCATTER} done as an independent loop:
\CODE
    !hpf$ independent, reduction(X)
    do i = 1, n
        X(index(i)) = X(index(i)) - f(i)
    enddo
\EDOC
The implementation will most likely make a local copy of an array of the
same shape as {\tt X}, and initialize it to zero.  Each iteration will
subtract the value of {\tt f(i)} from its local {\tt X(index(i))}.  To
create the final result, it will combine all the local arrays with the
initial value of {\tt X}.  The combining operator is addition, so that the
result is the sum of the initial value of {\tt X} and the local arrays.

\subsection{Restrictions on the 
lexical location of reduction statements, and on the
use of reduction variables}
\label{sec-location}

A variable that is updated by reduction statements in an independent
loop must be protected by explicit appearance in a reduction clause in
the INDEPENDENT directive for the outermost independent loop that
contains no NEW clause for that variable.

When a reduction operation is executed, some nest of {\tt DO} constructs will
be active; this is guaranteed by definition.
The reduction variable is protected while the loop in which it occurs in the reduction
clause of the INDEPENDENT directive is active.
If there are several nested INDEPENDENT DO loops surrounding the reduction statements
in which the variable is updated, which one is the right one to get the reduction
clause?  The answer is the outermost one, subject to the constraint that a reduction
variable may not appear in a NEW clause for that loop and may not occur at all
in that loop outside of reduction statements, in particular, it may not get the NEW
attribute.

On the other hand, if a statement that has the form of a reduction statement occurs
while an independent loop is active, but the updated variable is not a protected reduction
variable, then the programmer is guaranteeing that no two iterations of the independent
loop will update the same location.  For example:
\CODE
    !hpf$ independent
    do i = 1, n
        ! I know there are no repeated values in indx(1:n)
        x(indx(i)) = x(indx(i)) + f(i)
    end do
\EDOC

HPF requires that any reduction statement occur in the range of the {\tt DO}
construct in which its variable is protected: it must be lexically
contained in the body of the outermost loop in which it is protected.

A protected reduction variable may only occur in reduction statements,
and only in references to the left of the assignment operator and the
matching reference to the right of the assignment operator.  Outside of
the loop in which it is protected, however, it is treated as an ordinary
variable.

\subsection{Example loop nests}

When loops are nested, a reduction variable may need to be protected in an
independent outer loop even though the reduction operations in which it
occurs are nested inside an inner loop.  Moreover, the inner loop and
any intervening loops may or may not be independent.
\CODE
!   Nested Loop Example 1.  Inner loop is sequential

      x = 10
      outer: do while (x < 1000)   ! This loop is sequential
         !hpf$ independent, new(j), reduction(x)
         middle: do i = 1, n       
           inner: do j = 1, m
             x = x + j
             !  Note that it would be incorrect
             !  to refer to x here, except in a reduction statement
           enddo inner
           !  Note that it would be incorrect
           !  to refer to x here, except in a reduction statement
         enddo middle
         print *, x
      enddo outer
\EDOC
Since the variable {\tt x} occurs in a reduction clause for loop {\tt middle},
it is a protected reduction variable
throughout that loop, including inside the inner loop.   
If {\tt inner} was an independent loop,
it would be incorrect to include {\tt x} in a reduction clause of its INDEPENDENT
directive, either in addition to or instead of loop {\tt middle}.   

The outermost
loop is not {\tt INDEPENDENT}, and so {\tt x} need not and cannot be protected 
in that part of its range outside the middle loop.

This remains true if the inner loop is independent.   On the other
hand, a variable that occurs in a {\tt NEW} clause may not be a reduction
variable in the range of a loop, although it may be in the range of a
contained loop:
\CODE
!   Nested Loop Example 2.  Outer loop NEW clause.

      !hpf$ independent, new(i)
      outer: do k = 1, 100
         !hpf$ independent, new (j,x)
         middle: do i = 1, n
            x = 10
            !hpf$ independent, reduction(x)
            inner: do j = 1, m
              x = x + j**2
              !  Note that it would be incorrect
              !  to refer to x here, except in a reduction statement
            enddo inner
            y(i) = x
         enddo middle
      enddo outer
\EDOC
Here, {\tt x} is a protected reduction variable only in the inner loop.

The reduction statement may not occur in a subprogram called from the range of the
independent loop.   To get this effect, return the updating value from the
subprogram and perform the update with a reduction statement in the lexical loop
body.
\CODE

        program example
        real y, z(10)
        integer x, xinc, i
        
        x = 5
        !hpf$ independent, new(y), reduction(x)
        do i = 1, 10
            x = x + update_with_side_effect(i)
        enddo

        contains
        integer function update_with_side_effect(i)
          integer i, local_b, j, k

          !  Must return the increment to x
          !  A reduction statement for x may not appear inside sub
          !  local_b is implicitly NEW on each outer iteration.
          !
          local_b = (i-1)**2
          xinc = i * local_b

          !hpf$ independent, reduction(y)
          do j = 1, 10
             y = y + 1
          enddo

          !  Sequentially update a different part of z
          !  No reduction or new clause needed for z
          !
          do k = 1, 10
             z(i) = z(i) + k*y
          enddo
        end subroutine sub
	end program
\EDOC
\subsection{Reduction operators may not be mixed}

In most cases, only one operator will be used in the reduction
statements (if there are more than one) that update a given reduction
variable.   It is sensible, however, to use different operators, such
as + and -, together on the same reduction variable: mathematically,
subtraction is just addition of the additive inverse.  The same is true
for multiplication (*) and division (/).  No other mixing of operators
is allowed, however.

\section{Implementation and semantics}
\label{sec-implem}

The result of the Fortran 95 intrinsic function  {\tt SUM}  is defined to
be a processor dependent approximation to the sum of the elements of
its argument array.  On exit from the loop in which it is
protected, the value of a reduction variable is likewise not completely
specified by HPF.  One possible value is that which would have been
computed by sequential execution of the loop, but other,
processor-dependent approximations to this value may be produced.

While HPF does not explicitly define
the set of possible values, it does give a precise implementation
mechanism that serves to define this set.

Since no reference to a protected reduction variable can occur except in 
a reduction statement, it is not necessary to define
the values that these variables may have while protected.

In this discussion, the term processor means a single physical
processor or a group of physical processors that together sequentially
execute some or all of the iterations of an independent loop.

We first give a simple implementation mechanism that applies to
commutative reduction operations.  
On entry to an independent loop,
every executing processor allocates a local
reduction variable associated with each variable in the reduction clause on the 
{\tt INDEPENDENT} directive, and initializes it to the identity element for the
intrinsic reduction operator.  The local reduction variable has the
same type and kind type parameter
as the global reduction
variable.  Its shape is determined by the shape of the global reduction variable,
and is identical to it, except as noted below.

The identity elements for the intrinsic operators are as
follows:

\begin{table}
\begin{verbatim}
        Operator      Identity Element
        -------       ----------------
        +                0
        -                0
        *                1
        /                1
        .and.            .true.
        .or.             .false.
        .eqv.            .true.
        .neqv.           .false.

        //               ''

\end{verbatim}
\caption{Identity elements for intrinsic reduction operators.}
\end{table}

The following intrinsics may be used as reduction-functions, with the indicated
identity elements:

\begin{table}
\begin{tabular}{l@{~~~}p{3in}}
        Operator    &   Identity Element \\
\hline
        {\tt IAND(I,J)}   &   (all ones)\\
        {\tt IOR(I,J)}    &   0         \\
        {\tt IEOR(I,J)}   &   0         \\
        {\tt MIN(X,Y)}    &   {\tt -HUGE(X)} where {\tt X} 
                              has the same type and kind type
                              parameter as the reduction variable\\
        {\tt MAX(X,Y)}    &   {\tt +HUGE(X)} where {\tt X} has the same 
                        type and kind type
                        parameter as the reduction variable\\
        {\tt MATMUL(A,B)} &   {\tt MATMUL} is noncommutative; 
                        it is the only allowed
                        associative and noncommutative intrinsic reduction function.
                        As a consequence of the form {\tt A = MATMUL(A, B)}
                        of the reduction operation, it must be
                        true that {\tt SIZE(A,2) = SIZE(B,1) = SIZE(B,2)}.
                        As a consequence of the alternate form 
                        {\tt A = MATMUL(B, A)}
                        of the reduction operation, it must also be
                        true that {\tt SIZE(A,1) = SIZE(B,2) = SIZE(B,1)}.  Thus,
                        the identity is the identity matrix of shape 
                        {\tt (/N, N/)},
                        of the same type and kind type parameters as the reduction
                        variable; for the left identity, {\tt N = SIZE(A,1)}, for
                        the right identity, {\tt N = SIZE(A,2)}.   \\
\end{tabular}
\caption{Identity elements for intrinsic reduction functions.}
\end{table}

Let {\tt L} be the local reduction variable, {\tt A} the global reduction variable.
For array operators other than {\tt MATMUL}, 
the shape of {\tt L} will be the same as the shape of {\tt A},
except that if all the expr's in reduction statements for {\tt A} are 
scalar, then a
scalar {\tt L} may be used.   
(In the case of {\tt MATMUL}, {\tt L} will have shape
determined by that of {\tt A}, as described in Table 2.)

As an example, consider:
\CODE
        double precision y(100)
        real x(100), z(100)

        z = 5.
        x = (/ (real(i), i = 1:100) /)
        !hpf$ independent, new(y), reduction(x, z)
        do i = 1, 10
           z = i + z
           y = ...
           x = x - y
        enddo
\EDOC
The local reduction variables will be real arrays {\tt local_x} and 
{\tt local_z} of shape {\tt (/ 100 /)}.   All of these are initially zero.
An optimizing compiler will use a scalar {\tt local_z},
one per processor.
The double precision  values (of {\tt local_x + y}) will be rounded to single precision 
after every local addition.    

Each processor performs a subset of the loop iterations;  when it
encounters a reduction statement, it updates its own copy of the
reduction variable.    A processor is free to perform its loop
iterations in any order; furthermore, it may start an iteration,
suspend work on it, do some or all of the work of other iterations, and
resume work on the suspended iteration.   However, any update of a
local reduction variable occurs through the execution of a reduction
statement, and reduction statements {\em are} executed atomically.

The final value of the reduction variable is computed by combining the
local values with the value of the global reduction variable on entry
to the loop, using the combining operator.  
The ordering of this
reduction is language processor dependent, just as it is for the
intrinsic reduction functions ({\tt SUM}, etc.)   
Thus, in the example above, the final value of {\tt z} will be (approximately)
60.

Next, we give an implementation for the noncommutative operator {\tt MATMUL}.

First, assume that a given processor executes a contiguous set of the
loop iterations -- i.e. there is a block mapping of loop iterations to
processors.  In this case, the processor allocates a local reduction
variable to accumulate all updates applied from the right  (as in 
{\tt X = MATMUL(X, <expr>))} and initializes it to the right identity element.  it
allocates another local reduction variable to accumulate all updates
applied from the left  (as in {\tt X = MATMUL(expr, X))} and initializes it
to the value of the left identity element.  It may omit either of these
if there are no corresponding updates in the loop.

When all iterations of the loop have been completed, the local
reduction variables are combined to produce the final result.   The
combination is an evaluation of the expression
\CODE
(1)   lrv(P) <op> lrv(P-1) <op> ... <op> lrv(1) 
                     <op> A <op> 
                            rrv(1) <op> ... <op> rrv(P)
\EDOC
where A is the initial value of the reduction variable, lrv(k), k=1,
..., P is the value of the local reduction variable that has
accumulated the updates from the left due to the kth group of loop
iterations (on processor k, in the case of a block mapping of
iterations to processors), and rrv(k), k = 1, ..., P is the value of
the local reduction variable that has accumulated the updates from the
right due to the kth group of loop iterations.  In making this
evaluation, the language processor is free to employ any associated
rearrangement of the given expression; i.e. it is not required to
evaluate the expression left to right, but may use any parenthesization
of it.   It may not interchange the order of the expressions, since
the operator is not commutative.

If the mapping of iterations to processors is not a block mapping, then
a left, a right, or both a left and a right reduction variable must be
created and used for every contiguous group of loop iterations, with
updates accumulated and combined as above; in expression (1), P is the
number of such contiguous groups of iterations.

In the endorsed extension to HPF for noncommutative intrinsic reduction, updates
may occur on either the left or the right, but not both, in a given independent loop.

\section{Examples}

The reduction variable may be a scalar array element:
\CODE
      !hpf$ independent, reduction(x)
      do i = 1, n
          x(indx(i)) = x(indx(i)) + <expr>
      enddo
\EDOC
It may be to a section of an array:
\CODE
      !hpf$ independent, reduction(x)
      do i = 1, n
          x(i: i+4) = x(i: i+4) + <expr>   ! as many as 5 updates to elements of x
      enddo
\EDOC
Note that the compiler, if it distributes iterations of this loop in a block-wise
manner, will not need to make a private copy of the entire array x.

\section{Extension to user-defined operators}
This section is a separate and detachable part of the proposal.

\subsection{{\tt COMMUTATIVE} attribute for functions}

First, allow {\tt COMMUTATIVE} as an addition attribute of functions:
\CODE
        PURE, COMMUTATIVE FUNCTION FOOBAR(A,B)
\EDOC
The user asserts that the value of {\tt FOOBAR(A,B)} is equal to the value
of {\tt FOOBAR(B,A)} for all inputs A and B.

If it is desired to make no change to the language, but simply add an assertion directive,
then this can be done instead.

\subsection{Initial value for local variables in reduction clause}

The reduction clause is allowed to designate an initial value
for the local reduction variable:
\CODE
!hpf$ independent, new(B), reduction (A = MY_OP_IDENT(SIZE(B,1)))
do i = 1, n
   B = ...
   A = MY_OP(A, B)
enddo
\EDOC
The reason to assign the identity to the reduction directive rather than
make it an attribute of the function is the difficulty of restricting
the kind of expression allowed in the later case.  In this proposal, any
expression whose value is conformable with the reduction variable is
permitted.  Note that the expression is evaluated only once on entry to the
loop nest, and its value is conceptually replicated, so that all
processors initialize their local reduction variables to this value.
Thus, the expression could use a non-pure function.  We might disallow
this, in order to permit concurrent, one-per-processor evaluations of
the identity expression.  

In this proposal, all restrictions that apply to intrinsic reduction apply equally.
Thus, reduction functions may not be mixed.   If a reduction function is applied
``from the right" (as in the example) then it may not be applied from the left in the same
loop.

\subsection{Examples of derived type reduction}
First, one may use extended precision arithmetic:
\CODE

    interface operator (+)
       pure commutative type (bignum) function big_plus(x,y)
       type (bignum) x,y
       end function big_plus
    end interface

    type (bignum) b, big_zero
    external big_zero

    !hpf$ independent, reduction(b = big_zero())
    do i = 1, n
        b = b + <bignum-expr>
    enddo
\EDOC

In the following example the user decides to accumulate sparse local updates to
a distributed dense matrix.

\CODE
      module sparsity
      type real_sparse_matrix
          integer nrows, ncols, nnz ! shape, number of nonzeros
          real, pointer :: vals(:)
          integer, pointer :: row_index(:), col_index(:)
      end type real_sparse_matrix

      interface operator (+)
        pure function full_plus_sparse(x,y)
        real x(:,:), full_plus_sparse(size(x,1), size(x,2))
        type (real_sparse_matrix) :: y
        end function full_plus_sparse

        !  THIS IS THE REDUCTION FUNCTION
        pure commutative function sparse_plus_sparse(x,y)
        type (real_sparse_matrix) :: x, y, sparse_plus_sparse
        !hpf$ commutative, &
              identity (spzeros(x%nrows, x%ncols)) :: sparse_plus_sparse
        end function sparse_plus_sparse
      end interface

      contains
      pure function spzeros(nr, nc) result(res)
        type (real_sparse_matrix) res
        integer nr, nc
        res%nrows = nr;   res%ncols = nc;   res%nnz = 0
        allocate(res%vals(0))
        allocate(res%row_index(0))
        allocate(res%col_index(0))
      end function spzeros

      pure function elt(i, j, x) result(res)
        type (real_sparse_matrix) res
        integer i,j
        real x
        res%nrows = nr;   res%ncols = nc;   res%nnz = 1
        allocate(res%vals(1))
        allocate(res%row_index(1))
        allocate(res%col_index(1))
        res%vals(1) = x
        res%row_index(1) = i
        res%col_index(1) = j
      end function elt
      end module

      program test
      use sparsity
      real gs(n,n)
      type (real_sparse_matrix) :: s

      s = spzeros(n, n)
      !hpf$ independent, new(row, col, val), reduction(s = spzeros(n, n))
      do i = 1, num_updates
         s = s + elt(row, col, val) ! adds in val at position (row, col)
      enddo
      gs = full_plus_sparse(gs, s)
      end program

\EDOC

---------------------------------------------------------------------------
To (un)subscribe to this list, send mail to [email-masked].
Leave the subject line blank, and in the body put the line
(un)subscribe <email-address>
---------------------------------------------------------------------------