Subject: User-defined elemental procedures

From: John Merlin <[email-masked]>
Date: Sat, 27 Jun 92 16:58:46 BST

I'd like to propose extending the Fortran 90 concept of 'elemental 
procedures' for the purpose of HP Fortran -- the extension being to 
allow programmers to define such procedures, which in F90 are restricted 
to a subset of the intrinsic procedures.  I believe there are several 
grounds for introducing such procedures in HPF (and indeed in F90 
generally), namely: enhanced expressiveness and elegance; efficiency; 
and the possibility of expressing a limited form of MIMD parallelism 
in a simple and elegant way.

(Incidentally, I've heard on the grapevine, and in one of Chuck's 
messages, that this proposal was made at the last HPF meeting, 
so this is probably more of a case of 'seconding' than 'proposing'.
However, since I don't know anything about that proposal I still have 
an excuse for giving my two-pennies worth:-)

In this message I'll outline the proposal, explain the constraints,
and summarise what I see as the uses and advantages of this feature.

Introduction
------------

Fortran 90 introduces the concept of 'elemental procedures',
which are defined for scalar arguments but may also be applied to
conforming array-valued arguments.  For an elemental function,
each element of the result, if any, is as would have been obtained by
applying the function to corresponding elements of the arguments.
Examples are the mathematical intrinsics, e.g SIN(X).

Unfortunately, Fortran 90 restricts the application of elemental
procedures to a subset of the intrinsic procedures---the programmer
cannot write his own.  This note therefore proposes the extension of 
allowing the programmer to define elemental procedures.

Informal proposal
-----------------
To define an elemental procedure, the subroutine or function statement 
is prefixed with the new keyword 'ELEMENTAL', e.g.:

	ELEMENTAL SUBROUTINE S (X)

	ELEMENTAL REAL FUNCTION F (X)

Henceforth, most of the discussion is in terms of elemental functions;
the properties of elemental subroutines follow in an obvious fashion.

For a function to be valid for use as an elemental function, it must 
obey the following constraints:

	1. Its arguments must have intent '(IN)'.
	2. Local arrays cannot be SAVEd.
	3. It cannot define any global variable.
	4. Its arguments, result and local variables must not be distributed.
	5. It can only reference a global variable if it is *not distributed*.

('local' and 'global' are used here in the normal programming sense, 
rather than referring to distribution!).

Incidentally, I've chosen not to insist that the dummy arguments and 
result are scalar---that's another extension to the F90 definition!  
Rather, each dummy argument and the result can be an array of
constant shape, and their shapes need not be related.  An elemental
invocation of the function corresponds to passing it actual arguments
conforming with the corresponding dummy arguments, in the normal way.
However, the function can also be invoked such that each actual argument
is an array of shape (s1,...sn) of the corresponding dummy argument
-- where (s1,...sn) is the same for all arguments.  This is equivalent 
to an array of elemental function invocations, which are conceptually 
computed in parallel, and whose result is an array of shape (s1,...sn) 
of 'elemental' results, each of which is the same as would have been 
obtained by applying the function to the corresponding 'elements' of 
the arguments.   (This generalisation, if acceptable, would mean that 
the proposed functions are a combination of 'elemental' and 
'transformational' functions in the F90 terminology.  I don't foresee 
any problems with this generalisation of the F90 concept -- except for 
one restriction on its use -- but I may be wrong!  I'll say why I think 
this extension may be useful, and mention the restriction, below).

This can be expressed more formally as follows:

	In an invocation of an user-defined elemental function,
	the shape (ai_1,ai_2,...ai_ni) of each actual argument Ai 
	must be related to the shape (di_1,di_2,...di_mi) of the 
	corresponding dummy argument Di as follows:

	For all arguments i:
		ni >= mi
		ai_k = di_k,   k = [1, mi]

	(i.e. the lowest dimensions of each actual argument must conform 
	with the corresponding dummy argument).

	For all argument pairs i,j:
		(ni - mi) = (nj - mj)
		ai_k = aj_k',    k > mi,  k' > mj,  (k - mi) = (k' - mj). 

	(i.e. the 'extra' dimensions of all arguments must conform).

	The shape of the actual function result, (f_1,f_2,...f_n), is
	then related to that of the dummy function result, (r_1,r_2,...r_m)
	as follows:
		(n - m) = (ni - mi)
		f_k = r_k,    k = [1, m]
		f_k = ai_k,   k > m
	
	(i.e. the lowest dimensions of the actual result conform with the 
	dummy result, and the 'extra' dimensions of the result conform with
	the 'extra' dimensions of each argument).

Comments
--------
The constraints are obviously designed so that the function can
be invoked concurrently at each of a set of grid points.  (By 'grid 
points' I really mean elements of an underlying array, but with with the 
understanding that the data objects at each grid point can themselves 
be arrays).

Constraints 1-3 ensure that the function has no side effects.
(Actually not quite -- I should also add the rule that the function 
performs no I/O!).  The last two concern the use of the function in 
an environment with distributed data, and ensure that it does not 
perform any data communications.  Since by definition an elemental 
function is invoked independently in different processes (and may be
invoked in some processes but not others), this constraint is essential, 
for otherwise the function could not be used safely in environments 
where both the sending and receiving processes must participate in a 
communication (e.g. most message-passing machines).

These rules show that a compiler cannot deduce from a procedure's 
interface whether it can validly be used as an elemental procedure, 
as that depends on its local declarations and internal operations.  
Hence, it is necessary to use a specifier like 'ELEMENTAL' in the 
procedure interface to identify such procedures.  The compiler can 
check that the procedure satisfies all the necessary constraints 
when it compiles the procedure itself.

Uses and MIMD aspects
---------------------

1. Array expressions

Obviously such functions can be used in array expressions with the 
same interpretation as F90 elemental intrinsic functions.
The function result must be scalar or conform with the expression.
In the former case, the result is effectively 'broadcast' to an array
of the required shape (although in the HPF context no data communications 
are required, as the function would be invoked identically in all 
processes where the result is required).  In the latter case, it's as 
if the function is invoked elementally at each element of the 'underlying' 
array (i.e. the array left after omitting the first 'm' dimensions, 
where 'm' is the rank of the result).

I've allowed the generalisation of array-valued dummy arguments and 
result to handle applications that are data-parallel over an underlying 
grid, but where the data objects at each grid point are arrays rather 
than scalars.  An example is QCD (a physics application), which is
defined on a 4d grid, and where the data objects are a vector of
length 3 at each grid point and a 3*3 array on each link.  The vectors
and matrices can be initialised independently at each grid point
(viz. with an elemental function), but each vector and matrix must be
initialised as a whole as their individual elements are not independent
(e.g. the vectors must have length 1).

Incidentally, I should mention that elemental procedures can be used 
to perform an operation at a grid point even if the operation is a 
function of values at *other* grid points;  this can be done in parallel 
as long as the old values at the other sites are required.  In the 
implementation I envisage, all data communications are done *outside* 
the elemental procedure prior to the call, and the required values are 
received by the procedure locally via its argument list.  E.g. in QCD,
the updating of the matrix on a link is a function of the (old values of)
the matrices on the 3 neighbouring links forming a 'staple', e.g:

                        Ub
                     --------
                     |      |
                  Ua |      | Uc
                     |      |
                     ........
                      U_new

This can be computed using an elemental function, e.g.:

	U_new = link_update (Ua, Ub, Uc)

where the arguments and result are arrays of shape (3,3).

It can be done in parallel over the whole grid in an array assignment.
(Here I update the links in direction 1 of a 2d plane for simplicity.
'U1' and 'U2' are arrays containing the link matrices in directions
1 and 2 respectively.  They have dimensions (3,3,n1,n2) -- the first 
two dimensions being the link matrix and the last two the lattice indices):

	U1 = link_update (U2,                  &    ! Ua
	                  CSHIFT (U1, 1, 3),   &    ! Ub--shifted in dim 3
	                  CSHIFT (U2, 1, 4))        ! Uc--shifted in dim 4

The data communications are done (by 'CSHIFT()') prior to calling 
the elemental function.

2. WHERE statement & construct

Because of its elemental property, an elemental function can also be used 
in masked array assignments in WHERE statement and constructs (cf. normal 
array-valued functions cannot).

Constraint:  An elemental function used in a WHERE statement or construct
must have a scalar result.  (This is because the 'mask' array can select 
individual elements of the assignment for execution.  I don't think that 
any restriction on the shape of the dummy arguments is necessary though).

In the HPF context, this may be one way of obtaining MIMD parallelism, e.g.:

	REAL x (10,10)
	LOGICAL edges (10,10)
	INTERFACE
	  ELEMENTAL REAL FUNCTION f_egde (x)
	     REAL x
	  END FUNCTION f_edge
	  ELEMENTAL REAL FUNCTION f_interior (x)
	    REAL x
	  END FUNCTION f_interior
	END INTERFACE

	...  initialise mask array 'edges'

	WHERE (edges)
	  x = f_egde (x)
	ELSE WHERE
	  x = f_interior (x)
	END WHERE

(Of course, MIMD parallelism is only obtained if the compiler
can establish that the two assignments are independent and doesn't
force a synchronisation at the ELSEWHERE, which is probably a reasonable 
assumption in simple cases).

3. FORALL statement & construct

These functions can also be used in a FORALL.  Because a 'forall-assignment'
may be an 'array-assignment' (in most definitions anyway) the elemental
function can have an array result.  E.g. if a certain problem
is data-parallel over a 2d grid, and the data structure at each grid
point is a vector of length 3 (2d QCD?), we could have:

	REAL  v (3,10,10)
	INTERFACE
	  ELEMENTAL FUNCTION f (x)
	    REAL, DIMENSION(3) :: f, x
	  END FUNCTION f
	END INTERFACE
	...
	FORALL (i=1:10, j=1:10)  v(:,i,j) = f (v(:,i,j)) 

If 'IF' constructs are allowed in FORALL (a thorny subject at the moment!)
one has more opportunities for MIMD parallelism, e.g.:

	FORALL (i=1:n, j=1:n)
	  IF (i==1 .OR. i==n .OR. j==1 .OR. j==n)
	    x (i,j) = f_edge (x(i,j))
	  ELSE
	    x (i,j) = f_interior (x(i,j))
	  ENDIF
	END FORALL

(Incidentally, this can also be coded without an IF construct, using the 
F90 elemental intrinsic function MERGE(), thus:

	FORALL (i=1:n, j=1:n)                                        &
	    x (i,j) =  MERGE (f_edges (x(i,j)),                      &
	                      f_interior (x(i,j)),                   &
	                      (i==1 .OR. i==n .OR. j==1 .OR. j==n))
)

I would propose that only elemental functions (both user-defined and
intrinsic) are allowed in FORALL, just as in array expressions.  
Elemental subroutine calls could also be allowed in FORALL, with a very 
similar interpretation to elemental functions.  The main reason for using 
them would be to allow the return of multiple results.  
(Incidentally, arguments returning a result can have intent '(IN)' or 
'(INOUT)' as in F90).

4. MIMD parallelism

Perhaps the most direct way of obtaining MIMD parallelism is by
means of branches within an elemental function which depend on argument 
values.  These branches can be governed by content-based or index-based 
conditionals (the latter in a FORALL context).  For example:

	ELEMENTAL FUNCTION f (x, i)
	  IF (x > 0) THEN     ! content-based conditional
	    ...
	  ELSE IF (i==1 .OR. i==n) THEN    ! index-based conditional
	    ...
	  ENDIF
	END FUNCTION

	...
	FORALL (i=1:n)  x (i) = f(x(i), i)
	...

Content-based conditionals can be exploited generally, including in
array assignments, which may sometimes obviate the need for WHERE 
statements and constructs with their potential synchronisation overhead. 

Other advantages
----------------

There are other advantages to user-defined elemental procedures, apart 
from their MIMD potential and the ability to use them in WHERE and FORALL.

They would be a very convenient programming tool, as the same elemental 
procedure can be applied to arguments of any rank.

In addition, the implementation of an elemental function array in an
array expression is likely to be more efficient than that of an 
equivalent array function.  One reason is that it requires less temporary
storage for the result (i.e. storage for a single result versus storage 
for the entire array of results).  Another is that it saves on
looping if an array expression is implemented by sequential iteration
over the component elemental expressions (as may be done for the 'segment'
of the array expression local to each process).  This is because,
in the sequential version, the elemental function can be invoked
elementally in situ within the expression.  The array function,
on the other hand, must be executed before the expression is evaluated, 
storing its result in a temporary array for use within the expression.
Looping is then required during the execution of the array function
body as well as the expression evaluation.

                        John Merlin.

(P.S. I'm away for the next two weeks so I won't be able to answer your 
comments till mid July.  Sorry about that!)
-----------------------------------------------------------------------
John H. Merlin                               email: [email-masked]
Dept. of Electronics and Computer Science,   tel:   +44 703 593368
University of Southampton,                   fax:   +44 703 593045
Southampton S09 5NH,  U.K.

---

From: [email-masked] (Joel Williamson)
Date: Mon, 29 Jun 92 7:49:06 CDT

John Merlin writes:

	...stuff deleted...
> Informal proposal
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
> rather than referring to distribution!).

It seems that the above constraints preclude the possibility of
elemental subroutines unless constraint 1 is relaxed.
	...stuff deleted...
> '(INOUT)' as in F90).                                         ^^
                                                               OUT? 

This seems to be the needed relaxation, but allowing modification of
subroutine arguments probably opens yet another proverbial can of worms,
since a pathological and insidious programmer can find ways to introduce
innumerable side effects via this mechanism.  I think we'll be far
better off restricting ourselves to elemental functions only.

In any case, if we choose to allow function references within FORALL,
this proposal appears to define it in a clean, implementable way.

Best regards,

Joel Williamson