Clémentin Tayou Djamegni

Patrice Quinton

Sanjay Rajopadhye

Tanguy Risset

Maurice Tchuenté

A reindexing based approach towards mapping of DAG with affine schedules onto parallel embedded systems

Keywords: parallel embedded systems, timing function, timing surface, allocation function, space-optimality, space-time complexity, reindexing, dependence graph, systolic arrays

We address the problem of optimally mapping uniform DAGs to systolic arrays, given an affine timing function. We introduce an automatic allocation method based on a preprocessing by reindexing that transforms the initial DAG into a new one that enables the well known projection method to minimize the number of processors along a number of directions. We demonstrate its superiority to other methods, and establish the space-optimality of the proposed method. We also show an upper bound on the number of processors that corresponds to the best space complexity that both the projection method and the so-called grouping method can give for the initial DAG. We also describe how the new allocation method can be implemented in tools.

Introduction

Embedded systems are information processing devices that are integral parts of complex applications: consumers buy them for the functionality they provide, not because there is a computer inside. Embedded systems are becoming more and more important from both the economic and technical side. They are tailored to do a particular task that can be squeezed into an extremely compact space, and are often specialized to a given task. The proliferation of computational devices have pushed for specialization to meet constraints on speed, power, energy, and area. For example, a mobile phone contains a specialized embedded architecture consisting of several different processors such as a DSP for wireless channel processing, a micro-controller that is responsible for the user interface and communication protocol, as well as several custom hardware components e.g., in the receiver/transmitter. The performance of this heterogeneous architecture exceeds that of any existing PC at lower power consumption, size, weight, and production cost. It is specialization that makes the mobile phone feasible and economical. Embedded systems are also ubiquitous in automotive electronics, robotics, aerospace applications and video game platforms, to name just a few prominent examples.

Many candidates for embedded systems in signal processing and multimedia applications are compute-bound. In order to meet higher performance computing speeds, an important option is to use highly parallel architectures (VLSI processor arrays, FPGA, etc.) This is shown by the relative maturity achieved by academic research systems such as MMAlpha [START_REF] Guillou | Hardware Design Methodology with the Alpha Language, Forum on Design Languages[END_REF], Compaan [START_REF] Kienhuis | Compaan: Deriving process networks from Matlab for embedded signal processing architectures[END_REF], VASS [START_REF] Yeh | VASS a VLSI array system synthesizer[END_REF], PARO [START_REF] Teich | Partitioning of processor arrays: a piecewise regular approach[END_REF][START_REF] Bednara | Interface Synthesis for FPGA Based VLSI Processor Arrays[END_REF] as well as industrial tools such as PICO-N system [START_REF] Schreiber | High-Level Synthesis of Nonprogrammable Hardware Accelerators[END_REF] developed at HP Labs, and now commercialized as PICO Express by Synfora. In these systems, the target (virtual) architecture is usually a systolic array, i.e., a locally interconnected grid of elementary processors called Processing Elements (PEs) or cells. In addition to full custom ASICs, FPGAs and other reconfigurable platforms also serve as physical target architectures for such PE arrays.

In designing systolic arrays, one usually first seeks the fastest schedule, and then chooses an allocation function that maps the computation to the PEs. The problem of minimizing the number of PEs, for a given optimal schedule, has been widely studied in the literature. The techniques proposed have been either linear allocation functions [START_REF] Moldovan | On the analysis and synthesis of VLSI algorithms[END_REF][START_REF] Quinton | Automatic synthesis of systolic array from uniform recurrence equations[END_REF][START_REF] Quinton | The mapping of linear equations on regular arrays[END_REF][START_REF] Shang | On mapping of uniform dependence algorithms into lower dimensional processors arrays[END_REF][START_REF] Wong | Space-optimal linear processor allocation for systolic arrays synthesis[END_REF], called projection methods, as well as non-linear allocation functions, namely grouping [START_REF] Bu | Processor clustering for the design of optimal fixed-size systolic arrays[END_REF][START_REF] Clauss | Calculus of space-optimal mappings of systolic algorithms on processor arrays[END_REF][START_REF] Darte | Synthesizing systolic arrays: Some recent developments[END_REF], instruction shifts [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF], piling [START_REF] Benaini | Space-minimal systolic arrays for Gaussian elimination and the algebraic path problem[END_REF][START_REF] Clauss | Calculus of space-optimal mappings of systolic algorithms on processor arrays[END_REF] and partition [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Tchuenté | Parallel Computation on regular arrays[END_REF]. Of these, the last three are special cases of piecewise linear allocation functions, where different linear allocations are applied to different regions of the domain. It has been conjectured that the optimal allocation function is always piecewise linear, but this has never been proved. The projection method, the grouping method and the instruction shift method can, and have been implemented in automatic synthesis tools, but the others have not, primarily because there is no clear way of describing the "pieces." Moreover, all the existing methods require human insight and are therefore not automatic. None of the methods yields provably space-optimal arrays in general [START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF][START_REF] Djamegni | Contribution to the synthesis of optimal algorithms for regular arrays[END_REF]. Nevertheless, they have all been used to obtain space-optimal arrays for many specific problem instances. Of these, the piling method usually produces arrays with spiral interconnections, leading to somewhat poorer locality. For many problems, the instruction shift method and the partition method [START_REF] Benaini | Space-minimal systolic arrays for Gaussian elimination and the algebraic path problem[END_REF][START_REF] Bermond | Parallelization of the Gaussian elimination on systolic arrays[END_REF][START_REF] Bermond | Parallelization of the Gauss elimination on systolic arrays[END_REF][START_REF] Cappello | A processor-time minimal systolic array for cubical mesh algorithms[END_REF][START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF][START_REF] Djamegni | Contribution to the synthesis of optimal algorithms for regular arrays[END_REF] have proved to be superior to the piling method, yielding arrays with fewer PEs, although the power of these methods have not been formally compared previously.

In this paper, we resolve the standing open problem of deriving spaceoptimal arrays for a given uniform dependence DAG with a given affine schedule. Our method is constructive and consists of two steps applied recursively, dimension by dimension. The steps are: (i) a sequence of re-indexations applied to specific pieces of the domain (our method describes how to determine the pieces); and (ii) projections along canonical directions. Because of this, we are able to describe how the method may be automated and incorporated in a tool such as MMAlpha [START_REF] Guillou | Hardware Design Methodology with the Alpha Language, Forum on Design Languages[END_REF]. The proof of space optimality de facto proves the superiority of this method over the previous ad hoc techniques. We also systematically compare our method with the projection and grouping methods. The idea behind this allocation method was originally introduced, for specific problem instances, by Tayou Djamegni [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF][START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF].

The rest of the paper is organized in the following way. Section 2 illustrates our method (without the dimension by dimension recursion) on a simple sorting example. In Section 3, we introduce the notation, summarize how parallel algorithms are synthesized from recurrence equations, and present previous allocation methods. Section 4 describes our method, and in Section 5 we illustrate it on a number of examples. In Section 6, we compare the new allocation method with other methods. Section 7 describes how to resolve an important problem if one is to implement the method in tools like MMAlpha [START_REF] Guillou | Hardware Design Methodology with the Alpha Language, Forum on Design Languages[END_REF], and we conclude in Section 8. (1) We first illustrate the main idea of the reindexing scheme through a simple sorting example. It has been deliberately chosen so that the recursion is not needed (the problem is two-dimensional). Figure 1.a shows the bubble sort algorithm and its triangular shaped dependence graph, together with the optimal timing function, t(i, j) = i + j. The maximum number of active points is n 2 , so an optimal allocation function should have n 2 PEs. All projections (including the standard ones-orthogonal along the axes, or along the diagonal) yield arrays with at least n PEs, and are suboptimal. However, note that the diagonal projection yields an array where every PE is active only on alternate clock cycles, and two adjacent PEs can be grouped into a single PE, with no slowdown. Thus, the projection method combined with grouping yields an optimal array. However, we do not know how to discover this projection automatically.

Illustrative Example

(i,j) (i,i+j) (b) Mapping the timing vector (1,1) into (0,1) 1 x1 2 . x2 3 . .
Another optimal allocation is obtained if we partition the triangle into a family of L-shaped lines and allocate all nodes along each such line to a single PE: the left and bottom boundary of the triangle would be computed on the leftmost PE, and so on. The corners of these L-shaped lines all lie on the line, i + j = n, and this allocation corresponds to the piecewise linear allocation function

a(i, j) = i + j ≤ n : i i + j > n : n -j
This is the essence of the partition method, and the key difficulty is that discovering the i + j = n line to cut the domain requires designer insight.

For our example, the instruction shift method and the piling method yield another array which is also optimal (both of them follow very similar reasoning for this example and produce the same array). To explain these methods, it is useful to first produce a suboptimal array, say by the canonical projection (as shown in Figure 1.b). Every vertical line in this "skewed" triangle is the temporal activity of a PE. We observe that PEs have different starting and ending times. In fact, the i-th PE starts at 2i -1 and finishes at n + i -1. This means that the i -th PE in the "right half of the array" (i.e., one whose coordinates are n 2 + i ) starts just as the i -th PE (in the left half) stops. Thus they can both be executed by a single physical PE.

The problem is that all these approaches need some human intuition. Our method finds the optimal allocation as follows. We first choose some arbitrary allocation, such as the one in Figure 1.b (the method is insensitive to this choice). Each row in this triangle represents the set of simultaneous activity of some subset of PEs. We translate each such line (i.e., reindex the domain) so that it is flush with the vertical axis. This yields the graph shown in Figure 1.c, and corresponds to our optimal array (Figure 1.d). The partition method also gives us the same array.

The remainder of this paper formally proves that this approach is guaranteed to produce arrays with the minimal number of PEs, and that the method can be generalized to higher dimensions.

Background

Algorithms that are good candidates for systolic implementations are characterized by regular and repetitive operations on large sets of data. They are usually written in the form of nested loops [START_REF] Peir | Minimum distance:a method for partitioning recurrences for multiprocessors[END_REF] and more generally in the form of a System of Recurrence Equations (SREs) [START_REF] Karp | The organization of computations for uniform recurrence equations[END_REF][START_REF] Rao | Regular iterative algorithms and their implementation on processors arrays[END_REF]. These formalisms and the problem they involve were first developed by Lamport [START_REF] Lamport | The parallel execution of Do loops[END_REF] and by Karp et al. [START_REF] Karp | The organization of computations for uniform recurrence equations[END_REF].

The standard methodology for systolic array synthesis developed over the past 25 years [START_REF] Moldovan | On the analysis and synthesis of VLSI algorithms[END_REF][START_REF] Quinton | Automatic synthesis of systolic array from uniform recurrence equations[END_REF][START_REF] Quinton | The mapping of linear equations on regular arrays[END_REF][START_REF] Rajopadhye | Synthesizing systolic arrays with control signal from recurrence equations[END_REF][START_REF] Rajopadhye | Synthesizing systolic arrays from Recurrences Equations with Linear Dependencies[END_REF][START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Shang | On mapping of uniform dependence algorithms into lower dimensional processors arrays[END_REF] proceeds in four steps.

• Specify the problem as an SRE or as a loop program.

• Uniformize/localize the dependences in the specification [START_REF] Quinton | The mapping of linear equations on regular arrays[END_REF][START_REF] Rajopadhye | Synthesizing systolic arrays with control signal from recurrence equations[END_REF][START_REF] Rajopadhye | Synthesizing systolic arrays from Recurrences Equations with Linear Dependencies[END_REF][START_REF] Thiele | Systematic design of local processor arrays for numerical algorithms[END_REF][START_REF] Wong | Transformation of broadcasts into propagations in systolic algorithms[END_REF]. This leads to a so called System of Uniform Recurrence Equations (SURE) which has a uniform dependence graph: a Directed Acyclic Graph (DAG), G = (D, U ) where D is the set of tasks and U the set of edges. D is the set of integer points in a polyhedral region of Z n , and U is a set of constant n-vectors.

• Define a timing function (or schedule) t : D → N that gives the computation date t(v) of each task v of D, assuming that all tasks are of unit delay.

• Define an allocation function D → Z n-1 that assigns the tasks to PE coordinates so as to avoid conflict, i.e., no two tasks with the same execution date are assigned to the same PE.

For the purpose of this paper, we assume that the first three steps have been performed, and we seek optimal allocation function(s).

Definition 1 A timing function is said to be optimal with respect to a dependence graph, G iff it minimizes the number of time steps needed to execute all the tasks of G, i.e., the total execution time is equal to the length of the longest path of G. We restrict ourselves to affine timing functions, i.e., functions of the form t(z) = λ t z + α, where z ∈ Z n , λ = (λ 1 , λ 2 , . . . , λ n ) t and λ 1 , λ 2 , . . . , λ n , α are integral constants. The vector λ is called the timing vector or scheduling vector.

The schedule vector λ is orthogonal to (n -1)-dimensional hyperplanes of Z n whose points have the same execution instant. These hyperplanes are called the timing surfaces. The maximum number of tasks on any timing surface defines the potential parallelism. To satisfy the causal dependencies, any affine timing function must satisfy the following dependence constraint.

Definition 2 A vector v = (v 1 , v 2 , . . . , v n ) t is
said to be unimodular (or primitive) iff its components are relatively prime, i.e., gcd(v 1 , v 2 , . . . , v n ) = 1. An integral square matrix U is said to be unimodular iff |det(U )| = 1. The affine timing function t = [λ, α] is normalized iff its linear part λ is unimodular.

Proposition 1 An affine timing function specified by [λ, α] is valid if it satisfies λ t d z > 0 for any dependence vector d z associated with the problem.

A linear allocation is realized by projecting the iteration domain along a direction called allocation direction or projection direction. It is defined by a constant n-dimensional vector ξ.

Proposition 2 Let [λ, α] be a valid affine timing function. A constant ndimensional vector ξ defines a valid allocation direction if and only if λ t ξ = 0. Definition 3 Given a timing function, we call potential parallelism the maximum number of tasks having the same computation date. An allocation function is said to be optimal with respect to a dependence graph G and a timing function iff it minimizes the number of PEs, i.e., the PE count is equal to the potential parallelism of the timing function.

Definition 4

The minimum number of PEs over all possible allocation functions described by projections, defines the potential parallelism by projection.

Throughout this paper we assume without loss of generality that the affine timing function t = [λ, α] is normalized and that vector ξ is unimodular. We now provide a number of definitions that allow us to study the activities of PEs and arrays.

Definition 5 A working date of a PE is any date at which that PE executes a task of the iteration domain D. The working interval of a PE is the interval defined by the first and last working dates of that PE. An activation date of an array is a date at which at least one PE executes a task of D. The working rate of a PE is the ratio of the number of working dates of the PE and its working interval. The pipelining rate of an array is defined as the average of the working rates of its PEs.

The working rate of a PE is 100% if it is active at all dates within its working interval. The pipelining rate of an array is 100% if the working rate of all its PEs is 100%. This is not often achieved. For example, the pipelining rates of both optimal sorting arrays described in the previous section, is 50%.

Minimizing the number of PEs

We now summarize previous work related to the minimization of the number of PEs in regular arrays.

1. Projection method [START_REF] Moldovan | On the analysis and synthesis of VLSI algorithms[END_REF][START_REF] Quinton | Automatic synthesis of systolic array from uniform recurrence equations[END_REF][START_REF] Quinton | The mapping of linear equations on regular arrays[END_REF][START_REF] Shang | On mapping of uniform dependence algorithms into lower dimensional processors arrays[END_REF][START_REF] Wong | Space-optimal linear processor allocation for systolic arrays synthesis[END_REF]. This method corresponds to linear allocations. As stated earlier it is realized by projecting the dependence graph G along a direction ξ. All the tasks belonging to a same line of direction ξ are assigned to a same PE. The main drawback of this allocation technique is the low PE utilization occurring in the resulting array [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF]. Wong and Delosme [START_REF] Wong | Space-optimal linear processor allocation for systolic arrays synthesis[END_REF], Shang and Fortes [START_REF] Shang | On mapping of uniform dependence algorithms into lower dimensional processors arrays[END_REF], Ganapathy and Wah [START_REF] Ganapathy | Optimal synthesis of algorithm-specific lower-dimensional processor arrays[END_REF] use integer programming to get the best linear allocation possible. However, such a linear allocation may not correspond to an optimal allocation.

For instance, the minimal number of PEs obtained by projecting the dependence graphs of the Cholesky Factorization (henceforth called CF) [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF], Triangular Matrix Inversion (henceforth TMI) [START_REF] Louka | Triangular matrix inversion on systolic arrays[END_REF] and Optimal String Parenthesization by Dynamic Programming (DP) [START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF] is

N 2 /2 + Θ(N )
where N is the problem size. For square Matrix Multiplication (MM) the minimal number of PEs is N 2 . None of these arrays is space optimal with respect to the earliest or the latest optimal timing functions. The potential parallelism of CF and MM are, respectively, N 2 /8 + Θ(N ) [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF] and 3N 2 /4 + Θ(N ) [START_REF] Cappello | A processor-time minimal systolic array for cubical mesh algorithms[END_REF] for any optimal timing function. The potential parallelism of the TMI and DP are, respectively, N 2 /4 + Θ(N ) [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Djamegni | Scheduling of the DAG associated with pipeline inversion of triangular matrices[END_REF] and N 2 /6 + Θ(N ) [START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF] for the latest optimal timing functions and, N 2 /6 + Θ(N ) [START_REF] Louka | Triangular matrix inversion on systolic arrays[END_REF] and N 2 /14 + Θ(N ) [START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF] for the earliest optimal timing functions.

2. Grouping (or clustering) method [START_REF] Bu | Processor clustering for the design of optimal fixed-size systolic arrays[END_REF][START_REF] Clauss | Calculus of space-optimal mappings of systolic algorithms on processor arrays[END_REF][START_REF] Darte | Synthesizing systolic arrays: Some recent developments[END_REF]. This approach seeks to reduce the number of PEs of an initial array obtained by projection stemming from the analysis of PEs' activity. It consists in grouping |λ t ξ| neighboring PEs with distinct working dates. Typically, the PEs that are clustered are along a direction specified by a vector, γ. This approach permits to reduce the number of PEs by a factor of |λ t ξ|. In addition, the resulting solution is a locally connected array that has a 100% pipelining rate. In The grouping method cannot guarantee that the design is space-optimal, because a 100% pipelining rate in the array does not necessarily imply space-optimality. Another drawback of the grouping method is that it applies only if |λ t ξ| ≥ 2.

For the CF, TMI and DP the grouping method reduces the number of PEs by a factor of 2, i.e., from N 2 /2 + Θ(N ) to N 2 /4 + Θ(N ). For TMI this is optimal according to the latest optimal timing function [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Djamegni | Scheduling of the DAG associated with pipeline inversion of triangular matrices[END_REF].

For MM it is easy to prove that the grouping method cannot do better than N 2 PEs.

3. Instruction shift method [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF]. As in the grouping method, the starting point is an array obtained from the projection method. Then, the initial array is partitioned into a number of PEs segments parallel to a direction called partition direction. The tasks of each segment are reallocated so as to minimize the number of PEs on each segment. This amounts to minimizing the number of PEs on all intersection of the dependence graph with a plane generated by both the projection direction and the partition direction. This approach does not guarantee space-optimality. For the CF this technique reduces the number of PEs from N 2 /2 + Θ(N ) to N 2 /6 + Θ(N ) [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF]. Although this represents a significant improvement, this design is not space optimal as the potential parallelism of the CF is N 2 /8 + Θ(N ) [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF].

4.

Piling method [START_REF] Benaini | Space-minimal systolic arrays for Gaussian elimination and the algebraic path problem[END_REF][START_REF] Clauss | Calculus of space-optimal mappings of systolic algorithms on processor arrays[END_REF]. First, from a given dependence graph and affine schedule, a set M of tasks is found such that all tasks in the set are scheduled to be executed at the same time and the set size |M | is maximal. Second, an allocation method is applied to assign the tasks of the dependence graph to PEs. Any PE that has not been assigned

to execute a task of M is piled to a PE that executes a task of M and has disjoint working intervals. However, piling PEs results in long communications such as spiral links, and increases irregularity for the resulting arrays. In this paper, we seek to avoid this. Another drawback of this method is that it is too general: as stated it can be applied to any task graph; exploiting it for systolic array synthesis requires human insight.

5.

Partition method [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Tchuenté | Parallel Computation on regular arrays[END_REF]. The starting point is to partition the dependence graph G into a number of sub-graphs of lower dimension. Then, the tasks of each sub-graph are allocated to PEs, in an ad hoc manner, so as to minimize the number of PEs on each sub-graph. When the subgraphs are of dimension two, this corresponds to the instruction shift method. Thus, the partition method can be seen as a generalization of the instruction shift method.

In most of the cases where it has been applied in the literature [START_REF] Bermond | Parallelization of the Gaussian elimination on systolic arrays[END_REF][START_REF] Bermond | Parallelization of the Gauss elimination on systolic arrays[END_REF][START_REF] Djamegni | Contribution to the synthesis of optimal algorithms for regular arrays[END_REF][START_REF] Djamegni | Derivation of systolic algorithms for the algebraic path problem by recurrence transformations[END_REF][START_REF] Djamegni | A new dynamic programming algorithm on two dimensional arrays[END_REF][START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF][START_REF] Louka | Triangular matrix inversion on systolic arrays[END_REF][START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Tchuenté | Parallel Computation on regular arrays[END_REF] the sub-graphs correspond to twodimensional domains that are obtained by considering the intersections of the dependence graph with a set of parallel planes. This allocation heuristic leads to an array of N 2 /8 + θ(N ) PEs for the DP [START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF]. In [START_REF] Djamegni | A new dynamic programming algorithm on two dimensional arrays[END_REF] the number of PEs is reduced to N 2 /10 + θ(N ) by merging nodes of the dependence graph associated with the DP before applying the partition method. All these solutions are not space-optimal as they are based on the earliest optimal timing function whose potential parallelism is N 2 /14 + Θ(N ) [START_REF] Louka | Dynamic programming on two-dimensional systolic arrays[END_REF]. The partition method is also used in [START_REF] Djamegni | Derivation of systolic algorithms for the algebraic path problem by recurrence transformations[END_REF] to derive a space-optimal array for the Algebraic Path Problem (APP) and in [START_REF] Bermond | Parallelization of the Gaussian elimination on systolic arrays[END_REF][START_REF] Bermond | Parallelization of the Gauss elimination on systolic arrays[END_REF] by Bermond et al. to derive various arrays for the Gaussian Elimination (GE). However these GE arrays are not space-optimal. A Space-time optimal array for the GE is proposed in [START_REF] Benaini | Space-minimal systolic arrays for Gaussian elimination and the algebraic path problem[END_REF]. For CF the partition method leads to an array of N 2 /6+Θ(N ) PEs [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF]. For TMI, the partition method leads to space-time optimal arrays for both the earliest and latest optimal timing functions [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Louka | Triangular matrix inversion on systolic arrays[END_REF]. The partition method also provides space-optimal arrays for MM. The main drawback of the partition method comes from the fact that the tasks of each sub-graph are not allocated in a systematic manner. Because of this the partition method is not suitable for tools that automatically generate application specific VLSI PEs.

We see that only the projection method, the grouping method and the instruction shift method are suitable for tools that (semi) automatically generate application specific VLSI PEs. Moreover, the allocation methods presented above do not guarantee space-optimality. The weakest allocation method, in terms of reducing the number of PEs, seems to be the projection method, while the most promising ones are the grouping method, the instruction shift method and the partition method. The partition method is a generalization of the instruction shift method. According to results related to CF, TMI, DP, GE and MM the most powerful method seems to be the partition method. However, this has not yet been confirmed by a systematic comparison. In others words, we are not aware of any procedure based on the partition method that systematically leads to results that are better than (or as good as) the best result that the grouping method (respectively, the projection method) can give.

This brief survey of previous allocation methods shows that the problem of minimizing the number of PEs occurring in regular arrays remains open despite considerable research.

Notation related to the new allocation method

The allocation strategy introduced in section 4 is based on reindexing transformations. We now show how these transformations act on a uniform dependence graph associated with an SURE. Let S be an SU RE defined over a domain D ⊂ Z n and associated with a dependence graph G = (D, U ). Assume that S admits a valid affine timing function t = [λ, α]. Definition 6 A function g : D → g(D) is a valid reindexing transformation for S iff it is a bijection.

The reindexing function g transforms system S into a new one S that is is obtained by replacing in system S each index point z with its image by g, i.e., g(z). The dependence graph G becomes G = (D , U ) where D = g(D) and

U = {[g(z 1 ) → g(z 2 )] | [z 1 → z 2 ] ∈ U, with z 1 ∈ D and z 2 ∈ D}.
The equivalent timing function t corresponding to the reindexing g assigns to z ∈ D the computation date of its pre-image by g, i.e., g -1 (z ). This means that t = t • g -1 .

Proposition 3 When the reindexing g corresponds to an unimodular transformation, i.e., g(z) = U z + z 0 where U is a constant unimodular matrix and

z 0 is a constant, the schedule [λ, α]] is transformed to [U -t λ, α -λ t U -1 z 0 ].
Furthermore, a reindexing maintains the potential parallelism under the new schedule [START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF].

Definition 7 A segment seg(z 1 , z 2 ), z 1 = z 2 , with extremities z 1 and z 2 is the set {z ∈ Z n | z = α 1 z 1 + α 2 z 2 , α 1 , α 2 ≥ 0, α 1 + α 2 = 1}.
A domain D is said to be convex along an n-dimensional primitive vector s if its intersection with any segment parallel to s is either empty or a single segment. D is convex if it is convex along all directions.

If D is convex along s, the frontier

F (D, s) of D along s, is the subset of D defined by F (D, s) = {z ∈ D | z + s / ∈ D}.
Definition 8 Given two linearly independent vectors v 1 and v 2 , a v 1 , v 2slice of D is any intersection of D with a 2-hyperplane parallel to vectors v 1 and v 2 , called supporting 2-hyperplane.

Assume that D is convex along s. For some z ∈ D, say that z s,D (the "shadow" of z along s on the frontier of D) is the unique index point of F (D, s) that belongs to the line defined by vector s and point z. Similarly, let z s denote the projection of point z ∈ Z n along s on the (n-1)-hyperplane orthogonal to s and containing point (0, 0, 0, . . . , 0) T ∈ Z n .

Consider the array obtained by projecting domain D along ξ. N P E(ξ, D) denotes the number of PEs in this array. The number of PEs obtained by applying the clustering method along γ is denoted by N P E(ξ, D, γ), and the number of PEs obtained by applying the instruction shift method along γ is denoted by N B(ξ, D, γ).

The new allocation method

We are now ready to present the new allocation method. As stated earlier, it is based on a preprocessing by reindexing that transforms the original domain into a new one that is more suitable for the application of the projection method. Given that the projection method assumes that points having the same computation date belong to the same (n -1)-hyperplane, we will consider only reindexing functions that preserve this property. Throughout this section, S denotes an SU RE over a bounded domain D ⊆ Z n , with a valid timing function t = [λ, α], where λ is normalized (otherwise, t would not be optimal). Let e 1 , e 2 , e 3 , . . . , e n be the canonical basis of Z n . As in all related work we assume that the timing surfaces are convex.

In order to simplify the presentation, we first develop the method on twodimensional iteration domains. Then, results obtained from this simple case are generalized to n-dimensional domains.

Two-dimensional Domains

Here the timing surfaces are line segments. Depending on the direction of the timing vector λ, we consider two cases.

Case 1: λ is collinear to a canonical basis vector e l .

We present here the reindexing that enables the projection method to provide an optimal allocation. Let e h be the other canonical basis and choose s to be either e h or -e h . The key idea is to shift the timing segments along s so as to move the points of F (D, s) to the e l axis. To do so, we use the reindexing function P h,D defined as follows:

P h,D (z) = z ∈ F (D, s) : z s z / ∈ F (D, s) : z s + z -z s,D (1)
= z s + e t h (zz s,D )e h It is easy to prove that P h,D is a valid reindexing transformation. Clearly P h,D (z) = (e t 1 (zz s,D ), e t 2 z) t for h = 1 and P h,D = (e t 1 z, e t 2 (zz s,D )) t for h = 2. The reindexing P 1,D is illustrated in figure 2. In this figure, the projection of domain D = P 1,D (D) along e 2 leads to an optimal allocation. In this example, the reindexing P 1,D improves the potential parallelism by projection as the projection of the original domain along any direction does not lead to an optimal allocation. The following proposition shows that the whole domain D is located either at the left or at the right hand side of the e l axis.

Proof. By definition, a point z of D is allocated to PE [e t h z ]. Let m = max{e t h z | z ∈ D } -min{e t h z | z ∈ D } + 1.
Clearly, the number of PEs obtained by projecting D along e l is no more than m. Since all of D is to one side of the e l axis (by proposition 4), we know that m = max{|e t h z | | z ∈ D } + 1. Let z be a point of D for which |e t h z | = m -1. Since z, z s ∈ D , and D is convex along e h , Seg(z s , z) ⊂ D . This segment contains m integral points having the same computation date. Hence, the allocation is optimal. 2 Note that projecting D along direction e l amounts to allocating point z ∈ D to the PE numbered |e t h z | = |e t h (zz s,D )|, where z = P h,D (z). In the corresponding array, input and output operations are performed by a PE located at an extremity of the array if all the input/output points belong to F (D, s). This property is suitable for VLSI implementation [START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF][START_REF] Djamegni | Matrix product on modular linear systolic arrays for algorithms with affine schedule[END_REF].

Case 2: No vector of the canonical basis is collinear to λ.

We reduce the problem to the previous case. To do so we apply to the original domain D a reindexing function q that leads to a new timing vector that is parallel to a vector e l of the canonical basis. For this purpose, first note that e l = U -t λ for some unimodular matrix U of order 2 as the columns of e t l and those of λ t generate the same lattice [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]. Constructing matrix U is easy: one only needs to compute the Hermite form [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF] of λ t . From proposition 3 we can set q : z → U z. The number of PEs is minimized by applying to the new index domain D = q(D) the allocation strategy described for case 1. This is illustrated in figure 1.

As stated earlier in case 1, this mapping strategy allocates point z ∈ D to the PE numbered |e Proof. First note that s is the direction of the timing segments of D . This implies that the timing segments of D are of direction s = U -1 s. In addition, the timing segments of D are segments because the image of a segment by an unimodular transformation is also a segment. This means that D is convex along direction s. Now, assume that z = U z. We have z = z s,D + cs and z = z s ,D + c s , where c and c are two negative integers. Because z = U z we get:

z s,D = (c -c)s + U z s ,D (2)
As z s ,D ∈ D and q(D) = D we have U z s ,D ∈ D . Thus, (2) implies that (cc )s + z s,D ∈ D . Because z s,D + s / ∈ D we have (cc ) ≤ 0 as D is convex along direction s. Similarly, as z s,D ∈ D and D = q -1 (D ), [START_REF] Benaini | Space-minimal systolic arrays for Gaussian elimination and the algebraic path problem[END_REF] implies that U -1 ((c -c)s+U z s ,D ) ∈ D. This means that (c -c)s +z s ,D ∈ D. This implies that (cc) ≤ 0 as D is convex along direction s . Therefore c = c , and as a consequence z s,D = U z s ,D . 2

Generalizing to n-dimensional Domains

We now extend the method to higher dimensions. While the discussion above gives an insight on the method, it does not present a constructive approach, especially when we have more than two dimensions. This is done using "slice-based partitions" as described next.

Minimizing the number of PEs along one direction

Given two primitive and linearly independent vectors v 1 = x 1 e 1 +. . .+x n e n and v 2 = y 1 e 1 + . . . + y n e n , we show how the tasks of the iteration domain D can be allocated to PEs so as to minimize the number of PEs along direction v 1 , v 2 , i.e., so as to minimize the number of PEs on all v 1 , v 2 -slice of D. Indeed, as in the partition method [START_REF] Sakho | Méthode de conception d'algorithmes parallèles pour réseaux réguliers[END_REF][START_REF] Tchuenté | Parallel Computation on regular arrays[END_REF], our starting point is to partition D into v 1 , v 2 -slices denoted D 1 , D 2 , . . . , D r . However, we do not search for an allocation function that minimizes the number of PEs on all D i "all at once" but proceed dimension b y dimension. Instead of using ad hoc approaches as in the partition method, we proceed in a systematic way. We perform a preprocessing by reindexing that provides a projection direction that separately allocates the tasks of each D i so as to minimize the number of PEs on each D i .

Note that if v 1 and v 2 are parallel to the timing surfaces the minimal number of PEs required by each v 1 , v 2 -slice is equal to its cardinality. Because of this, we assume that the v 1 , v 2 -slices are not parallel to the timing surfaces, i.e., λ t v 1 = 0 or λ t v 2 = 0. In what follows, we assume that λ t v 2 = 0 without loss of generality. Depending on the direction of the timing vector λ we consider two main cases.

Case 1: λ is collinear to some vector e l of the canonical basis Sub-Case 1.1. The v 1 , v 2 -slices are parallel to some e h such that h = l. Let s be either e h or -e h . We now generalize the reindexing defined by equation [START_REF] Bednara | Interface Synthesis for FPGA Based VLSI Processor Arrays[END_REF]. The main idea is again to shift the timing segments of all D i along direction s so as to place the points of F (D i , s) on the (n -1)hyperplane of Cartesian equation e t h z = 0. To this end, we set:

P h,D (z) = P h,D (z s,D ) + (z -z s,D ) = z s + (e t h (z -z s,D ))e h = (
The reindexing P h,D maintains all points of D on their timing surface and on the supporting 2-hyperplane of its v 1 , v 2 -slice. The points of

F (D, s) = k=r k=1 F (D i , s)
are mapped into points belonging to the (n -1)-hyperplane of Cartesian equation e t h z = 0. This is illustrated in figure 3. The projection of domain D = P 1,D (D) along direction e 3 provides a space optimal array. In this example, the reindexing P 1,D improves the potential parallelism by projection as the projection of the original domain along any direction does not provide a space optimal array. The reindexing P h,D acts similarly on each v 1 , v 2 -slice. It maps the points of F (D i , s) into points belonging to the intersection of the supporting 2-hyperplane of D i and the (n-1)-hyperplane of Cartesian equation e t h z = 0. This implies that P h,D (F (D i , s)) is parallel to vector v 2y h e h . Hence, we have the following: Proposition 7 The projection of D = P h,D (D) along direction v 2y h e h separately allocates the points of all v 1 , v 2 -slice= e h , v 2y h e h -slice so as to minimize the number of PEs on all e h , v 2y h e h -slice.

Proof. The proof is similar to that of proposition 5. 2 Proposition 7 shows that the projection of D = P h,D (D) along direction v 2y h e h leads to an array in which the number of PEs is minimized along direction e h , i.e., for each line of PEs along direction e h there is an instant for which all PEs are simultaneously active. As stated earlier in the case where D is of dimension 2, the resulting array is such that input/output operations are performed by PEs located at the border of the array if all the input/output points belong to F (D, s), a property that is suitable for VLSI implementations.

Sub-Case 1.2. The v 1 , v 2 -slices are not parallel to some e h with h = l.

Here we perform a reindex that reduces the problem to sub-case 1.1. The idea behind this transformation is to locate all D i parallel to some e h with h = l so as to leave unchanged the timing vector λ. To this end, we consider a unimodular vector v 0 that is collinear to vector x l v 2y l v 1 . As the columns of v t 0 and those of e t h generate the same lattice we have e h = U n v 0 for some unimodular matrix U n of order n [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]. Because e t l v 0 = 0 and e t l e h = 0 we can assume that the l -th column and the l -th line of U n are respectively e l and e t l . For instance, if l = n and h = 1 we have

U n = U n-1 0 n-1 0 t n-1 1 ( 4
)
where U n-1 is an unimodular matrix of order n -1 and 0 n-1 is the null vector of Z n-1 . Note that U t n e l = e l . This implies that U -t n e l = e l . Stemming from proposition 3, this implies that the timing vector λ is left unchanged by the reindexing q : z → U n z. Because the potential parallelism is invariant under valid reindexing functions [START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF], the problem of minimizing the number of PEs on each v 1 , v 2 -slice of D can be translated into the the problem of minimizing the number of PEs on each q(v 1 ), q(v 2 ) -slice of D = q(D). This last problem can be solved with the allocation strategy proposed for sub-case 1.1.

Case 2:

No vector of the canonical basis is collinear to λ. We perform a reindexing that reduces the problem to case 1. To do so we consider a reindexing function q that leads to a new timing vector that is parallel to some e l . As e l and λ are unimodular vectors we have e l = U -t n λ for some unimodular matrix U n of order n. From proposition 3 we can set q : z → U n z. Now, we use the allocation strategy of case 1 to minimize the number of PEs on each q(v 1 ), q(v 2 ) -slice of D = q(D).

Space-optimality

The allocation technique presented in section 4.2.1 is based on reindexing transformations that compress the initial iteration domain along a given direction in order to construct a new DAG that leads to an array which is such that the number of PEs is minimized along the given direction. This technique can further be optimized stemming from propositions 8 and 9. The idea behind these propositions is that any array with a 100% pipelining rate whose number of PEs is minimized along any direction is space-optimal. Indeed, this idea can be exploited by successively compressing the initial iteration domain along a number of directions, using reindexing transformations introduced in section 4.2.1. As the new obtained iteration domain is compressed along many directions, it enables the projection method to provide an array whose number of PEs is minimized along many directions. Note that the 100% pipelining rate criterion is easy to satisfy. It is easy to see that if the iteration domain is convex along a projection vector ξ, that satisfies |λ t ξ| = 1, then the pipelining rate of the resulting array is 100%. Recall that the new allocation method applies the projection method with ξ = λ ∈ {e 1 , e 2 , e 3 , . . . , e n }. Therefore, propositions 8 and 9 guarantee the space-optimality of the new allocation technique.

Proposition 8 Any array with a 100% pipelining rate that is such that the intersection of working intervals of any pair of PEs is not empty is spaceoptimal.

Proof. Denote np the number of PEs, pp the potential parallelism and t a date at which the maximum number of PEs are working. Clearly np ≥ pp. Assuming that np > pp, let's denote C = {c 1 , c 2 , . . . , c pp } the set of cells scheduled to work at date t, f (c) (resp. l(c)) the first working date (resp. last working date) of cell c. There is then at least one PE p ∈ C that is not active at date t. Because the array has a 100% pipelining rate, each PE works at all activation date between its first and last working dates. This implies that either t < f (p) or l(p) < t. On the other hand we have

I(c) = [f (c), l(c)] ∩ [f (p), l(p)] = ∅ for all c ∈ C
as the intersection of the working intervals of any pair of PEs is not empty. This implies that f (p) ∈ I(c) for all c if t < f (p) and l(p) ∈ I(c) for all c if t < l(p). It follows that the PEs of C are simultaneously working at date f (p) or at date l(p). Therefore pp + 1 PEs are working at the same date. Hence a contradiction.

2

The proof of the following proposition is simple.

Proposition 9 In any array in which the number of PEs is minimal along all direction, the intersection of activation intervals of every pair of PEs is not empty.

Examples

Here we illustrate the method on Matrix Multiplication (MM) and the Cholesky Factorization (CF). We derive asymptotically space-time optimal arrays for MM and CF.

Matrix multiplication [START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF][START_REF] Djamegni | Exécution d'un graphe cubique de tâches sur un réseau 2D et asymptotiquement optimal[END_REF] The computation of the product C = AB of two square matrices A and B of order N can be defined by the following SURE:

Initialization: For 0 ≤ i ≤ N -1, 0 ≤ j ≤ N -1, 0 ≤ k ≤ N -1 A(i, -1, k) = A i,k B(-1, j, k) = B k,j C(i, j, -1) = 0 Computation: For 0 ≤ i ≤ N -1, 0 ≤ j ≤ N -1, 0 ≤ k ≤ N -1 A(i, j, k) = A(i, j -1, k) B(i, j, k) = B(i -1, j, k) C(i, j, k) = C(i, j, k -1) + A(i, j -1, k)B(i -1, j, k) Output: For 0 ≤ i ≤ N -1, 0 ≤ j ≤ N -1 C i,j = C(i, j, N -1)
The iteration space is

D (0) = {(i, j, k) ∈ Z 3 | 0 ≤ i ≤ N -1, 0 ≤ j ≤ N -1, 0 ≤ k ≤ N -1}.
The dependence vectors are the vectors of the canonical basis of Z 3 : e 1 , e 2 and e 3 . A corresponding optimal timing function 20

(0) D D (1) D (2) e e 1 3 1,D P q e 2 D (3)
P (2) ( Regarding the dependencies of (8) a uniform version can be obtained with e 1 , e 2 and e 3 as the dependence vectors [START_REF] Quinton | The mapping of linear equations on regular arrays[END_REF], and this without changing the domain D (0) of equations ( 8). An optimal timing function corresponding to such a uniformization is t(i, j, k) = i + j + k. By projecting the domain D (0) along vector e 1 + e 2 , we obtain a triangular orthogonally connected array of N (N + 1)/2 PEs. The number of PEs can be reduced as follows. First, the timing vector e 1 + e 2 + e 3 is mapped into e 1 with the reindexing q(i, j, k) (2) ) [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF]. Third, a asymptotically spacetime optimal array of N 2 8 + Θ(N ) [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF] PEs is obtained by projecting D (3) along vector e 1 . This improves the results obtained from the instruction shift method [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF] and the partition method [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF]. They lead to arrays of

= (i + j + k, j, k). Second, D (1) = q(D (0) ) ⊆ {(i, j, k) ∈ Z 3 | - i + 2j + k ≤ 0, i -j -k -N ≤ 0, -k ≤ 0, -j + k ≤ 0} is mapped into D (2) = P 3,D (1) (D (1) ) ⊆ {(i, j, k) ∈ Z 3 | j + k -N ≤ 0, -k ≤ 0, i - 2j + k -N ≤ 0, -i + 2j + k ≤ 0, -j + k ≤ 0 } and D (2) is in turn mapped into D (3) = P 2,D (2) (D
N 2
6 + Θ(N ) PEs. Note that the projection of D (2) along vector e 1 leads to an array of N 2 4 + Θ(N ) [START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF] PEs. Therefore the numbers of PEs is reduced by a factor of two at each of the last two steps of the reindexing process. Stemming from the results of section 7 we have:

P 3,D (1) (i, j, k) = (i, j, -i + j + k + N ) if (i, j, k) ∈ D (1) 1 (i, j, k) if (i, j, k) ∈ D (1) 2
where D

(1)

1 = {(i, j, k) ∈ D (1) | -i + j + N < 0} D (1) 2 = {(i, j, k) ∈ D (1) | -i + j + N ≥ 0} P 2,D (2) =        (i, -1 2 i + j -1 2 k + N 2 , k) if (i, j, k) ∈ D (2) 1,0 (i, -1 2 i + j -1 2 k + N +1 2 , k) if (i, j, k) ∈ D (2) 1,1 (i, j -k, k) if (i, j, k) ∈ D (2) 2
where

D (2) 1,0 = {(i, j, k) ∈ D (2) | -i + j + N < 0 ∧ (i + k -N )(mod 2) = 0} D (2) 1,1 = {(i, j, k) ∈ D (2) | -i + j + N < 0 ∧ (i + k -N )(mod 2) = 1} D (2) 2 = {(i, j, k) ∈ D (2) | -i + j + N ≥ 0}
6 Comparison with other methods

In this section, we compare our allocation method with the projection method, the grouping (or clustering) method, the instruction shift method and the partition method. Our method is sustained by a space-optimality result (section 4.2.2) while the methods mentioned above do not guarantee space-optimality. The space-optimality is achieved by compressing the initial iteration domain along all directions with reindexing transformations introduced in section 4.2.1. In what follows, we establish that the lowest space complexity that both the projection method and the grouping method can give for the initial iteration domain is an upper bound on the space complexity of our method, i.e., our method can do no worse.

Preliminary result

Let λ be a timing vector and let v 1 = x 1 e 1 + . . . + x n e n and v 2 = y 1 e 1 + . . . + y n e n ,v 1 be two unimodular linearly independent vectors such that v 1 , v 2 -slices are not parallel to timing surfaces.

We have the following result.

Proposition 10 If λ is collinear to some e l and v 1 , v 2 -slices are parallel to some e h with l = h then N P E(v 2y h e h , P h,D (D))

≤ N P E(v 2 , D, e h ) ≤ N P E(v 2 , D).
Proof. Because the v 1 , v 2 -slices are not parallel to timing surfaces, i.e., λ t v 1 = 0 or λ t v 2 = 0, we can assume without loss of generality that λ t v 2 = 0. It is easy to see that N P E(v 2 , D, e h ) ≤ N P E(v 2 , D). N P E(v 2 , D) = i=r i=1 N P E(v 2 , D i ) as points allocated to a same PE belong to a same line of direction v 2 , and obviously to a same v 1 , v 2 -slice. Denote ξ 0 = v 2y h e h and D i = P h,D (D i )). We have N P E(v 2 , D, e h ) = i=r i=1 N P E(v 2 , D i , e h ) and N P E(ξ 0 , D ) = i=r i=1 N P E(ξ 0 , D ). On the other hand, the reindexing P h,D retains the potential parallelism on each v 1 , v 2 -slice [START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF], i.e., the potential parallelism of D i is equal to that of D i . In addition, according to proposition 7, the potential parallelism of D i is equal to N P E(ξ 0 , D i ). This implies that N P E(ξ 0 , D i ) ≤ N P E(v 2 , D i ) and N P E(ξ 0 , D i ) ≤ N P E(v 2 , D i , e h ). Therefore N P E(ξ 0 , D ) ≤ N P E(v 2 , D, e h ) ≤ N P E(v 2 , D). 2

Two strong comparison results

Here, we establish two strong comparison results. We first establish a strong comparison result between our allocation method and the projection method. Then, we establish another strong comparison result between our allocation method and the clustering method.

To simplify the presentation, we assume without loss of generality that the timing vector λ = e n . Thus t = [e n , α].

Denote

D (0) = D, D (h) = P h,D (h-1) (D (h-1)
) with h ∈ {1, 2, . . . , n -1} and assume that D (h-1) is convex along direction e h . Consider an unimodular vector ξ = ξ 1 e 1 + ξ 2 e 2 + . . . + ξ n e n that corresponds to a valid projection vector. As λ = e n , we get ξ n = 0. Denote ξ (h) , h ∈ {0, 1, 2, . . . , n -1}, an unimodular vector that is collinear to vector ξ h+1 e h+1 +ξ h+2 e h+2 + . . .+ ξ n e n .

Corollary 1 We have N P E(D (n-1) , λ) ≤ N P E(D (0) , ξ) for all vector ξ such that λ t ξ = 0.

Proof. From proposition 10 we get

N P E(ξ (h) , D (h) ) ≤ N P E(ξ (h-1) , D (h-1) ), h ∈ {1, 2, 3, . . . , n -1}
This implies that N P E(e n , D (n-1) ) ≤ N P E(ξ, D (0) ). 2

As stated earlier, Wong and Delosme [START_REF] Wong | Space-optimal linear processor allocation for systolic arrays synthesis[END_REF], Shang and Fortes [START_REF] Shang | On mapping of uniform dependence algorithms into lower dimensional processors arrays[END_REF], Ganapathy and Wah [START_REF] Ganapathy | Optimal synthesis of algorithm-specific lower-dimensional processor arrays[END_REF], have developed an approach based on integer programming to get the best projection direction possible. But for large problem instances the time spent to solve such programs cannot be neglected [START_REF] Darte | Synthesizing systolic arrays: Some recent developments[END_REF]. Here we propose another approach based on reindexing transformations that permits to do better. A straightforward consequence of corollary 1 is that it provides a procedure based on a preprocessing by reindexing that systematically leads to arrays whose number of PEs is bounded from above by the number of PEs obtained from the best projection direction. We also have the following.

Theorem 1 We have N P E(D (n-1) , λ) ≤ 1 |λ t ξ| N P E(D (0) , ξ)+ h=n-2 h=0
Θ(nl h ) for all vector ξ such that λ t ξ = 0 where nl h denotes the number of distinct e h+1 , ξ (h) -slices of D (h) .

Proof. From the definition of ξ (h) , h ∈ {0, 1, 2, . . . , n -2}, there exists an integer number l h such that e t n ξ (h) = e t n ξ (h+1) l h . Denote A h the array obtained by projecting domain D (h) along direction ξ (h) . In this array, all PE works only once over |λ t ξ (h) | = |e t n ξ (h) | successive activation dates. Now, let's examine the working dates of PEs belonging to a same line of direction e h+1 denoted L h . All the tasks allocated to L h belong to a same ξ (h) , e h+1 -slice= e h+1 , ξ (h+1) -slice of D (h) . The computation dates d 1 and d 2 of two arbitrary tasks belonging to a same e h+1 , ξ (h+1) -slice are such that (d 1d 2 ) mod |e t n ξ (h+1) | = 0. This implies that there are |l h | activation dates of L h over |e t n ξ (h) | successive activation dates of A h as e t n ξ (h) = e t n ξ (h+1) l h . It follows that all PE of L h works only once over |l h | successive activation dates of L h . Thus by applying the clustering technique to any line of PEs parallel to vector e h+1 , the space complexity of array A h is reduced by a factor of |l h |, i.e.,

N P E(ξ (h) , D (h) , e h+1 ) = 1 |l h | N P E(ξ (h) , D (h) ) + Θ(nl h )
On the other hand from proposition 10 we get

N P E(ξ (h+1) , D (h+1) ) ≤ N P E(ξ (h) , D (h) , e h+1 )
for all h ∈ {0, 2, 3, . . . , n -2}. This implies that

N P E(ξ (h+1) , D (h+1) ) ≤ 1 |l h | N P E(ξ (h) , D (h) ) + Θ(nl h )
for all h ∈ {0, 1, 3, . . . , n -2}. It follows that

N P E(ξ (h+1) , D (h+1) ) ≤ 1 |l 0 l 1 l 2 . . . l h | N P E(ξ (0) , D (0) ) + i=h i=0 Θ(nl i )
for all h ∈ {0, 1, 3, . . . , n -2}. Moreover for h = n -2 we have |l 0 l 1 l 2 . . . l h | = |λ t ξ|. Hence the result. 2

A straightforward consequence of theorem 1 is that the new allocation method can be applied so as to systematically derive arrays whose space complexities are bounded from above by the lowest space complexities that the grouping method can give. Compared to other allocation methods, we are not aware of any similar strong comparison result between the instruction shift method and the grouping method, nor between the partition method and the grouping method. No similar strong comparison result exists even between the instruction shift method and the projection method, nor between partition method and the projection method.

Closed form of the allocation function

We now investigate the implementation of the new allocation method in tools for (semi) automatic derivation of VLSI arrays. In particular, we consider the MMAlpha [START_REF] Guillou | Hardware Design Methodology with the Alpha Language, Forum on Design Languages[END_REF] system based on the functional language Alpha [START_REF] Mauras | ALPHA: Un langage équationnel pour la conception et la programmation d'architectures parallèles et synchrones[END_REF]. Alpha programs are almost identical to the systems of equations that we have used in this paper. Each equation is defined over a possibly parameterized index domain that is a finite union of polyhedra [START_REF] Quinton | Extension of the ALPHA language to recurrences on sparse periodic domains[END_REF]. Such domains constitute an abstract data type with closure under certain operations, notably intersection, union, and image by bijective affine index transformations. Such an index transformation is valid if every point of the index domain to which it is applied has a unique integral image. In particular, any unimodular transformation is always a valid index transformation.

In order to implement our method in MMAlpha, we only need to show how to implement allocation function (3) since it is the only function that may not correspond to an obvious unimodular transformation. The function, repeated below for convenience, is a piecewise affine function, but specified implicitly.

P h,D (z) = (e t 1 z, . . . , e t h-1 z, e t h (zz s,D ), e t h+1 z, . . . , e t n z) t

Let s denote either e h or -e h , where e h is the h -th vector of the canonical basis. To simplify the presentation, we assume that D = Z n P where P is a single convex1 polytope, Q n defined by an irredundant set of integral inequalities {z ∈ Q n | Az ≤ b}. Therefore, any row A i of matrix A defines a facet F i = {z ∈ P | A i zb i = 0} [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF] of D. Among the F i 's, we choose those that satisfy A i s > 0. We assume without loss of generality that these facets are F 1 , F 2 , F 3 , . . . , F r for some integer r. Proposition 11 shows that such facets exist.

Proposition 11 There exists an integer i such that A i s > 0.

Proof. Suppose that A i s ≤ 0 for all i. This implies that A i (z + αs) ≤ b i for any positive rational number α and for any z belonging to P . This implies that z + αs ∈ P , which in turn implies that s is a ray of P . This contradicts the fact that P is bounded. 2

Step 1: We first partition D into dense sub-domains [START_REF] Quinton | Extension of the ALPHA language to recurrences on sparse periodic domains[END_REF]. Let σ i (z) denote the projection of point z on the supporting hyperplane of facet F i along direction s. Let P i = {z ∈ P | σ i (z) ∈ F i } and D i = Z n P i for i = 1, 2, . . . , r. The following proposition provides us a constructive method to determine the dense sub-domains. for some i. We have A(z + α z s) ≤ b and A i (z + α z s)b i = 0. This implies that σ i (z) = z + α z s ∈ F i , which in turn implies that z ∈ P i . Now consider two points z 1 and z 2 of a given P i . Let z 3 ∈ Seg(z 1 , z 2 ). Because P is convex we have z 3 ∈ P . In addition, a simple linear algebraic argument shows that σ i (z 3 ) ∈ Seg(σ i (z 1 ), σ i (z 2 )), and this implies that σ i (z 3 ) ∈ F i because F i is convex. Thus z 3 ∈ P i , and Seg(z 1 , z 2 ) ⊆ P i . Therefore P i is convex. 2

Step 2: Next, we partition the D i 's into sparse periodic domains.

Proposition 13 Let z be a point of P i and let α be a positive rational number. We have:

z + αs ∈ P i iff A i (z + αs) ≤ b i (9)
Proof. Let us assume that A i (z + αs) ≤ b i . From the definition of σ i we get σ i (z) = z + b i -A i s A i s s. This implies that z + αs ∈ Seg(z, σ i (z)) as 0 ≤ α ≤ b i -A i z A i s . We have Seg(z, σ i (z)) ⊆ P as domain P is convex along direction s. Thus z + αs ∈ P . In addition, σ i (z + αs) = σ i (z). Therefore z + αs ∈ P i . The necessary condition is easy to check. 2 Let z ∈ D i . Clearly, z s,D = z + αs for some integer number α. We have A i (z + αs) ≤ b i . From proposition 13, we get A i (z + (α + 1)s) > b i as z s,D + s / ∈ P i . Thus

b i -A i z A i s -1 < α ≤ b i -A i z A i s .
This implies that

α = b i -A i z A i s , which means that α = b i -A i z -j A i s (10)
where (b i -A i z) mod A i s = j.

Stemming from [START_REF] Djamegni | Contribution to the synthesis of optimal algorithms for regular arrays[END_REF], we consider a partition D i,j = {z ∈ D i | (b i -A i z) mod A i s = j}, j = 0, 1, 2, . . . , A i s -1 of D i . D i,j is a sparse domain if j = 1.

Step 3: Analytic expression of the re-indexing function P h,D Assuming that z ∈ D i,j , we have from (3)

P h,D (z) =
For h = 1 this is equivalent to

P h,D (z) = 1 B i 0 n-1 I n-1
z + (j-b i )(e t 1 s) A i s 0 n-1 where 0 n-1 is the null vector of Z n-1 , I n-1 the identity matrix of order n -1 and e t 1 s A i s A i = 1 B i .

This shows that the restriction of P h,D to each D i,j corresponds to an affine transformation. If A i s = 1 we have D i = D i,0 and the restriction of P h,D to D i corresponds to an unimodular transformation. If A i s > 1 then the D i,j 's are sparse periodic domains [START_REF] Quinton | Extension of the ALPHA language to recurrences on sparse periodic domains[END_REF] and the reindexing P h,D corresponds to a piecewise affine transformation. In this particular case, the resulting array may correspond to a piecewise regular array [START_REF] Teich | Partitioning of processor arrays: a piecewise regular approach[END_REF].

Conclusion

In this paper, we have formally introduced a new allocation method that is based on a preprocessing by reindexing that leads to a new index domain that permits to systematically derive arrays in which the number of PEs is minimized along a number of directions. In addition, the preprocessing permits to improve the potential parallelism by projection of the initial domain. It also permits to allocate input/output nodes to PEs located at the border of the array. This can permit to eliminate eventual additional steps for data loading and result unloading that may force a slowdown of the algorithm [START_REF] Djamegni | Derivation of systolic algorithms for the algebraic path problem by recurrence transformations[END_REF][START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Mapping rectangular mesh algorithms onto asymptotically space-optimal arrays[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Djamegni | Matrix product on modular linear systolic arrays for algorithms with affine schedule[END_REF]. This new allocation approach is simple, suitable for software tools and could be included in a language manipulating recurrence equations like the ALPHA language.

Figure 1 :

 1 Figure 1: Illustration of the allocation method on the bubble sort

Figure 2 :

 2 Figure 2: Illustration of the reindexing P 1,D . The timing vector is e 2

Proposition 4 2 Proposition 5

 425 For all z ∈ D, either e t h (zz s,D ) ≥ 0 or e t h (zz s,D ) ≤ 0.Proof. Assume that there exist two points z and z of D such that e t h (zz s,D ) > 0 and e t h (zz s,D ) < 0. This implies that z = z s,D + ce h and z = z s,Dc e h where c and c are two strictly positive integers. Because D is convex along e h , segments Seg(z, z s,D ) and Seg(z , z s,D ) are in D. This implies that z s,D + e h and z s,De h belong to D. This is in contradiction either with the definition of z s,D or of z s,D as s = e h or s = -e h . The projection of D along e l leads to an optimal allocation.

Proposition 6

 6 t h (zz s,D )| where s denotes either e h or -e h with h = l. Proposition 6 shows that this amounts to allocate point z ∈ D to the PE numbered |e t h U (zz s ,D )| with s = U -1 s. Given two points z, z ∈ D. If z = U z then z s,D = U z s ,D where s = U -1 s.

Figure 3 :

 3 Figure 3: Illustration of the reindexing P 1,D on a three-dimensional domain. The timing vector is e 3

Figure 4 :

 4 Figure 4: Reindexing process for matrix multiplication

  each P k is a convex polytope. As a result, D = k=r k=1 D k . Proof. Let z ∈ P . Denote α z = min{ b k -A k z A k s | k = 1, 2, . . . , r}, i.e., α z = b i -A i z A i s

  e t 1 z, . . . , e t h-1 z,

	e t h s A i s	(-b

i + A i z + j), e t h+1 z, . . . , e t n z t

In the general case, we would require a preprocessing step to first decompose the domain into a number of sub-domains.

is t(i, j, k) = i + j + k. The projection of D (0) along vector e 3 leads to an N × N array of N 2 PEs. To reduce the number of PEs we proceed as follows. First, the timing vector e 1 + e 2 + e 3 is mapped into e 3 with the reindexing q: q(i, j, k) = (i, j, t(i, j, k) = (i, j, i + j + k) [START_REF] Bu | Processor clustering for the design of optimal fixed-size systolic arrays[END_REF] This leads to a new iteration domain D (1) = q(D (0 [START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF][START_REF] Djamegni | Exécution d'un graphe cubique de tâches sur un réseau 2D et asymptotiquement optimal[END_REF]. This is illustrated in figure 4. Second, D (1) is compressed along direction e 1 with the reindexing P 1,D (1) to obtain [START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF][START_REF] Djamegni | Exécution d'un graphe cubique de tâches sur un réseau 2D et asymptotiquement optimal[END_REF]. Similarly, D (2) is compressed along direction e 2 with the reindexing P 2,D (2) to obtain [START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF][START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF][START_REF] Djamegni | Exécution d'un graphe cubique de tâches sur un réseau 2D et asymptotiquement optimal[END_REF]. Stemming from the results of section 7 we have: (3) is projected along direction e 3 to obtain an asymptotically spacetime optimal array of 3 4 N 2 + Θ(N ) [START_REF] Djamegni | Derivation of 2D space-optimal regular arrays for 3D cubical mesh algorithms[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF] PEs.

Cholesky factorization [START_REF] Clauss | Optimal mapping of systolic algorithms by regular instructions shifts[END_REF][START_REF] Djamegni | Synthesis of space-time optimal systolic algorithms for the Cholesky factorization[END_REF][START_REF] Djamegni | Méthodes d'optimisation pour la synthèse d'architectures régulières[END_REF] The CF is defined as follows: Given a N × N symmetric positive definite matrix A, the CF calculates a lower triangular matrix L such that A = LL t . It is defined by the following well known affine recurrence equations: