A Linear Algebra Framework for Static HPF Code Distribution

Corinne Ancourt  Fabien Coelho
François Irigoin  Ronan Keryell

Centre de Recherche en Informatique, Ecole Nationale Supérieure des Mines de Paris^{∗} 35, rue Saint Honoré, 77305 Fontainebleau cedex, France

Abstract

High Performance Fortran (hpf) was developed to support data parallel programming for simd and mimd machines with distributed memory. The programmer is provided a familiar uniform logic address space and specifies the data distribution by directives. The compiler then exploits these directives to allocate arrays in the local memories, to assign computations to elementary processors and to migrate data between processors when required. We show here that linear algebra is a convenient framework to encode hpf directives and to synthesize distributed code with space-efficient array allocation, tight loop bounds and minimal communications for INDEPENDENT loops. The systematic use of an affine framework makes it easy to prove the compilation scheme correct. The generated code includes traditional optimizations such as guard elimination, message vectorization, and aggregation, overlap analysis,...

1 Introduction

Distributed memory multiprocessors can be used efficiently if each local memory contains the right pieces of data, if local computations use local data and if missing pieces of data are quickly moved between processors.

Macro packages and libraries are available to ease the programmer’s burden but the level of details still required transforms the simplest algorithm, e.g. a matrix multiply, into hundreds of lines of code. This fact decreases programmer productivity and jeopardizes portability, as well as the economical survival of distributed memory machines.

Manufacturers and research laboratories, lead by Digital and Rice University, decided in 1991 to shift part of the burden onto compilers by providing the programmer a uniform address space to allocate objects and a (mainly) implicit way to express parallelism. Numerous research projects [15, 19, 35] and a few commercial products had shown that this goal could be achieved and the High Performance Fortran Forum was set up to select the most useful functionalities and to standardize the syntax. The definition of the new language, hpf, was frozen in May 1993 [13] and the first compilers are already available [7, 8, 11, 35, 24]. This quick delivery was made possible by defining a subset of the language which only allows static distribution of data and prohibits dynamic redistribution.

This paper deals with this hpf static subset and shows how simple changes of basis and affine constraints can be used to relate the global memory and computation spaces seen by the programmer to the local memory and computation spaces allocated to each elementary processor. These relations, which depend on hpf directives added by the programmer, are used to allocate local parts of global arrays and temporary copies which are necessary when non-local data is used by local computations. These constraints are also used in combination with the owner-compute rule to decide which computations are local to a processor, and to derive loop bounds. Finally they are used to generate send and receive statements required to access non-local data.

These three steps, local memory allocation, local iteration enumeration and data communication, are put together as a general compilation scheme for parallel loops, known as INDEPENDENT in hpf, with affine bounds and subscript expressions. This compilation scheme directly generates optimized code which includes techniques such as guard elimination [15], message vectorization and aggregation [35], and overlap analysis [15]. There are no restrictions neither on the kind of distribution, block or cyclic, nor on the rank of array references.

The memory allocation part is independent of parallel loops and can always be used. The relations between the global programmer space and the local processor spaces can also be used to translate sequential loops with a run-time resolution mechanism or with some optimizations.

The reader is assumed knowledgeable in hpf directives [13] and optimization techniques for hpf [15, 35].

The paper is organized as follows. Section 2 shows how hpf directives can be expressed as affine constraints and normalized to simplify the compilation process. Section 3 presents an overview of the compilation scheme and introduces the basic sets own, compute, send and receive that are used to allocate local parts of hpf arrays and temporaries, to enumerate local iteration and to generate data exchanges between processors. Section 4 refines these sets to minimize the amount of memory space allocated, to reduce the number of loops whenever possible and to improve the communication pattern. This is achieved by using different coordinates to enumerate the same sets. Examples are shown in Section 5 and the results are compared with previous work in Section 6.

∗{ancourt,coelho,irigoin,keryell}@cri.ensmp.fr.

2 HPF directives

In hpf, data distribution is specified in two steps:

1. array elements are mapped onto virtual processor arrays called templates by an ALIGN directive;

2. templates, and implicitly all array elements mapped on them, are distributed onto hpf processors.

Array elements aligned on the same template are allocated on the same elementary processor. Expressions using these elements can be evaluated locally, without inter-processor communications. This alignment step mainly depends on the algorithm.

These templates are packed in blocks to reduce communication and scheduling overheads without increasing load imbalance too much. The block sizes depend on the target machine, while the load imbalance stems from the algorithm.

Templates can be bypassed by aligning an array on another array, and by distributing array directly on processors. This does not increase the expressive power of the language but implies additional check on hpf declarations. Templates are systematically used in this paper to simplify algorithm descriptions.

Our framework deals with both stages and could easily tackle direct alignment and direct distribution. The next sections show that any hpf directive can be expressed as a set of affine constraints.

2.1 Alignment

The ALIGN directive is used to specify the alignment relation between one object and one or many other objects through an affine expression. To avoid denoting constant terms in affine expressions, homogeneous coordinates are used when necessary. The cartesian domains of Fortran array and hpf template declarations are expressed by linear homogeneous^{1} inequalities such as

D̃a ≤ 0, T̃t ≤ 0  (1)

where a and t describe the array and template elements respectively.

The alignment of an array A on a template T is an affine mapping from a to t. Alignments are specified dimension-wise with integer affine expressions as template subscript expressions. each time only one index of the template and once is used in an affine expression it can no longer be used in another expression of the same alignment [13]. Each array index can be used at most once in a subscript expression in any given alignment [13], and each subscript expression cannot contain more than one index. Thus the relation between A and T is given by

t = Ãa  (2)

where Ã is a matrix with at most one non-zero integer element per column and per row (without considering the last homogeneous column).

Equation (2) cannot express replication of elements of A on dimensions of T although it is allowed by hpf. A projection matrix R̃ is added to leave unspecified template dimensions unconstrained and Equation (2) becomes:

R̃t = Ãa  (3)

The same formalism can also deal with a chain of alignment relations (A is aligned on B which is aligned on C,...) as allowed in hpf.

The relation T̃t ≤ 0 is redundant with this last equation and Fortran array declaration inequalities because, according to hpf semantics, an array element cannot be aligned outside a template.

These notations are illustrated with four examples in a seismic 2d benchmark translated from Fortran 77 to hpf. The declaration

                                                   real V(3:200,2:100)
                                             !HPF$ template T(3:200,2:100)
                                             !HPF$ align V(I,J) with T(I,J)

means a perfect alignment translated in

D̃ = T̃ = ( \begin{matrix}−1 & 0 & 3 \\ 1 & 0 & −200 \\ 0 & −1 & 2 \\ 0 & 1 & −100 \end{matrix} )

R̃ = Ã = ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix} )

with

1Homogeneous matrices are labelled with “ ”.

t = ( \begin{matrix}t_{1} \\ t_{2} \\ 1 \end{matrix} ) , a = ( \begin{matrix}a_{1} \\ a_{2} \\ 1 \end{matrix} )

in the homogeneous coordinates. The second declaration

real UINT_W(10:40)
!HPF$ align UINT_W(I) with T(5*I-7,2)

is represented by

R̃ = ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix} ) , Ã = ( \begin{matrix}5 & −7 \\ 0 & 2 \end{matrix} ) , D̃ = ( \begin{matrix}−1 & 10 \\ 1 & −40 \end{matrix} )

The third one

real U(198,0:197,3)
!HPF$ align U(I,J,*) with T(I,J)

leads to

\begin{aligned}R̃ & = & ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix} ) , & Ã & = & ( \begin{matrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{matrix} ) \\ D̃ & = & ( \begin{matrix}−1 & 0 & 0 & 1 \\ 1 & 0 & 0 & −198 \\ 0 & −1 & 0 & 0 \\ 0 & 1 & 0 & −197 \\ 0 & 0 & −1 & 1 \\ 0 & 0 & 1 & −3 \end{matrix} ) \end{aligned}

since the third dimension of the array is collapsed on the template.

Finally the last example uses replications along the lines of the template and the second dimension of the array is collapsed:

real SISMO(60,1000)
!HPF$ align SISMO(I,*) with T(*,2*I+50)

which is represented by

\begin{aligned}R̃ & = & ( 0 1 0 ) , & Ã & = & ( 2 0 50 ) \\ D̃ & = & ( \begin{matrix}−1 & 0 & 1 \\ 1 & 0 & −60 \\ 0 & −1 & 1 \\ 0 & 1 & −1000 \end{matrix} ) \end{aligned}

To summarize, “*” in an array reference induces a zero column in Ã and “*” in a template reference a zero column in R. When no replication is specified, R is the identity matrix.

2.2 Distribution

Once optional alignments are defined in hpf, objects or aligned objects can be distributed on cartesian processor sets, declared like arrays and templates, and represented by inequality:

P̃p ≤ 0  (4)

Each dimension can be distributed by block, or in a cyclic way on the processor set or even collapsed on the processors. hpf allows a different block size for each dimension. The extents (n in BLOCK(n) or CYCLIC(n)) are stored in a diagonal matrix, C.

Such a distribution may be written as a linear relation between the processor coordinate p, the template coordinate t and two additional variables, ℓ and c. Vector ℓ represents the offset within one block in one processor. Vector c represents the number of wrap-arounds that must be performed to allocate blocks cyclically on processors.

Π̃t = CP c + Cp + ℓ  (5)

The projection matrix Π̃ is needed when some dimensions are collapsed on processors. These dimensions are thus orthogonal to the processor parallelism and can be eliminated. This matrix and the vector t are the only homogeneous objects in Equation (5).

The usual modulo and integer division operators are replaced by predefined constraints on p and additional constraints on ℓ. Vector c is implictly constrained by array declarations and, redundantly, by the TEMPLATE declaration.

The diagonal matrix P describes the processors geometry and each element of its diagonal describes the number of processors in each dimension. It is related to the homogeneous matrix P̃ which declares the processors geometry and offset. Since P is used instead of P̃ in Equation (5), the lower bounds are taken into account at the level of Π̃.

Specifying BLOCK(n) in a distribution instead of CYCLIC(n) is equivalent but tells the compiler that the distribution will not wrap around and the Equation (5) can be simplified by zeroing the c coordinate for this dimension. This additional equation reduces the number of free variables and, hence, remove a loop in the generated code. CYCLIC is equivalent to CYCLIC(1) and BLOCK without parameter is equivalent to BLOCK(⌈\frac{e_{t}}{e_{p}} ⌉) where e_{t} and e_{p} are the template and processor extents in the matching dimensions. As BLOCK, the CYCLIC declaration is translated into an additional equation, ℓ = 0. The variables are bounded by:

P̃p ≤ 0, T̃t ≤ 0, 0 ≤ C^{−1}ℓ < 1

For example, the hpf declaration

!HPF processors PE(5:20,6:10)
!HPF distribute T(block,block) onto PE

with the same template as in the first example de-

real A(0:42)
!HPF$ template T(0:127)
!HPF$ processors P(0:3)
!HPF$ align A(i) with T(3*i)
!HPF$ distribute T(cyclic(4))
A(0:42:3) = ...

Figure 1: Example in Chatterjee et al. [10]

fines the matrices

\begin{aligned}P̃ & = & ( \begin{matrix}−1 & 0 & 5 \\ 1 & 0 & −20 \\ 0 & −1 & 6 \\ 0 & 1 & −10 \end{matrix} ) , & Π̃ & = & ( \begin{matrix}1 & 0 & 3 \\ 0 & 1 & 4 \end{matrix} ) \\ C & = & ( \begin{matrix}13 & 0 \\ 0 & 20 \end{matrix} ) , & P & = & ( \begin{matrix}16 & 0 \\ 0 & 5 \end{matrix} ) \end{aligned}

since ⌈\frac{200−3+1}{20−5+1}⌉ = 13 and ⌈\frac{100−2+1}{10−6+1}⌉ = 20. The constant terms of Π are chosen to exactly align the lower bounds of the processors and the templates.

In this example, the distribution is of type BLOCK in every dimension, the CP c is useless and may be discarded. But the addition of equation c = 0 is preferred because it leads to a unique algorithm for any kind of distribution.

2.3 Normalization

To clarify the exposition of our compilation scheme, hpf declarations are normalized. Every declaration of array, processor set or template is translated to have 0 as lower bound. This obviously implies no loss of generality for the method. This normalization lets us get rid of homogeneous matrices and replace them by diagonal matrices.

Another normalization is performed. If a dimension of a template is collapsed on the processors, the corresponding coordinate is useless because even if it is used to distribute the dimension of an array, this dimension will be collapsed on the processors. Thus it can be discarded from the equations. These two normalization steps simplify the inequalities used by our algorithms and allow removing the Π̃ = I matrix from Equation (5).

All the declaration of distributions, arrays, processors and templates can thus be rewritten as

0 ≤ C^{−1}ℓ < 1 , 0 ≤ D^{−1}a < 1
0 ≤ P^{−1}p < 1 , 0 ≤ T^{−1}t < 1
Rt = Aa + t_{0} , t = CP c + Cp + ℓ

where t_{0} is the coordinate in the template of the first element of the aligned array. Strict rational inequalities are in fact replaced by non-strict integer ones by multiplying them by the determinant

!HPF$ independent

forall(i, 0 ≤ L^{−1}i < 1)

        X(S_{X}i + X_{0})=f(Y(S_{Y}i + Y_{0}),Z(S_{Z}i + Z_{0}),. . .)

Figure 2: Generic Loop

Δ of the non-inverted matrix. The non-strict upper bound is Δ − 1. For instance the example given in [10] (Figure 1) is already normalized and can be translated in the following set of affine constraints

D = (43), T = (128), A = (3), t_{0} = (0), R = (1), C = (4), P = (4).

3 Overview of the compilation scheme

The compilation scheme is based on the hpf declarations, the owner-compute rule and the spmd paradigm, and deals with INDEPENDENT loops. Loop bound expressions and array subscript expressions are assumed affine. Such a parallel loop independently assigns elements of some lhs X with some rhs expression f on variables Y, Z,... (Figure 2).

The loop is normalized to start with iteration 0 and to use unit steps. The upper bounds are encoded in matrix L. Non-unit steps are transferred into the subscript matrices S_{X}, S_{Y}, S_{Z},… Since the loop is independent, S_{X} must be injective.

3.1 Own Set

The subset of X allocated on processor p according to hpf directives must be declared in the output program, up to a one-to-one mapping. Let own_{X}(p) be this subset of X and X′ the local array mapping this subset. Each element of own_{X}(p) must be mapped to one element of X′ but some elements of X′ may not be used. Finding a good mapping involves a tradeoff between the memory space usage and the access function complexity. This is studied in Section 4.3. Subset own_{X}(p) is derived from hpf declarations, expressed in an affine formalism (see Section 2), as:

\begin{matrix}own_{X}(p) = \{x | ∃t ∃c s.t. & Rt = Ax \\ ∧ & t = CP c + Cp + ℓ \\ ∧ & 0 ≤ DX^{−1}x < 1 \\ ∧ & 0 ≤ P^{−1}p < 1 \\ ∧ & 0 ≤ C^{−1}ℓ < 1\} \end{matrix}

Subset own_{X}(p) is a bounded polyhedral set. It can be mapped onto a Fortran array by projection on c and ℓ, or on any equivalent set of variables, i.e. up to a one-to-one mapping.

3.2 Compute set

Using the owner-compute rule, the set of iterations local to p, compute(p), is directly derived from the previous set, own_{X}(p):

compute(p) = {i | S_{X}i ∈ own_{X}(p)}

Note that compute(p) and own_{X}(p) are equivalent sets if S_{X} it the identity function. This is used in Section 4.3 to study the allocation of local arrays as a special case of local iterations.

Since i is linked to p, c and ℓ, local iterations can be enumerated by scanning the i polyhedral set with DO loops over p, c and ℓ or any equivalent set of variables, i.e. up to a unimodular change of basis. When X is replicated, its new values are computed on all processors having a copy. If the number of physical processors is reasonably large, this is faster than a combination of non-redundant computations followed by broadcasts.

3.3 View set

The subset of Y (or Z,...) that is used by the loop is induced by the set of local iterations:

view_{Y}(p) = {y | ∃i ∈ compute(p) s.t. y = S_{Y}i}

Note that matrix S_{Y} is not constrained and cannot always be inverted. If the intersection of this set with own_{Y}(p) is non-empty, the usual optimization is to enlarge the local copy Y′ so as to include its neighborhood, when it is not too large^{2}. This is known as overlap analysis [15].

When Y (or X, or Z,...) is referenced many times in the input loop or in the input program, these references must be clustered according to their connexity in the dependence graph [4]. Input dependencies are taken into account as well as usual ones (flow-, anti- and output-dependencies). If two references are independent, they access two distant regions and two different local copies should be allocated to reduce the total amount of memory allocated. A unique copy should be as large as the convex hull of these distant regions. If the two references are dependent, only one local copy should be allocated to avoid any consistency problem between copies. If Y is a hpf array, its local elements must be taken into account as a special reference.

The definition of view is thus altered to take into account regions [34, 20] of arrays. The equations due to subscript expression S_{Y} are replaced by a parametric polyhedral subset of Y which can be automatically computed. For instance, the set of references to Y in:

                                                 do I1 = 1, N
                                                    do I2 = 1, N
                                                       X(I1,I2) = Y(I1-1,I2) + Y(I1+1,I2)
                                                              + Y(I1,I2-1) + Y(I1,I2+1)
                                                    enddo
                                                 enddo

is summarized by the parametric region on y = (y_{1}, y_{2}):

{y | y_{1} + y_{2} ≤ i_{1} + i_{2} + 1
∧ y_{2} + i_{1} ≤ y_{1} + i_{2} + 1
∧ i_{1} + i_{2} ≤ y_{1} + y_{2} + 1
∧ y_{1} + i_{2} ≤ y_{2} + i_{1} + 1}

Subset view_{Y} is still polyhedral and array bounds can be derived by projection, up to a change of basis. If regions are not precise enough because convex hulls are used to summarize multiple references, it is possible to use additional parameters to exactly express a set of references [2]. This might be useful for red-black sor.

If the region of Y includes hpf local elements, additional equations must be added to follow the standard allocation scheme and to access all elements in a unique way from processor p at cycle c:

t′ = Ay
t′ = CP c + Cp + ℓ′

Vectors c and p are constrained as usual but vector ℓ′ is no more constrained by C. It is indirectly constrained by y and the Y region. New variables t′ and ℓ′ are added to avoid clashes with t and ℓ which appear in compute(p).

3.4 Send and Receive Sets

The set of elements that are owned by p and used by other processors and the set of elements that are used by p but owned by other processors are both intersections of previously defined sets:

send_{Y}(p) = {(p′, y) | y ∈ own_{Y}(p) ∩ view_{Y}(p′)}

receive_{Y}(p) = {(p′, y) | y ∈ view_{Y}(p) ∩ own_{Y}(p′)}

These two sets cannot always be used to derive the data exchange code. First of all p does not need to exchange data with itself. This cannot be expressed directly as an affine constraint and must be added as a guard in the output code^{3}. Secondly, replicated arrays would lead to useless communication: each owner would try to send data to each user. An additional constraint should be added to real X′((c, ℓ) ∈ view_{X}(p)),

2The neighborhood is usually considered not too large when it is bounded by a numerical constant.

3This could be avoided by exchanging data first with processors p′ such that p′ < p and then with processors such that p′ > p, using the lexicographic order. But this additional optimization is not likely to decrease the execution time much.

Y′((c, ℓ) ∈ view_{Y}(p)), Z′((c, ℓ) ∈ view_{Z}(p))

forall(U′ ∈ (Y′, Z′, . . .))

if send_{U}(p) ≠ ∅

forall((p′, ℓ, c) ∈ send_{U}(p))

if p ≠ p′

send(p′, U′(ℓ, c))

if receive_{U}(p) ≠ ∅

forall((p′, ℓ, c) ∈ receive_{U}(p))

if p ≠ p′

U′(ℓ, c) = receive(p′)

if compute(p) ≠ ∅

forall ((ℓ, c) ∈ compute(p))

X′(S_{X}′(ℓ, c)) = f(S_{Y}′(ℓ, c), S_{Z}′(ℓ, c),...)

Figure 3: The output spmd code.

restrict the set of communications to needed ones. For instance, the cost function could be the minimal Manhattan distance between p and p′ or the lexicographically minimal vector p′ − p if the interconnection network is a |p|-dimension grid. A combined cost function might even be better by taking advantage of the Manhattan distance to minimize the number of hops and of the lexicographical minimum to insure uniqueness. These problems can be cast as linear parametric problems and solved [12].

When no replication occurs, elementary data communications implied by send_{Y} and receive_{Y} can be parametrically enumerated in basis (p′, u), where u is a basis for Y, the local part of Y. send and receive are polyhedral sets and algorithms in [2] can be used. If the last component of u is allocated contiguously, vector messages can be generated.

3.5 Output SPMD code

The generic output spmd code, parametric on p the local processor, is shown in Figure 3.

Communications between processors can be executed in parallel when asynchronous send/receive calls and/or parallel hardware is available as in Transputer-based machines or in the Paragon. The guards send_{U},... can be omitted if unused hpf processors are not able to free physical processors for other useful tasks. The loop bound correctness ensures that no spurious iterations or communications occur.

The parallel loop on U′ can be exchanged with the parallel inner loop on p′ included in the forall. This other ordering makes message aggregation possible.

Note that the set of local iterations compute(p) is no longer expressed in the i frame but in the (c, ℓ) frame. Additional changes of basis or changes of coordinates must be performed to minimize the space allocated in local memory and to generate a minimal number of loops. Each change must exactly preserve integer points, whether they represent array elements or iterations. This is detailed in the next section.

4 Refinement

The pseudo-code shown in Figure 3 is still far from Fortran. Additional changes of coordinates are needed to generate proper Fortran declarations and loops.

4.1 Iterations enumeration

A general method to enumerate local iterations is described below. It is based on searching a good lattice basis to scan them in an appropriate order, and on using unimodular transformations to switch from the user visible frame to the local frame.

The subscript function S is supposed to be affine. All object declarations are supposed normalized with 0 as lower bound. The hpf array mapping is also normalized: Π = I (Section 2.3). Furthermore, the algorithm is presented when no replication is requested (R = I) to make it clearer.

Under these assumptions and according to Section 2, hpf declarations and Fortran references are represented by the following set of equations:

\begin{matrix}alignment & t = Aa + t_{0} \\ distribution & t = CP c + Cp + ℓ \\ affine reference & a = Si + a_{0} \end{matrix} \} (6)

where p is the processor vector, ℓ the local offset vector, c the cycle vector, a the accessed element and i the iteration. In order to simplify the following discussion, it is assumed that every dimension of the template has a matching array dimension which is aligned to that dimension, so a is decomposed as a_{|} for the aligned dimensions, and a_{+} for not aligned ones (corresponding to null columns of A). Let |v| denotes the number of variables of vector v. S is of rank |i|, since the iterations are independent. Under these assumptions, we have |t| = |c| = |p| = |ℓ| = |a_{|}|.

Let x be the following set of variables:

x = (p, ℓ, c, a_{+}, a_{|}, t, i)  (7)

They are needed to describe the access. Equation (6) can be rewritten as F x = f_{0}. The set of interest is the lattice F = {x|F x = f_{0}}. The goal is to find a parametric solution for each component of x, where p would remain unconstrained. This allows each processor to scan its part of F with minimal control overhead.

Let x′ be the subset of x:

x′ = (p, ℓ, c, a_{+})  (8)

x′ is a set of free variables in F. Indeed, t can be computed with the distribution equation and x′. Then a_{|} is derived with the alignment equation, if it is non-empty. Finally a_{|} and a_{+} and the affine subscript expression define i. Moreover these variables are independent because F is of rank 2|t| + |i| (F is non-degenerated), with 5|t| + |a_{+}| + |i| variables, so there are at most 3|t| + |a_{+}| = |x′| free variables.

Set F is implicitly defined but a parametric definition is needed to enumerate its elements. Parameters can be found by computing a Hermite form [32] of integer matrix F. This form associates to F three matrices, H, P and Q, such that H = P F Q. P is a permutation, H a lower triangular integer matrix and Q a unimodular change of basis. Since F is non-degenerated, no permutation is needed and P = I. By definition, H is a lower triangular matrix which can be decomposed as (H_{L} 0). Let v = Q^{−1}x be a parameter vector in a new basis. Vector v can also be decomposed in two components like H: v = (v_{0}, v′), where |v′| = |x′|. A solution x_{0} of F x = f_{0} is obtained by computing v_{0} = H_{L}^{−1}f_{0} and x_{0} = Q_{0}v_{0}. By construction, H does not constrain v′ and Q can also be decomposed like v as Q = (Q_{0} F′). The parametric solutions are x = x_{0} + F′v′ and set F can be expressed in a parametric way:

F = {x|∃v′ s.t. x = x_{0} + F′v′}  (9)

However the generated code has to be parametrized by the processor vector p. This can be achieved by another change of basis. F′ is nondegenerated since it links the free variables of v to x. Under this condition, its Hermite form H′ = P′F′Q′ is obtained without permutations and P′ = I. Let u = Q^{′−1}v′ define a new basis, then

x − x_{0} = F′v′ = H′u  (10)

The new generating system of F is based on a triangular transformation between x and u (Equation (10)). Since H′ is a lower triangular integer matrix, and processor vector p is the first set of variables of x, p is only defined by the first |p| variables of u.

This is designed to parametrize the scanning of F for each processor p_{0}. The first |p| variables of u are functions of p_{0} through the upper left |p| × |p| submatrix. The remaining part of H′ is then used to generate the elements of F that are on p_{0}, by enumerating u values.

Let’s explain why the u’s defining x in F are easy to scan with DO loops. Let Kx ≤ k be the polyhedron that constrains x. These constraints are coming from declarations, hpf directives and normalized loop bounds. They are:

0 ≤ P^{−1}p < 1 , 0 ≤ C^{−1}ℓ < 1
0 ≤ T^{−1}t < 1 , 0 ≤ D^{−1}a < 1
0 ≤ L^{−1}i < 1

Using (10) the constraints on u can be written K(x_{0} + H′u) ≤ k, that is K′u ≤ k′, where K′ = KH′ and k′ = k − Kx_{0}.

Algorithms presented in [3] can be used to generate the loop nest enumerating the local iterations. When S is of rank |a|, optimal code is generated because no projections are required. Else, the quality of the control overhead depends on the accuracy of integer projections [29], but the correctness does not.

The correctness of this enumeration scheme stems from the unimodularity which preserves the one-to-one relationship between integer points and iterations, from the independent parallelism of the loop which allows any enumeration order, and from the choice of p as first component of x, which parametrizes the solution by the processor identity. Replication is handled by a similar algorithm. The method is applied on aligned template dimensions. Then the generated spmd code is not parametrized by the processor identity on the replicated dimension, thus every such processor also scans its local data, replicated or not. If arrays are mapped on one element along a template dimension, guards have to be inserted to select the proper processors.

4.2 Symbolic resolution

The previous method can be applied in a symbolic way, when the components of the variables are not coupled and thus can be dealt with independently. This is the case when array sections are used as in [10, 18, 33]. Equations (6) then become for a given dimension:

\begin{matrix}alignment & t = αa + t_{0} \\ distribution & t = πγc + γp + ℓ \\ affine reference & a = σi + a_{0} \end{matrix} \} (11)

where π is the number of processors (a diagonal coefficient of P), and γ the block size (idem C). In order to simplify the symbolic resolution, variables a and t are eliminated. The matrix form of the system is then f_{0} = (αa_{0} + t_{0}), x = (p, ℓ, c, i) and F = ( γ 1 πγ −ασ ).

∀p_{0} ∈ [0 … π − 1], u_{1} = p_{0} do u_{2} = \frac{−b_{u}−γp_{0}+g−1}{g}, \frac{−b_{l}−γp_{0}}{g}

do u_{3} = \frac{−ωu_{2}+ \frac{πγ}{g} −1}{\frac{πγ}{g}}, \frac{β−1−ωu_{2}}{\frac{πγ}{g}}

x = H′u + x_{0}

Figure 4: Generated code

The  Hermite  normal  form  is

H = ( 1 0 0 0 ) = P F Q, with P = I and:

Q = ( \begin{matrix}0 & 1 & 0 & 0 \\ 1 & −γ & −πγ & ασ \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} )

Let g, μ and ω be such that g is gcd(πγ, ασ) and g = πγμ − ασω is the Bezout identity. The Hermite form H′ of the three rightmost columns of Q is such that x − x_{0} = H′u with:

H′ = ( \begin{matrix}1 & 0 & 0 \\ −γ & −g & 0 \\ 0 & μ & \frac{ασ}{g} \\ 0 & ω & \frac{πγ}{g} \end{matrix} ) , x_{0} = ( \begin{matrix}0 \\ αa_{0} + t_{0} \\ 0 \\ 0 \end{matrix} )

This links the three unconstrained variables u to the elements x of the original lattice F. Variables a and t can be retrieved using Equations (11).

The translation of conses traints on x to u gives a way to generate a loop nest to scan the polyhedron. Under the assumption α > 0 and σ > 0, assuming that loop bounds are rectangular^{4}, and using the constraints in K:

0 ≤ p ≤ π − 1 , 0 ≤ ℓ ≤ γ − 1
0 ≤ a ≤ s − 1 , 0 ≤ i ≤ β − 1

the constraints on x (the a one is redundant if the code is correct) are:

( \begin{matrix}−1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & −1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & −1 \\ 0 & 0 & 0 & 1 \end{matrix} ) ( \begin{matrix}p \\ ℓ \\ c \\ i \end{matrix} ) ≤ ( \begin{matrix}0 \\ π − 1 \\ 0 \\ γ − 1 \\ 0 \\ β − 1 \end{matrix} )

and can be translated as constraints on u:

( \begin{matrix}−1 & 0 & 0 \\ 1 & 0 & 0 \\ γ & g & 0 \\ −γ & −g & 0 \\ 0 & −ω & −\frac{πγ}{g} \\ 0 & ω & \frac{πγ}{g} \end{matrix} ) ( \begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix} ) ≤ ( \begin{matrix}0 \\ π − 1 \\ −b_{l} \\ b_{u} \\ 0 \\ β − 1 \end{matrix} )

[figure: c; ℓ; u_{2}; u_{3}; u′_{2}; u′_{3}]

Figure 5: Packing of elements

where b_{l} = −αa_{0} − t_{0}, b_{u} = γ − 1 + b_{l}.

The resulting generic spmd code for an array section is shown in Figure 4. As expected, the loop nest is parametrized by p_{0}, the processor identity. Integer divisions with positive remainders are used in the loop bounds.

4.3 HPF array allocation

The previous two sections can be used to allocate parts of hpf arrays. If the subscript function S is the identity function, local iterations and local elements are strictly equivalent. Thus the constraints on local iterations can be reused as constraints on local elements. However, Fortran array declarations are not as general as Fortran loop bounds, and it is interesting to add another change of frame to fit Fortran specifications for array declarations and to reduce the amount of allocated memory.

The local array elements are packed onto local memories dimension by dimension. The geometric intuition of the packing scheme for one dimension, in the example in [10] displayed Figure 1, is shown in Figure 5 for the first processor of the dimension. The first grid is the part of the template local to the processor, in the cycle and offset coordinate system. The second grid shows the same grid through unimodular transformation from x to u. The third one is the final packed form.

The components in x′ (Definition (8)) related to the current dimension are considered. If the coordinate is in a_{+}, the dimension is collapsed, and no packing is needed, so the initial declaration is preserved. If the dimension is aligned, there are three corresponding coordinates (p, ℓ, c). For every processor p, local (ℓ, c) pairs have to be packed onto a smaller area. This is done by first packing up the elements along the columns, then by removing the empty ones. A dual method is to pack first along the lines, then removing the empty ones. Formulae are derived below to perform this task.

Let’s consider the result of the above symbolic resolution, when the subscript expression is the identity (σ = 1, a_{0} = 0). The equation between u and the free variables of x is obtained by selecting 0 ≤ p ≤ π − 1, 0 ≤ ℓ ≤ γ − 1

4This assumption is, of course, not necessary for the general algorithm described in Section 4.1.

⌊\frac{t_{0}}{γπ} ⌋ ≤ c ≤ ⌊\frac{α(s − 1) + t_{0}}{γπ} ⌋

array A'(0:⌈\frac{c_{up}−c_{low}+1}{\frac{α}{g}}⌉,0:⌈\frac{γ}{g} ⌉)

Figure 6: Local new declaration

the first |x′| lines of H′. If H′′ = (I 0)H′ is the selected matrix, we have x′ − x′_{0} = H′′u, i.e.:

( \begin{matrix}p \\ ℓ − t_{0} \\ c \end{matrix} ) = ( \begin{matrix}1 & 0 & 0 \\ −γ & −g & 0 \\ 0 & μ & \frac{α}{g} \end{matrix} ) ( \begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix} )

Variables (u_{2}, u_{3}) could be used to define an efficient packing, since holes are removed by the change of basis (Figure 5) but the non diagonal terms of H′′ induce poor array declaration bounds. The space has to be skewed in order to match a rectangle as closely as possible, if space is at a premium, and if more complex, non-affine access functions are acceptable. Let G be the positive diagonal integer matrix composed of the absolute value of the diagonal coefficients of H′′.

alloc : u′ = ⌊G^{−1}(x′ − x′_{0})⌋ = ⌊G^{−1}H′′u⌋ (12)

G provides the right parameters to perform the proposed packing scheme. To every u a vector u′ is associated through formula 12. This formula introduces an integer division. Let’s show why u′ induces a better mapping of array elements on local memories than u. Since H′′ is lower triangular, formula 12 can be rewritten:

∀i ∈ [1 … |u|], u′_{i} = ±u_{i} + \frac{∑_{j=1}^{i−1} h_{i,j}u_{j}}{|h_{i,i}|}  (13)

This shows that alloc(u) is bijective. If u^{a} and u^{b} are different vectors, and i is the first dimension for which they differ, formula 13 shows that u_{i}^{′a} and u_{i}^{′b} will also differ. The function is also surjective, since the property allows to construct a vector that matches any u′ by induction on i.

Two of the three components of x′, namely p and ℓ, are explicitly bounded in K. Implicit bounds for the third component, c, are obtained by projecting K on c. These three pairs of bounds, divided by G, are used to declare the local part of the hpf array. Figure 6 shows the resulting declaration for the local part of the array, in the general symbolic case, for one dimension of the array.

The packing scheme induced by u′ is better than the one induced by u because there are no nondiagonal coefficients between u′ and x′ that would

α = 3, t_{0} = 0, σ = 3, a_{0} = 0, γ = 4, π = 4, β = 15

∀p_{0} ∈ [0 … 3] do u′_{2} = 0, 3

lb_{3} = \frac{16p_{0}+4u′_{2}+8}{9}, ub_{3} = \frac{16p_{0}+4u′_{2}+7}{9} lb′_{3} = \frac{9lb_{3}−4u′_{2}−16p_{0}}{3} do u′_{3} = lb′_{3}, lb′_{3} + (ub_{3} − lb_{3})

…

Figure 7: Simplified code example

introduce a waste of space, and u′ is as packed as u because contiguity is preserved by formula 13. The row packing scheme would have been obtained by choosing x′ = (p, c, ℓ) instead of x′ = (p, ℓ, c).

These packing schemes define two parameters (u′_{2}, u′_{3}) to map one element of one dimension of the hpf array to the processor’s local part. The local array declaration can be linearized with the u′_{3} dimension first, if the Fortran limit of 7 dimensions is exceeded.

4.4 Properties

The proposed iteration enumeration and packing scheme has several properties that should be emphasized. This scheme is compatible with efficient cache exploitation and overlap analysis. Moreover, some improvements can statically enhance the generated code.

If loop u_{2} in Figure 4 is reversed, the iteration enumeration order and the packing scheme are such that accesses are contiguous. Thus the local cache, if any, is efficiently used.

The packing scheme is also compatible with overlap analysis techniques [15]. Local array declarations are extended to provide space for border elements that are owned by neighbor processors, and to simplify accesses to non-local elements. The overlap is induced by relaxing constraints on ℓ, which is transformed through the scheme as relaxed constraints on u′_{2}. This allows overlaps to be simply considered by the scheme. Moreover a translated access in the original loop leading to an overlap is transformed into a translated access in the local spmd generated code.

The generic code proposed in Figure 4 can be greatly improved in many cases. Integer division may be simplified, or performed efficiently with shifts, or even removed. Node splitting and loop invariant code motion should be used to reduce the control overhead. Constraints may also be simplified, for instance if the concerned elements just match a cycle. Moreover, it is possible to generate the loop nest directly on u′, when u is not used in the loop body. For the main example in [10], such transformations produce the code shown in Figure 7.

In the general resolution (Section 4.1) the cycle variables c were put after the local offsets ℓ. The induced inner loop nest is then on c. It may be interesting to exchange ℓ and c orders when the block size is larger than the number of cycles: the loop with the larger range would then be the inner one. This may be useful when elementary processors have some vector capability.

4.5 Allocation of temporaries

Temporary space must be allocated to hold nonlocal array elements accessed by local iterations. For each loop L and array X, this set is infered from compute_{L}(p) and from subscript function S_{X}. For instance, references to array Y and Z in Figure 10 require local copies. Because of Y’s distribution, overlap analysis is fine but another algorithm is necessary for Z.

Let’s consider the generic loop in Figure 2, and assume that local elements of Y cannot efficiently be stored using overlap analysis techniques. First of all, if all or most of local elements of the lhs reference are defined by the loop nest, and if S_{Y} is of rank |i|, temporaries can be stored as X. If S_{Y} is of rank |i|, and if the access function S_{X} use one and only one index in each dimension, the resulting access pattern is still hpf like, so the result of Section 4.3 can be used to allocate the temporary array. Otherwise, another multi-stage change of coordinates is required to allocate a minimal area.

The general idea of the temporary allocation scheme is first to reduce the number of necessary dimensions for the temporary array, then to use this new basis for declaration, or a related one through a new change of coordinates if it is possible and desired to reduce the allocated space.

The set of local iterations, compute(p), is now defined by a new basis and new constraints, such that x − x_{0} = H′u (Equation (10)) and K′u ≤ k′. Some lines of H′ define iteration vector i, which is part of x (Equation (7)). Let H′_{i} be this sub-matrix: i = H′_{i}u.

Let H_{Y} = P_{Y}S_{Y}Q_{Y} be S_{Y}’s Hermite normal form. Let v = QY^{−1}i be the new parameters, then S_{Y}i = PY^{−1}H_{Y}v. Vector v can be decomposed as (v′, v′′) where v′ = Πv contribute to the computation of i and v′′ belongs to H_{Y}’s kernel. If H_{YL} is a selection of the non-zero columns of H_{Y} = (H_{YL}0), then we have:

a_{Y} − a_{Y0} = PY^{−1}H_{YL} v′
v′ = ΠQY^{−1} i

and by substituting i:

v′ = ΠQY^{−1}H′_{i}u

input: {y = f.w, W w ≤ 0}

output: {y′ = f′.w}

initial system: f^{0} = f for i = 1 … (|w| − 1)

g_{i} = gcd({f_{j}^{i−1}, j > i})

s_{i} = extent(∑_{j≤i} f_{j}^{i−1}w_{j}, W)

if (s_{i} ≥ g_{i})

then f^{i} = f^{i−1}

else ∀j ≤ i, f_{j}^{i} = f_{j}^{i−1} and ∀j > i, f_{j}^{i} = \frac{f_{j}^{i−1}}{g_{i}} s_{i} end for f′ = f^{|w|−1}

Figure 8: Heuristic to reduce allocated space

It is sometime possible to use x′ instead of u′. For instance, if the alignments and the subscript expressions are translations, i is an affine function of x′ and v′ simply depends on x′. The allocation algorithm is based on u but can also be applied when vxp is used instead.

Since the local allocation does not depend on the processor identity, the first |p| components of u should not appear in the access function. This is achieved by decomposing u as (u′, u′′) with |u′| = |p| and ΠQY^{−1}H′_{i} in (F_{P} , F_{Y}) such that:

v′ = F_{P} u′ + F_{Y}u′′

and the local part v′′ is introduced as:

v′′ = F_{Y}u′′

Then the first solution is to compute the amount of memory to allocate by projecting constraints K onto v′′. This is always correct but may lead to a waste of space because periodic pattern of accesses are not recognized and holes are allocated. In order to reduce the amount of allocated memory, a heuristic, based on large coefficients in F_{Y} and constraints K is used.

Each dimension, i.e. component of v′′, is treated independently. Let FY^{i} be a line of F_{Y} and y the corresponding component of v′′: y = fv′′. A change of basis G is applied to v′′ to reduce f to a simpler equivalent linear form f′ where each coefficient appears only once and where coefficients are sorted by increasing order. For instance:

f = ( 16 3 −4 −16 0 )

is replaced by:

f = f′G
G = ( \begin{matrix}0 & 1 & 0 & 0 & 0 \\ 0 & 0 & −1 & 0 & 0 \\ 1 & 0 & 0 & −1 & 0 \end{matrix} )
f′ = ( 3 4 16 )

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; i; a_{X}; a_{Y}; xX^{′}; xY^{′}; aX^{′}; aY^{′}; u; aY^{′′}; F_{i}; F_{x′}; S_{X}; S_{Y}; hpf; hpf; aY^{′′′}; alloc_{X}; alloc_{Y}; M′′; M′′′]

Figure 9: Relationships between frames

and v′′ is replaced by w = Gv′′. Constraints K on u are rewritten as constraints W on w by removing the constraints due to processors and by using projections and G.

Linear form f is then processed by the second part of the heuristic, shown in figure 8. It reduces the extent by \frac{s_{i}}{g_{i}} at each step, if the constraints W shows that there is a hole. A hole exists when g_{i} is larger than the extent of the partial linear form being built, under constraints W. Linearity of access to temporary elements is preserved.

The correctness of the scheme is shown by proving that a unique location y′ is associated to each y. Further insight on this problem, the minimal covering of a set by interval congruences, can be found in Granger [17] and Masdupuy [26].

4.6 Data movements

The relationships between the bases and frames defined in the previous sections are shown in Figure 9. Three areas are distinguished. The top one contains user level bases for iterations, i, and array elements, a_{X}, a_{Y},... The middle area contains the bases and frames used by the compiler to enumerate iterations, u, and to allocate local parts of hpf arrays, aX^{′}, aY^{′},... as well as the universal bases, xX^{′} and xY^{′}, used to define the lattices of interest, F and F′. The bottom area shows new bases used to allocate temporary arrays, aY^{′′} and aY^{′′′}, as defined in Section 4.5.

Solid arrows denote affine functions between different spaces and dashed arrows denote possibly non-affine functions built with integer divide and modulo operators. A quick look at these arrows shows that u, which was defined in Section 4.1 is the central basis of the scheme, because every other coordinates can be derived from u, along one or more paths.

Two kinds of data exchanges must be generated: updates of overlap areas, and initializations of temporary copies. Overlap areas are easier to handle because the local parts are mapped in the same way on each processor, using the same alloc function. Let’s consider array X in Figure 9. The send and receive statements are enumerated from the same polyhedron, up to a permutation of p and p′, with the same basis u. To each u corresponds only one element in the user space, a_{X}, by construction^{5} of u. To each a_{X} corresponds only one aX^{′} on each processor^{6}. As a result, data exchanges controled by loop on u enumerate the same element at the same iteration.

Because the allocation scheme is the same on each processor, the inner loop may be transformed into a vector message if the corresponding dimension is contiguous and/or the send/receive library supports constant but non-unit strides. Block copies of larger areas also are possible when alloc is affine, which is not the general case from a theoretical point of view but should very often happen in real applications.

Temporary arrays like Y′′ and Y′′′ are more difficult to initialize because there is no such identity between their allocation function as local part of an hpf array, alloc_{Y}, and their allocation functions as temporaries, t alloc_{Y}′′ and t alloc_{Y}′′′. However, basis u can still be used. When the temporary is allocated exactly like the master array, e.g. Y′′′ and X′, any enumeration of elements u in U enumerates every element aY^{′′′} once and only one because the input loop is assumed independent. On the receiver side, aY^{′′′} is directly derived from u. On the sender side, a_{Y} is computed as S_{Y}F_{i}u and a non-affine function is used to obtain xY^{′} and aY^{′}. Vector messages can be used when every function is affine, because a constant stride is transformed into another constant stride, and when such transfers are supported by the target machine. Else, calls to packing and unpacking routines call be generated.

When the temporary is allocated with its own allocation function M, it is more difficult to find a set of u enumerating exactly once its elements. This is the case for copy Y′′ in Figure 9. Elements of Y′′, aY^{′′}, must be used to compute one related u among many possible one. This can be achieved by using a pseudo-inverse of the access matrix M:

u = M^{t}(MM^{t})^{−1}aY^{′′}  (14)

Vector u is the rational solution of equation aY^{′′} = Mu which has a minimum norm. It may well not belong to U, the polyhedron of meaningful u which

5As explained in Section 4.3, allocation is handled as a special case by choosing the identity for SX.

6Elements in overlap area are allocated more than once, but on different processors.

      real X(0:42,0:42),
     &     Y(0:42,0:42), Z(0:42)

!HPF$ template T(0:150,0:150)
!HPF$ align X(I,J) with T(3*I,3*J)
!HPF$ align Y with X
!HPF$ align Z(I) with T(0,0)
!HPF$ processors P(0:3,0:3)
!HPF$ distribute                      &
!HPF$   T(cyclic(4),cyclic(4)) onto P

!HPF$ independent(I,J)
      do I = 1, 14
         do J = 0, 14
            X(3*I,3*J) = Y(3*I,3*J) +
     &        Y(3*I-1,3*J) + Z(3*I)
         enddo
      enddo
!HPF$ independent(I)
      do I = 0, 14
         Y(I,I) = 1.0
      enddo

Figure 10: Code example.

are linked to a user iteration. However, M was built to have the same kernel as S_{Y}F_{i}. Thus the same element a_{Y} is accessed as it would be by a regular u. Since only linear transformations are applied, these steps can be performed with integer computations by multiplying Equation (14) with the determinant of MM^{t} and by dividing the results by the same quantity. The local address, aY^{′}, is then derived as above. This proves that sends and receives can be enumerated by scanning elements aY^{′′}.

5 Examples

The different algorithms presented in the previous section can be used to distribute the contrived piece of code of Figure 10.

This is an extension of an example given in [10] showing that allocation of hpf arrays may be nontrivial. The reference to X in the first loop requires an allocation of X′ and the computation of new loop bounds. It is not a simple case because the subscript function is not the identity and because a block-cyclic distribution is specified. The two references to Y in the first loop are similar but they imply some data exchange between processors and the allocation of an overlap area in Y′. The reference to Z implies the allocation of a temporary array, Z′, which is independent from the local part of Z. The second loop shows that the compute(p) set may have fewer dimensions than the array referenced and that fewer loops have to be generated.

                                                                c     Local declaration of Y(0:42,0:42):
                                                                      real X_hpf(0:7,0:1,0:7,0:1)
                                                                      real Y_hpf(0:7,-1:1,0:7,0:1)
                                                                c     For each processor (0<=p1<=3,0<=p2<=0) :
                                                                c     do I = 1, 14
                                                                c        do J = 1, 14
                                                                c     ...
                                                                c     Send view(Y) to compute(X):
                                                                      do u3 = 0, min(p1+1, 3)
                                                                        if(u3.ne.p1) then
                                                                          do c1 = 0, 7
                                                                            do c2 = 0, 7
                                                                              do u6 = max(2, divide(2+4*p1+16*c1,3)),
                                                                     &                min(42, divide(3+4*p1+16*c1,3))
                                                                                if (divide(u6+2,3).le.divide(u6+1,3)) then
                                                                                  do u7 = max(0,divide(8+4*p2+16*c2,9)),
                                                                     &                    min(14, divide(3+4*p2+16*c2,9))
                                                                                    l1 = 3 * u6 - 16 * c1 - 4 * p1
                                                                                    l2 = 9 * u7 - 16 * c2 - 4 * p2
                                                                                    send(u3, p2, Y_hpf(c1, l1, c2, l2))
                                                                                  enddo
                                                                                endif
                                                                              enddo
                                                                            enddo
                                                                          enddo
                                                                        endif
                                                                      enddo
                                                                c     ...
                                                                c     do I=0,14 Y(I,I)=1.0
                                                                      do u3 = 0, 2
                                                                        do u4 = 0, 2
                                                                          do u5 = max(max(0,divide(2 + 4*p1 + 16*u3,3)),
                                                                     &                divide(2 + 4*p2 + 16*u4,3)),
                                                                     &            min(min(14,divide(3 + 4*p1 + 16*u3,3)),
                                                                     &                divide(3 + 4*p2 +16*u4,3))
                                                                            w5 = divide(3*u5 - 4*p1 - 16*u3,3)
                                                                            Y_hpf(u3, w5, u4, u6) = 1.0
                                                                          enddo
                                                                        enddo
                                                                      enddo
                                                                      end
                                                                c     Integer Fortran divide has not the good semantics:
                                                                      integer function divide(i,j)
                                                                      integer i,j
                                                                      if(i.ge.0) then
                                                                         divide = i/j
                                                                      else
                                                                         divide = -(-i+j-1)/j
                                                                      endif
                                                                      end

Figure 11: Output code (excerpts).

The output code is too long to be printed here. Interesting excerpts are shown in figure 11 and commented below.

5.1 HPF Declarations

hpf declarations already are normalized. Templates are fully utilized and Π = I. Other key parameters are:

P = ( \begin{matrix}4 & 0 \\ 0 & 4 \end{matrix} ) , T = ( \begin{matrix}150 & 0 \\ 0 & 150 \end{matrix} ) , C = ( \begin{matrix}4 & 0 \\ 0 & 4 \end{matrix} )

A_{X} = A_{Y} = ( \begin{matrix}3 & 0 \\ 0 & 3 \end{matrix} ) , R_{X} = R_{Y} = ( \begin{matrix}1 & 0 \\ 0 & 1 \end{matrix} )

t_{0X} = t_{0Y} = ( \begin{matrix}0 \\ 0 \end{matrix} ) , A_{Z} = ( \begin{matrix}0 \\ 0 \end{matrix} ) , t_{0Z} = ( \begin{matrix}0 \\ 0 \end{matrix} ) ,

5.2 Allocation of X′

Array X is referenced only once, as a lhs. Thus its local part is directly derived from the set own_{X}(p) using the algorithm described in section 4.3.

own_{X}(p) = { t_{1} = 3x_{1}, t_{2} = 3x_{2}

t_{1} = 16c_{1} + 4p_{1} + ℓ_{1} t_{2} = 16c_{2} + 4p_{2} + ℓ_{2} 0 ≤ x_{1}, x_{2} ≤ 42, 0 ≤ p_{1}, p_{2} ≤ 3 0 ≤ ℓ_{1}, ℓ_{2} ≤ 3, 0 ≤ t_{1}, t_{2} ≤ 150 }

The equations are parametrically solved and x is replaced by u in the constraints according to the following change of basis:

\begin{matrix}u_{1} = p_{1}, & 3u_{5} − 16u_{3} − 4u_{1} = ℓ_{1} \\ u_{2} = p_{2}, & 3u_{6} − 16u_{4} − 4u_{2} = ℓ_{2} \\ u_{3} = c_{1}, & u_{5} = x_{1} \\ u_{4} = c_{2}, & u_{6} = x_{2} \end{matrix}

Then, frame u is replaced by frame u′ to avoid non-diagonal terms in the constraints. Frame u′ is linked to the user basis x by:

\begin{matrix}u′_{1} = p_{1}, & u′_{2} = p_{2}, & u′_{3} = c_{1} \\ u′_{4} = c_{2}, & 3u′_{5} = ℓ_{1}, & 3u′_{6} = ℓ_{2} \end{matrix}

and constraints K are projected on each component of u′ to derive the declaration bounds of X′:

\begin{matrix}0 ≤ u′_{1} ≤ 3, & 0 ≤ u′_{2} ≤ 3, & 0 ≤ u′_{3} ≤ 7 \\ 0 ≤ u′_{4} ≤ 7, & 0 ≤ u′_{5} ≤ 1, & 0 ≤ u′_{6} ≤ 1 \end{matrix}

Coordinates u′_{1} and u′_{2} are the processor coordinates. The other coordinates, u′_{3} to u′_{6}, are the local array coordinates. The Fortran declaration of X′ is shown in Figure 11.

16×16 elements are allocated while 11×11 would be enough. The efficiency factor would be better in the (p, ℓ, c) basis. It would also be better with less contrived aligments and subscript expressions and/or larger blocks.

5.3 Sending Y′

Before the first loop nest is executed, each processor (p_{1}, p_{2}) must send some elements of Y to its right neighbor (p_{1} + 1, p_{2}), with a wrap-around for the rightmost processors, (3, p_{2}). The pairs of communicating processors are not accurately represented by convex constraints and the u3 loop is not as tight as it could. This suggests that central and edge processors should be handled differently when overlap is observed.

The bounds for the inner loops are complicated but they could mostly be evaluated iteratively if node splitting and loop invariant code motion are used. Also, the test on u6 can be simplified as a modulo operation. Finally l1 and l2 can also be evaluated iteratively. Note that the dimensions of Y′ should be permuted to reduce the stride for a possible vector message along u7.

5.4 Local iterations for the second loop

The second loop cannot be represented as an array section assignment, which makes it harder to handle than the first one. The set of local iteration is defined by:

compute_{Y}(p) = { y_{1} = i, y_{2} = i
t_{1} = 3y_{1}, t_{2} = 3y_{2}
t_{1} = 16c_{1} + 4p_{1} + ℓ_{1}
t_{2} = 16c_{2} + 4p_{2} + ℓ_{2}
0 ≤ i ≤ 14, 0 ≤ y_{1}, y_{2} ≤ 42
0 ≤ p_{1}, p_{2} ≤ 3, 0 ≤ ℓ_{1}, ℓ_{2} ≤ 3
0 ≤ t_{1}, t_{2} ≤ 150 }

The loop bounds were retrieved by projection on u. It was found useless to generate a guard on the processor identity because non-involved processor have nothing to do and because this does not delay the other ones. The guard is hidden in the inner loop bounds.

Note that an extra-loop is generated. Y diagonal can be enumerated with only two loops and three are generated. This is due to the use of an imprecise projection algorithm but does not endanger correctness [2]. Further work is needed in this area.

5.5 Integer divide

One implementation of the integer divide is finally shown. The divider is assumed strictly positive, as is the case in all call sites. It necessary because Fortran remainder is not positive for negative numbers. It was added to insure the semantics of the output code.

6 Related work

Techniques to generate distributed code from sequential or parallel code using a uniform memory space have been extensively studied since 1988 [9, 28, 39, 30, 15, 22, 16, 6, 25, 2, 21, 27, 31, 23, 5, 19, 4, 8, 11, 38].

The most obvious, most general and safest technique is called run-time resolution[9, 28, 30]. Each instruction is guarded by a condition that is true for processors that must execute it. Each memory address is checked before it is referenced to decide whether the address is local and the reference is executed, whether it is remote, and a receive is executed, or whether it is remotely accessed and a send is executed. This rewriting scheme is easy to implement [11] but very inefficient at run-time because guards, tests, sends and receives are pure overhead. Moreover every processor has to execute the whole control flow of the program, and even for parallel loop, communications may sequentialize the program at run-time [9].

Many optimization techniques have been introduced to handle specific cases. Gerndt introduced overlap analysis in [15] for block distributions. When local array parts are allocated with the necessary overlap and when parallel loops are translated, the instruction guards can very often be moved in the loop bounds and the send/receive statements are globally executed before the local iterations, i.e. the loop nest with run-time resolution is distributed into two loop nests, one for communications and one for computations. The communication loops can be rearranged to generate vector messages.

Tseng [35] presents lots of additional techniques (message aggregation and coalescing, message and vector message pipelining, computation replication, collective communication…). He assumes affine loop bounds and array subscripts to perform most optimizations. He only handles block and cyclic(1) distributions and the alignment coefficient must be 1.

Three recent publications tackle any alignment and distribution but restrict references to array sections. Each dimension is independent of the others as was assumed in Section 4.2. Chatterjee et al. [10] developed a finite state machine approach to enumerate local elements. No memory space is wasted and local array elements are ordered by Fortran lexico-graphic order exactly like user array elements. They are sequentially accessed by while loops executing the fsm’s, which may be a problem if vector units are available. Moreover, accesses to an auxiliary data structure, the fsm transition map, add to the overhead. Note that the code generated Figure 4 may be used to compute the fsm. In fact the lower iteration of the innermost loop is computed by the algorithm that constructs the fsm.

Two papers by Stichnoth [33] on one hand and Gupta et al. [18] on the other present two similar methods to solve the same problem. They use array sections but compute some of the coefficients at run-time. Gupta et al. solve the block distribution case and use processor virtualisation to handle cyclic distributions. Arrays are densely allocated as in [10] and the initial order is preserved but no formulae are given. Stichnoth uses the dual method for array allocation as in [11], that is blocks are first compressed, and the cycle number is used as a second argument.

7 Conclusion

Although the usual affine assumptions are made for loop bounds and subscript expressions, our compilation scheme simultaneously lifts several restrictions.

Array references are not restricted to array sections. General hpf alignments and distributions are supported, and the same algorithms also provide optimal code for classical block distributions. Memory allocation is almost 100 % efficient on large blocks and performs quite well on small ones when strange alignments are used. We believe that this slight memory waste is more than compensated by the stride-1 vector load, store, send and receive which can be performed on the copies and which are necessary for machines including vector units. These contiguous accesses also perform well with a cache. The overlap analysis and allocation is integrated to the basic allocation scheme. Finally, most computations are performed at compile-time and no auxiliary data structures are used.

Many partial optimization techniques are integrated in our direct synthesis approach: message vectorization, and aggregation [19], overlap analysis [15]. A new storage management scheme is also proposed. Moreover other optimizations techniques may be applied to the generated code such as vectorization [37], loop invariant code motion [1] and software pipelining [14, 36].

Future work includes an implementation in our hpf compiler [11], extensions to optimize sequential loops, including i/o loops, to overlap communication and computation, and to handle indirections.

Acknowledgments

We are thankful to Pierre Jouvelot for his careful reading and the many improvements he suggested.

References

[1] Alfred V. Aho, Ravi Sethi, and Jeffrey D. Ullman. Compilers Principles, Techniques, and Tools. Addison-Wesley Publishing Company, 1986.

[2] Corinne Ancourt. Génération automatique de codes de transfert pour multiprocesseurs à mémoires locales. PhD thesis, Université Paris VI, March 1991.

[3] Corinne Ancourt and François Irigoin. Scanning Polyhedra with DO Loops. In Symposium on Principles and Practice of Parallel Programming, April 1991.

[4] Corinne Ancourt and François Irigoin. Automatic Code Distribution. In Third Workshop on Compilers for Parallel Computers, July 1992.

[5] Françoise André, Olivier Chéron, and Jean-Louis Pazat. Compiling Sequential Programs for Distributed Memory Parallel Computers with Pandore II. Technical report, IRISA, May 1992.

[6] Françoise André, Jean-Louis Pazat, and Henry Thomas. Pandore: A System to Manage Data Distribution. Publication Interne 519, IRISA, February 1990.

[7] Thomas Brandes. Efficient Data Parallel Programming without Explicit Message Passing for Distributed Memory Multiprocessors. Internal Report AHR-92 4, High Performance Computing Center, German National Research Institute for Computer Science, September 1992. ftp ftp.gmd.de,

gmd/adaptor.

[8] Thomas Brandes. Adaptor: A Compilation System for Data Parallel Fortran Programs. Technical report, High Performance Computing Center, German National Research Institute for Computer Science, August 1993. ftp ftp.gmd.de,

gmd/adaptor/docs.

[9] David Callahan and Ken Kennedy. Compiling Programs for Distributed-Memory Multiprocessors. The Journal of Supercomputing, 2:151–169, 1988.

[10] Siddhartha Chatterjee, John R. Gilbert, Fred J. E. Long, Robert Schreiber, and Shang-Hua Teng. Generating Local Addresses and Communication Sets for Data-Parallel Programs. In Symposium on Principles and Practice of Parallel Programming, 1993.

[11] Fabien Coelho. Étude de la Compilation du High Performance Fortran. Master’s thesis, Université Paris VI, September 1993. Rapport de DEA Systèmes Informatiques, stage effectué au CRI, EN-SMP. Rapport EMP-CRI E/178.

[12] Paul Feautrier. Parametric Integer Programming. RAIRO Recherche Opérationnelle, 22:243– 268, September 1988.

[13] High Performance Fortran Forum. High Performance Fortran Language Specification. Version 1.0 CRPC-TR 92225, Center for Research on Parallel Computation, Rice University, Houston, USA, May 1993.

[14] Franco Gasperoni and Uwe Scheiegelshohn. Scheduling Loop on Parallel Processors: A Simple Algorithm with Close to Optimum Performance. Lecture Note, INRIA, 1992.

[15] Hans Michael Gerndt. Automatic Parallelization for Distributed-Memory Multiprocessing Systems. PhD thesis, University of Vienna, 1989.

[16] Hans Michael Gerndt and Hans Peter Zima. Optimizing Communication in Superb. In CONPAR90, pages 300–311, 1990.

[17] Philippe Granger. Analyses Sémantiques de Congruence. PhD thesis, École Polytechnique, July 1991.

[18] S. K. S. Gupta, S. D. Kaushik, S. Mufti, S. Sharma, C.-H. Huang, and P. Sadayappan. On compiling Array Expressions for Efficient Execution on Distributed-Memory Machines. In International Conference on Parallel Processing, pages II–301–II–305, August 1993.

[19] Seema Hiranandani, Ken Kennedy, and Chau-Wen Tseng. Compiling Fortran D for MIMD Distributed-Memory Machines. Communications of the ACM, 35(8):66–80, August 1992.

[20] François Irigoin. Interprocedural Analyses for Programming Environment. In Jack J. Dongara and Bernard Tourancheau, editors, Environments and Tools for Parallel Scientific Computing, pages 333–350, Saint-Hilaire-du-Touvet, September 1992. North-Holland, Amsterdam, NSF-CNRS.

[21] Charles Koelbel and Piyush Mehrotra. Compiling Global Name-Space Parallel Loops for Distributed Execution. IEEE Transactions on Parallel and Distributed Systems, 2(4):440–451, October 1991.

[22] Charles Koelbel, Piyush Mehrotra, and John Van Rosendale. Supporting Shared Data Structures on Distributed Memory Architectures. Technical Report ASD 915, Purdue University, January 1990.

[23] Oded Lempel, Shlomit S. Pinter, and Eli Turiel. Parallelizing a C Dialect for Distributed Memory MIMD Machines. In Language and Compilers for Parallel Computing, August 1992.

[24] John Levesque. FORGE 90 and High Performance Fortran. Applied Parallel Research, Inc., 1992. xHPF77 presentation.

[25] J. Li and Marina Chen. Compiling communication-efficient programs for massively parallel machines. IEEE Transactions on Parallel and Distributed Systems, 2(3):361–376, July 1991.

[26] François Masdupuy. Array Indices Relational Semantic Analysis using Rational Cosets and Trapezoids. PhD thesis, École Polytechnique, November 1993. Technical Report EMP-CRI A/247. To be defended.

[27] John Merlin. Techniques for the Automatic Parallelisation of ‘Distributed Fortran 90’. SNARC 92 02, University of Southampton, 1992.

[28] Ravi Mirchandaney, Joel S. Saltz, Roger M. Smith, David M. Nicol, and Kay Crowley. Principles of Runtime Support for Parallel Processors. In International Conference on Supercomputing, pages 140–152, July 1988.

[29] William Pugh. A Pratical Algorithm for Exact Array Dependence Analysis. CACM, 35(8):102–114, August 1992.

[30] Anne Rogers and Keshav Pingali. Process Decomposition Through Locality of Reference. In ACM SIGPLAN International Conference on Programming Language Design and Implementation, June 1989.

[31] Anne Rogers and Keshav Pingali. Compiling for distributed memory architectures, June 1992.

[32] Alexander Schrijver. Theory of linear and integer programming. Wiley, New-York, 1986.

[33] James M. Stichnoth. Efficient Compilation of Array Statements for Private Memory Multicomputers. CMU-CS-93 109, School of Computer Science, Carnegie Mellon University, February 1993.

[34] Rémi Triolet, Paul Feautrier, and François Irigoin. Direct Parallelization of Call Statements. In Proceedings of the ACM Symposium on Compiler Construction, 1986.

[35] Chau-Wen Tseng. An Optimising Fortran D Compiler for MIMD Distributed Memory Machines. PhD thesis, Rice University, Houston, Texas, January 1993.

[36] Jian Wang and Christine Eisenbeis. Decomposed Software Pipelining: A New Approach to Exploit Instruction Level Parallelism for Loop Programs. In Working Conference on Architectures and Compilation Techniques for Fine and Medium Grain Parallelism, January 1993.

[37] Hans Zima and Barbara Chapman. Supercompilers for Parallel and Vector Computers. ACM Press, 1990.

[38] Hans Zima and Barbara Mary Chapman. Compiling for Distributed-Memory Systems. Proceedings of the IEEE, February 1993.

[39] Hans Peter Zima, H. J. Bast, and Hans Michael Gerndt. SUPERB : A Tool for Semi-Automatic MIMD/SIMD Parallelization. Parallel Computing, 6:1–18, 1988.