Subject: hpff-task: reductions

From: Rob Schreiber <[email-masked]>
Date: Tue, 2 Jan 1996 12:13:26 -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.
---------------------------------------------------------------------------

Here is a new rev of the reduction proposal.   Note that except for the last
section, only intrinsic reduction is considered.

    Proposed HPF Extension for Reduction Operations in Independent Loops

    Rob Schreiber and Chuck Koelbel, and Joel Saltz

    January 2, 1996

0.  Changes Since the Last Version.

RSS is now 46 years old.

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.

1.  Rationale

It is often the case that a data parallel computation cannot be
expressed in HPF Version 1.1 as an 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:

        do i = 1, 1000000000
            x = x + complicated_function(i)
        end do

This loop can run in parallel as long as its iterations make their
modifications to the shared variable (X) in an atomic manner.
Alternately, the loop can be run in parallel by making updates to
temporary copies of 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 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 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.

2.   Technical Proposal

2.1 Add a REDUCE Directive

H499 reduce-directive         is  !HPF$ REDUCE

Constraint:  The first non-comment, non-blank line following a
reduce-directive must be an assignment-stmt of a special form, called a
reduction statement.

2.2 Reduction Statment

A reduction statement is an assignment statement of a special form, that
occurs in the range of an independent DO loop.

H498 reduction-stmt           is  assignment-stmt

Any statement immediately preceded by a REDUCE directive is a reduction
statement.  An assignment to a whole scalar or array variable in the
range of an independent DO loop must be a reduction statement and
therefore must be preceded by a reduce directive.

Constraint 1: The form of a reduction-stmt must be

     reduction-variable = expr reduction-op reduction-variable

or

     reduction-variable = reduction-variable reduction-op expr
or
     reduction-variable = &
       reduction-function(reduction-variable, expr)

or
     reduction-variable = &
       reduction-function(expr, reduction-variable)

A variable that appears on the left hand side of a reduction statement
is said to be a reduction variable.  The operator or function that
precedes or follows the reduction variable reference after the = sign
is the reduction operator.

H498.1  reduction-variable    is  variable

H498.2  reduction-op          is  intrinsic-op 

H498.3  reduction-function    is  MAX
                              or  MIN
                              or  IAND
                              or  IOR
                              or  IEOR
                              or  MATMUL

Constraint 2:  The two reduction-variable references that occur in a
reduction statement must be lexically identical.

Constraint 3:  A reduction-op may not be a rel-op; and it may not be
power-op (i.e. **).

Constraint 4:  The following two forms of the reduction-stmt are not allowed:

    reduction-variable = expr - reduction-variable

    reduction-variable = expr / reduction-variable

Constraint 5:  A reduction-stmt must occur in the range of an
INDEPENDENT DO loop in which the reduction variable does not have the
NEW attribute.    Additional constraints on the location of a reduction
statement are covered in Section 2.5.

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 3).  All are commutative except for
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:   A - B means A + (-B).)

According to Constraint 5 above, a reduction statement occurs in an
INDEPENDENT loop in which the reduction variable is not NEW.  In fact,
it may occur in a nest of such loops.   It is said to be *protected* in
the outermost of these.  While this loop or loop nest is active, the
reduction variable may occur only as a reduction variable; it does not
have a defined value until control passes to the statement following
the loop.   This issue is discussed fully in Section 2.5.

2.3 Relaxation of the HPF Standard's requirement for 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) unless all assignments in the loop to
the object occur in reduction statements.  (end added text)

Thus, the user's assertion implied by INDEPENDENT about the behavior of
the program is weakened, and no longer applies to reduction variables
that are protected in the DO construct.  On exit from the outermost
loop in which it is protected the value of a reduction variable is
processor dependent, although it is constrained by the model
implementation mechanism given in Section 3.

2.4  Example.

Here is an ADD_SCATTER done as an independent loop:

    !hpf$ independent
    do i = 1, n
        !hpf$ reduce
        X(index(i)) = X(index(i)) - f(i)
    enddo

The implementation will most likely make a local copy of an array of the
same shape as X, and initialize it to zero.  Each iteration will
subtract the value of f(i) from its local X(index(i)).  To
create the final result, it will combine all the local arrays with the
initial value of X:  the combining operator is addition, so that the
result is the sum of the initial value of X and the local arrays.

2.5  Restrictions on the Lexical Location of Reduction Statements, and on the
Use of Reduction Variables

In this discussion, we use the term *range* (of a DO construct) and
*active* (an active DO construct) in their technical sense as defined
by Fortran 90.  The range of a DO construct consists of the statements
lexically bounded by the DO statement and the terminal statement of the
construct.  Thus, a subprogram invoked from inside a DO loop is not in
its range, although the CALL statement or function reference that
invokes it is.    On the other hand, the DO construct is active from
the time the 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 INDEPENDENT loops (see page 302 of the
Handbook and Section 4.4 of the HPF Standard.)

When a reduction operation is executed, some nest of DO constructs will
be active; this is guaranteed by the Constraint 5 in Section 2.2.  The
reduction variable is said to be protected in some, and possibly all of
the active nest; the rule is this:  it is protected in the outermost DO
construct that is INDEPENDENT, has no NEW clause that names the
reduction variable, and does not *contain* an independent directive in
which the reduction variable occurs in a NEW clause (see Example 2,
below).   It is protected in all of this DO construct, including any
other control structures nested within it.

HPF requires that any reduction statement occur in the range of the 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.

2.5.1 Example Loop Nests

When loops are nested, a reduction variable may 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.

!   Nested Loop Example 1.  Inner loop is sequential

      x = 10
      outer: do while (x < 1000)
         !hpf$ independent, new(j)
         middle: do i = 1, n
	   inner: do j = 1, m
             !hpf$ reduce
	     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

Since the variable x occurs in a reduction statement in the range of
the independent middle loop, it is a protected reduction variable
throughout that loop, including inside the inner loop.   The outermost
loop is not INDEPENDENT, and so x is not protected in that part of its
range not included in the middle loop.

This remains true if the inner loop is independent.   On the other
hand, a variable that occurs in a 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:

!   Nested Loop Example 2.  Outer loop NEW clause.

      !hpf$ independent
      outer: do k = 1, 100
         !hpf$ independent, new (x)
         middle: do i = 1, n
            x = 10
            !hpf$ independent
            inner: do j = 1, m
              !hpf$ reduce
              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

Here, 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.

	program example
        integer x, xinc, i
	
	x = 10
        !hpf$ independent, new(xinc)
	do i = 1, n
            call sub(i, xinc)
            !hpf$ reduce
	    x = x + xinc
	enddo

	subroutine sub(i, xinc)
	integer i, xinc, local_b
	  local_b = (i-1)**2
          xinc = i * local_b
	end subroutine sub

2.6  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.

3.  Implementation and Semantics  

The result of the Fortran 90 intrinsic function  SUM  is defined to
be a processor dependent approximation to the sum of the elements of
its argument array.  On exit from the outermost 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 outside of
a reduction statement in the loop range, it is not necessary to define
the values that these variables may have inside the loop.

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.  Each processor allocates a local
reduction variable 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:

	Operator	Identity Element
	-------		----------------
	+		0
	-		0
	*		1
	/		1
	.and.		.true.
	.or. 		.false.
	.eqv.		.true.
	.neqv.          .false.

	//		''

        Table 1.   Identity elements for intrinsic reduction operators.

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

	Operator	Identity Element
	-------		----------------
	IAND(I,J)	(all ones)
	IOR(I,J)        0
	IEOR(I,J)       0
	MIN(X,Y)	-HUGE(X) where x has the same type and kind type
			parameter as the reduction variable
	MAX(X,Y)	+HUGE(X) where x has the same type and kind type
                        parameter as the reduction variable
        MATMUL(A,B)     MATMUL is noncommutative; it is the only allowed
                        associative and noncommutative intrinsic reduction function.
			As a consequence of the form A = MATMUL(A, B)
			of the reduction operation, it must be
			true that SIZE(A,2) = SIZE(B,1) = SIZE(B,2).
			As a consequence of the alternate form A = MATMUL(B, A)
			of the reduction operation, it must also be
			true that SIZE(A,1) = SIZE(B,2) = SIZE(B,1).  Thus,
			the identity is the identity matrix of shape (/N, N/),
                        of the same type and kind type parameters as the reduction
                        variable; for the left identity, N = SIZE(A,1), for
			the right identity, N = SIZE(A,2).

        Table 2.   Identity elements for intrinsic reduction functions.

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

As an example, consider:

        integer, external :: intfun
	double precision y(100)
	real x(100)

  	!hpf$ independent, new(y)
	do i = 1, n
           y = ...
           !hpf$ reduce
	   x = intfun(i) + x
           !hpf$ reduce
	   x = x - y
        enddo

The local reduction variable will be a real array of shape (/ 100 /).
Thus, the double precision  values (of L + 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 *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 (SUM, etc.)

Next, we give an implementation for the noncommutative operator 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 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 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

(1)   lrv(P) <op> lrv(P-1) <op> ... <op> lrv(1) 
                     <op> A <op> 
                            rrv(1) <op> ... <op> rrv(P)

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.

4.   Examples.

The reduction variable may be a scalar array element:

      !hpf$ independent
      do i = 1, n
          !hpf$ reduce
          x(indx(i)) = x(indx(i)) + <expr>
      enddo

It may be to a section of an array:

      !hpf$ independent
      do i = 1, n
          !hpf$ reduce
          x(i: i+4) = x(i: i+4) + <expr>   ! as many as 5 updates to elements of x
      enddo

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.

5.   Alternative syntax proposed.

In the proposal above, the programmer explicitly identifies the reduction statements,
and only implicitly the loop in which each reduction variable is protected.
The alternative reverses these decisions.  Thus, we have:

      x = 10
      outer: do while (x < 1000)
         !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

Now, the loop in which a reduction variable is protected is labeled explicitly.
Each and every occurrence of X encountered while the loop is active must be
as a reduction variable in a reduction statement of the form specified in Section
2.2 above, that occurs in the lexical range of the loop.

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).

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

6.1  COMMUTATIVE attribute for functions.

First, allow COMMUTATIVE as an addition attribute of functions:

        PURE, COMMUTATIVE FUNCTION FOOBAR(A,B)

The user asserts that the value of FOOBAR(A,B) is equal to the value
of FOOBAR(B,A) for all inputs A and B.

6.2  Initial value for local variables in reduction clause

The reduction clause is allowed to designate an initial value
for the local reduction variable:

!hpf$ independent, new(B), reduction (A = MY_OP_IDENT(SIZE(B,1)))
do i = 1, n
   B = ...
   A = MY_OP(A, B)
enddo

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 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.

---------------------------------------------------------------------------
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>
---------------------------------------------------------------------------