Subject: Independent do loops

From: [email-masked] (Geoffrey Fox)
Date: Fri, 21 Aug 92 12:55:16 EDT

"INDEPENDENT" DO LOOPS
This memo has evolved from discussions with Chuck Koebel. It describes two classes of
"independent" DO loops which are not easy in current HPF. The first class is actually
hardest as has reductions. The second solely is independent calculations. 

This note describes a class of problems which is probably to date the largest 
academic use of parallel computers over the last 15 years and further is the 
class for which systems such as PVM, Linda and "network-express" are currently
attracting the greatest attention in industry.  My understanding is that HP 
Fortran will not support them.  This surprises me.

The simplest example is the Monte Carlo integral.  We have a set of points 
x(i).  We need to calculate function values f(x(i)) and global sums such as 
sum(i) f(x(i)) and f(x(i))**2.  The parallel implementation is straightforward.
The f(x(i)) are calculated independently in separate processors and the global
averages by some suitable combine operation.  The original (high performance) 
Fortran code would be:

	av = avsq = cnt = o
	do i = 1, forever
	get x = x(i) (from tape or random number)
	calculate f = f(x)
	av = av + f
	avsq = avsq + f * f
	cnt = cnt + 1
	end do

        You can argue that this is trivial in HPF by declaring x f av avsq to be
        vectors. Unfortunately this is not possible in many real examples where
        these variables generalize to data structures with values of "forever" such
        that vectors would take up terabytes of memory. In natural MIMD implementation
        one would only be aware of N copies of x f av avsq at a time where N is number 
        of processors.

Real (and hence more complicated) examples would be
1)	Calculate string theory of gravity
	x(i) represent mesh of points on a two dimensional surface
        which is generated randomly initially and evolved independently.
	f(x(i)) are quantum field values on surface
	av (and error avsq) are integrals of "action."
	One plots av as a function of coupling constant to look for phase 
transitions:

2)	Analyze data from accelerators.

This application has, for last 15 years, used MIMD "processor farms" in 
production.  Ferimilab has been leader recently.

Here x(i) are read from tape and specify an interaction with from 2 to many 
hundreds of particles.  There are a huge number of f(x(i)) specifying feature 
of the physics and response of apparatus.  The averages are now accumulated 
more sophisticatedly in histograms and scatterplots.

In this case, the code inside the DO i=1, forever could well be 100,000 lines 
long.  The size of histogram and scatterplot storage needed may or may not fit
on a single node.  The variation in processing time from event to event is 
many orders of magnitude.

3)	In 1) and 2) one has a simple source of "events" x(i).  In particle 
shower or cascade calculations, the independent x(i) are not naturally 
specified by a DO LOOP.  Examples are the calculation of the response of a 
uranium plus scintillator detector to an incident photon which produces a 
cascade shower; the rays in a classic ray tracing graphics rendering  package 
(these may or may not spawn additional rays in their evolution through image);
alternatively, we have multiplication of neurons meandering through a doomsday
world.  
  Here the x(i) are particle (ray) parameters and f(x(i)) the result of 
tracking particle through medium.  However, the result of the tracking is one 
or more additional particles which must recursively be tracked.

In 1) and 2), we have a rather simple master-slave formalism where master just 
doles out x(i) from loop iterations.  In 3), the slaves will produce the x(i) 
which are either handed back to the master or distributed by the slave itself.

Other implementation issues include parallel random-numbers; distributed 
histograms; distributed data sets to specify the world in 3).  Initially, one 
could assume that each node had enough memory to store a full set of 
histograms.  One would accumulate independently into these during the loop 
over x(i) and then combine histograms at the end.

Language Issues:

Neither the FORALL construct nor the INDEPENDENT assertion as they
currently appear in HPF are suitable for these problems for the following
reasons.

1) FORALL - This is an array assignment, and no data arrays are stored. 
Rewriting the code in vector format would be very difficult, if not
impossible, due to the complexity of the f()'s used.  Even if this were
possible, the sizes of the arrays needed would be prohibitive.  Finally,
there is still the problem of the data reductions, not all of which can be
done directly with HPF intrinsics.

2) INDEPENDENT - This asserts that no memory location is written by one
iteration and accessed by another.  This is not the case with the original
Fortran, where all of the variables are reused from iteration to iteration.
 Even with a "private" declaration, the data reductions by their nature
must both read and write the same locations on many iterations.

Actually there appeaers to be a problem with an even simpler class of applications
which are typified by
a)Chemical Calculations where you first find independently matrix elements
and then manipulate matrices.
b)Circuit Simulations(spice, electrical grid, gas supply, phone etc.)
First independently calculate very sparse matrix elements (e.g. model components)
then solve sparse matrix.

In both a),b) one can see serious load balance issues as different matrix elements
have different calculational complexities.
The INDEPENDENT component is
Do I=1,numrow
do j=1,nonzero elements in column
Long complicated set of calculations with "private" reused variables (different values on
each node of parallel machine)
M(i,j)=answer
enddo
 We need again to declare variables used inside loop as "private"(replicated variables with same
name but different value on each node of MIMD parallel machine). It is perhaps controversial but
my experience would be that vast majority of all nondecomposed arrays will be best thought of
as Private not global. i.e. private should be default for all variables not appearing in a
decomposition statement.