Array-OL is a high-level specification language dedicated to the definition
of intensive signal processing applications. Several tools exist for
implementing an Array-OL specification as a data parallel program. While
Array-OL can be used directly, it is often convenient to be able to deduce part
of the specification from a sequential version of the application. This paper
proposes such an analysis and examines its feasibility and its limits

Feautrier, Paul

English

arXiv

Array-OL is a high-level speciﬁcation language dedicated to the definitionof intensive signal processing applications. Several tools existfor implementing an Array-OL speciﬁcation as a data parallel program.While Array-OL can be used directly, it is often convenient to be able todeduce part of the speciﬁcation from a sequential version of the application.This paper proposes such an analysis and examines its feasibilityand its limits.Array-OL est un système de spéciﬁcation de haut niveau spécialisé dans la déﬁnition d’application de traitement du signal intensif. Il existe plusieurs ateliers qui transforment une spéciﬁcation Array-OL en un programme à parallélisme de données. Bien que Array-OL puisse être utilisé tel quel, il est souvent intéressant de pouvoir déduire ses paramètres d’une version séquentielle de l’application. Ce rapport propose une telle analyse et en examine la faisabilité et les limites

Hal-Diderot

Elementary transformation analysis for array-OL

HAL Descartes

HAL-ENS-LYON

Array-OL is a high-level specification language dedicated to the definition of intensive signal processing applications. Several tools exist for implementing an Array-OL specification as a data parallel program. While Array-OL can be used directly, it is often convenient to be able to deduce part of the specification from a sequential version of the application. This paper proposes such an analysis and examines its feasibility and its limits

Elementary transformation analysis for Array-OL

Archive Ouverte en Sciences de l'Information et de la Communication

HAL Id: hal-02102502
https://hal-lara.archives-ouvertes.fr/hal-02102502
Submitted on 17 Apr 2019
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Elementary transformation analysis for array-OL
Paul Feautrier
To cite this version:
Paul Feautrier. Elementary transformation analysis for array-OL. [Research Report] LIP RR-2007-28,
Laboratoire de l’informatique du parallélisme. 2007, 2+7p. ￿hal-02102502￿
Laboratoire de l’Informatique du Parallélisme
École Normale Supérieure de Lyon
Unité Mixte de Recherche CNRS-INRIA-ENS LYON-UCBL no 5668
Elementary transformation analysis for
Array-OL
Paul Feautrier May 2007
Research Report No 2007-28
École Normale Supérieure de Lyon
46 Allée d’Italie, 69364 Lyon Cedex 07, France
Téléphone : +33(0)4.72.72.80.37
Télécopieur : +33(0)4.72.72.80.80
Adresse électronique : lip@ens-lyon.fr
Elementary transformation analysis for Array-OL
Paul Feautrier
May 2007
Abstract
Array-OL is a high-level specification language dedicated to the def-
inition of intensive signal processing applications. Several tools exist
for implementing an Array-OL specification as a data parallel program.
While Array-OL can be used directly, it is often convenient to be able to
deduce part of the specification from a sequential version of the applica-
tion. This paper proposes such an analysis and examines its feasibility
and its limits.
Keywords: Array-OL, multidimensional signal processing, program analysis
Résumé
Array-OL est un système de spécification de haut niveau spécialisé dans
la définition d’application de traitement du signal intensif. Il existe plu-
sieurs ateliers qui transforment une spécification Array-OL en un pro-
gramme à parallélisme de données. Bien que Array-OL puisse être utlisé
tel quel, il est souvent intéressant de pouvoir déduire ses paramètres
d’une version séquentielle de l’application. Ce rapport propose une telle
analyse et en examine la faisabilité et les limites.
Mots-clés: Array-OL, traitement du signal multidimensionnel, analyse de programme
Analysis for Array-OL 1
1 Introduction
In the Array-OL formalism [1, 2], a program is a network of processes which communicate
through shared arrays. A process is made of one or more parallel loops. At each iteration
of these loops, a task (or elementary transform) is executed. The elementary transform may
contain one or more loops, which are executed sequentially.
The execution of an elementary task can be decomposed into three steps:
• Move portions of the input array(s) (regions) to the local memory of the processor
executing the task.
• Execute the elementary transform and generate portions of the output array(s).
• Move the results to the output array(s).
In order to simplify code generation, the input and output regions must move uniformly
across the shared arrays. It is admissible that each elementary transform use only a subset
of regularly spaced entries in the input and output regions. In the present version of the
software, regions must not overlap, as this would precludes parallel execution of the outer
loops. The useful elements of a region are collected in a pattern, which must be a rectangular
parallelepiped of fixed size.
The Array-OL formalism may be used directly. The programmer is responsible for con-
structing the elementary transform, identifying the input and output regions, checking par-
allelism and specifying the regions parameters. Another possibility is to infer the Array-OL
specification from a sequential version of the program. This requires the solution of three
problems:
• Rewriting the sequential program in such a way that the outer loops have no depen-
dences.
• Deducing the shape and size of the regions from an analysis of the array subscript
functions.
• Rewriting the sequential code by substituting pattern accesses to the original array
accesses.
This note is dedicated to a proposal for the solution of the second and third problems. The
assumption is that one is given the sequential code, together with a list of input and output
arrays, and an indication of which loop(s) are to be considered as the outer (repetition)
loop(s).
2 Paving
Let A be an input or output array and let its occurences in the sequential code be numbered
from 1 to N . Let r be the counter(s) of the repetition loop(s), and let jk be the counter(s) of
the inner loop(s) that surround occurence k of A. Let ek(r, jk) be its subscript function. ek
is a vector function whose dimension is the rank of A.
To be amenable to an Array-OL implementation, the subscript function ek must be affine
in r and jk. A convenient way of checking this property consists in computing the two
Jacobian matrices:
2 P. Feautrier
P k = (
∂ekα
∂rβ
) Bk = (
∂ekα
∂jkβ
),
checking that they do not depend on r or jk, and verifying the identity:
ek(r, jk) = P kr + Bkjk + ek(0, 0).
In Array-OL terminology, P k is the paving matrix, and ek(0, 0) is the origin of the paving.
The elements of these entities may be numbers, or they may depend on constants, which
must be given numerical values just before code generation. References with different paving
matrices may be separated by arbitrary distance in the source or target array; it is not possible
to group them efficiently; they must be implemented as separate channels.
In the following example:
myTE( in[][], out[]){
for(i=0;i<7; i++) // boucle TE
{
for ( k=0;k<11;k++)
{
S=0;
for(j=0;j<100;j++)
{
S+= in[0][j+11] * in[i+1][k+j];
}
out[ i][k]=S;
}
}
there are two references to in with repective subscript functions e1(i, k, j) =
(
0
j + 11
)
and e2(i, k, j) =
(
i + 1
k + j
)
. The corresponding paving matrices are P 1 =
(
0
0
)
and
P 2 =
(
1
0
)
. Hence, the two accesses must be handled separately.
In the following, I assume that accesses to A have been partitioned according to their
paving matrix, and consider only one partition at a time. The size of the repetition space
is deduced simply from the bound(s) of the elementary transform loop(s). In the Spear/DE
implementation of Array-OL, there may be further constraints on the paving matrix (e.g.
that it be a permutation of a diagonal matrix).
3 Pattern and fitting
A pattern is a compact specification of all the elements of an array that are accessed, with
references having the same paving matrix, in one iteration of the external loop(s).
When discussiong patterns, one has to consider three frames of reference (see Fig. 1).
The first one is the original (input or output) array. Its dimension is the rank of the array,
Analysis for Array-OL 3
array
iteration
domain
pattern
footprint
local
subscript
function
fitting
modified subscript
function
Figure 1: Data access in Array-OL
noted |A|, and its coordinates are called subscripts. The shape of an array is always a (hyper-)
rectangle.
The second frame of reference is the iteration space of the inner loops of the elementary
transform. Its dimension is the number of loops enclosing the reference, noted dk, and its
coordinates are called loop counters. There may be as many iteration domains as there
are references, or several references may share the same iteration domain. The shape of an
iteration domain is arbitrary. The only requirement in the present context is to be able to
construct its vertices, either because the iteration domain is rectangular, or because it can be
expressed as a convex polyhedron with parameters in the constant terms only. The iteration
domain of reference k will be denoted as Dk in what follows.
The third frame of reference is the pattern. According to Boulet [1] the pattern is always of
rectangular shape. The pattern associated to reference k is denoted by T k and its dimension
is pk. The associated fitting matrix, F k, connects the pattern space to the array space and
its dimension, accordingly, is |A| × pk.
The relation of these objects are as follows. Firstly, the local subscript function fk(jk) =
Bkjk +ek(0, 0) = ek(0, jk) gives the coordinates of an array cell relative to the reference point
P k.r which moves according to the paving matrix.
Next, the image fk(Dk) is the footprint of reference k. Its shape is arbitrary. The images
of the vertices of Dk by fk form a superset of the vertices of the footprint; a representation
as a convex polyhedron can be recovered by one application of the Chernikova algorithm [3].
Lastly, the image of the pattern by the fitting matrix must enclose the footprint, and it
must be feasible to retrieve a datum from the pattern instead of the original array. This
implies that there exists a function φk from Dk to T k such that for every iteration vector
jk ∈ Dk, fk(jk) = F kφk(jk). In the text of the elementary transform, φk must be substituted
to ek in reference k to A.
As one may see from this discussion, while the iteration domain and footprint are fixed once
the sequential program is given, the choice of the pattern and fitting matrix are somewhat
4 P. Feautrier
arbitrary. There are two obvious solutions: in the first one, the pattern is the smallest
rectangular box enclosing the footprint, the fitting matrix is the identity, and the subscript
function is not changed. In the second solution, the pattern is isomorphic to the iteration
domain (provided it is a parallelepiped), Bk is the fitting matrix, and the new subscript
function is the identity.
In signal processing applications, it is often the case that several references to the same
array have similar subscript functions; constructing only one pattern for several references
is an interesting optimization. However, this should not be obtained at the cost of a large
overhead in the size of the pattern. In other word, the number of useless elements in the
pattern must be minimized. Useless elements come from two sources:
• A subscript matrix which is not of full row rank: the pattern will have more dimensions
than the footprint.
• A subscript matrix whose determinant is not of modulus one: there will be holes (un-
used elements) in the footprint. The inverse of the determinant gives an asymptotic
evaluation of the ratio of useful elements.
The next section presents a method for computing a pattern and a fitting matrix in the
general case (many references). This method can only be applied if all elements of the matrices
Bk and the vectors bk have known numerical values. Section 5 presents fail-soft solutions for
cases in which these elements depend on unknown parameters.
4 The General Case
The basic observation is that a conservative estimate of the footprint can be obtained by
computing the projection of each iteration domain by the associated subscript function, then
constructing a convenient superset of the union of these projections. One practical method
consists in projecting the vertices of the iteration domains. One then gathers all such projec-
tions, and constructs their convex hull by familiar (e.g., Chernikova’s) algorithms.
To reduce the size overhead, one should notice that a useful point for reference k also
belongs to the lattice which is generated by the column vectors of Bk. Hence, Bk, properly
simplified (see later) could be used as the fitting matrix. However, in the case of several
references, we have to combine several lattices into one, since each pattern has only one
fitting matrix. As an illustration of this construction, consider the one dimensional case.
A one-dimensional lattice is simply a set of regularly spaced points. Combining two lattices
generates a lattice whose spacing is the gcd of the component spacings. The many-dimensional
equivalent of the gcd is the construction of the Hermite normal form of the subscript matrices.
Let Λ(B, b) be the lattice generated by B with origin b, i.e. the set of points {Bx+ b | x ∈
INd}. Let L1 = Λ(B1, b1) and L2 = Λ(B2, b2) be two such lattices. I claim that the union of
L1 and L2 is included in the lattice L = Λ([B1B2(b2 − b1)], b1).
Proof Let B1.x + b1 be a point of L1. We have:
B1.x + b1 = B1.x + B2.0 + (b2 − b1).0 + b1
hence B1.x + b1 is in L. Similarly:
B2.y + b2 = B1.0 + B2.y + (b2 − b1).1 + b1.
Analysis for Array-OL 5
I conjecture that L is the smallest lattice which includes L1 and L2. The proof
is obvious if the bs are null. The general case is left for future work.
The construction can be extended to any number of component lattices. The resulting
matrix is [B1 . . . BN (b2 − b1) . . . (bN − b1)] and the origin is b1. Furthermore, b1 can be moved
to the origin of the paving and hence taken as 0 when computing the fitting.
In case where B has been obtained by mixing many references, it must be simplified before
being used for an Array-OL specification.
The starting point of this simplification is the row echelon form of B. One can show (see
the appendix) that there exists two unitary matrices P and U such that:
B = P
[
H 0
C 0
]
U,
where H is a square upper triangular matrix of size r × r with positive diagonal coefficients,
C is arbitrary, and both 0 represent null matrices of appropriate sizes. r is the row rank of
B. Furthermore, U can be partitioned, row wise, in two matrices of size r×d and (d− r)×d,
U =
[
U ′
U ′′
]
.
Let j be a point in the iteration domain of the inner loops. The corresponding point in
the footprint is:
Bj = P
[
H 0
C 0
] [
U ′
U ′′
]
j (1)
= P
[
H
C
]
(U ′j) (2)
One possible interpretation of this formula is that the pattern for the current reference is the
image of its iteration domain by U ′, and that the corresponding paving matrix is P
[
H
C
]
.
In the body of the elementary transform, accesses to Bj in the input or output array have to
be replaced by accesses to U ′j in the pattern. It may be that the pattern computed in this
way is not rectangular, in which case it must be “boxed” by computing the component-wise
minima and maxima of its extreme points. The dimension of the pattern is r.
It is interesting to notice that this general solution reduces to one of the approximate
methods above in special cases. If B is unitary, then its row echelon form is the unit matrix.
In that case, the pattern is the footprint, eventually extended to a rectangular box and the
fitting matrix is the identity. Conversely, if B is already in row echelon form, P and U are
identities. The pattern is isomorphic to the iteration space, and B is the fitting matrix.
5 The Parametric Case
Parameters occurs mostly in loop bounds. They may also appear as strides and, more seldom,
in the coefficients of subscript functions.
In the Array-OL formalism, the repetition loops must be square. Hence, their bound may
be extracted diretcly from the program text. The extraction of the paving matrix is a simple
derivative computation, which is an easy task for a competent computer algebra system.
6 P. Feautrier
Similarly, the Bk matrices are the result of a derivation, and may contain parameters.
There are no restrictions on the inner loops. For the construction of the pattern, one
needs to know the vertices of the inner iteration domain. There are three cases:
• The bounds are constant: they can be extracted even if parametric.
• The bounds are affine expressions in other loop counters and parameters: the vertices
can be computed with the help of the polylib.
• In other cases, there is no way of computing vertices, but the user may supply a bounding
box.
The computation of the row echelon form can be done only if the matrix is known numer-
ically, except in two cases: the matrix is 1 × 1 (it is its own normal form) or 2 × 2.
The row echelon form of
(
a b
c d
)
is
(
gcd(a, b) 0
cu + dv |(ad − bc)|/ gcd(a, b)
)
where u et v
are the integers such that au + bv = gcd(a, b) whose existence is guaranteed by Bezout
identity.
If none of these circumstance applies, the solution of last resort is to use one of the
approximate schemes above. For instance, if the vertices of the inner iteration domain are
available, it is possible, whatever the B matrix, to compute the vertices of the footprints and
to enclose them in a rectangular box. The paving matrix is then the identity.
6 Extensions
The Syntol tool computes dependences; it is thus possible to check that the repetition loops
are actually parallel. One must take care that Syntol will find dependences if temporary
scalars are used in the code of the elementary transforms. These scalars must be expanded
or privatized at code generation time.
Overlap between patterns (or, rather, between footprints) is another concern. For input
arrays, overlap is just a cause of inefficiency, since some arrays cells will be copied several
times to processors. Overlap for output arrays are more dangerous since they may induce
non-determinism. The existence of overlap may be tested provided one stays inside the
polytope model (affine loop bounds and indexing functions, with numerical coefficients and
linear parameters). In the same context, it is possible to quantify the overhead by comparing
the size of the pattern and the size of the real footprint using the barvinok library [4].
A Computing the row echelon form of a matrix
For more details, see [3]. Let B be an arbitrary matrix of size p × q.
1. At any stage of the computation, we have constructed two unitary matrices P and U
such that:
B = PB′U, B′ =
[
H 0
C D
]
where H is lower triangular with positive diagonal coefficients. Initially, P and U are
identity matrices, H and C are empty and D = B. Let i be the index of the first row
of C and D.
Analysis for Array-OL 7
2. If D is null, the process stops.
3. If not, let j be the index of some non zero row of D. Let πij be the unitary matrix that
permutes rows i and j of B′. Since πij is its own inverse, one can write:
B = (Pπij)(πijB′)U,
and the new D has a non zero first row.
4. Let k be the index of a negative element in the first row of D. Let σk be the unit matrix
with the k-th diagonal element set to −1. Since σk is its own inverse, one can write:
B = P (B′σk)(σkU),
and element k in the first row of D is now positive.
5. If all elements in the first row of D are positive, let l be the index of the smallest element,
and let πil be the matrix that interchange columns i and l of B′. Again:
B = P (B′πil)(πilU)
and now the first element of the first row of D is smallest.
6. Let m > i be the index of some nonzero element in the first row of D. Set α = B′im÷B′ii.
By construction, α > 0. Let κim(α) be the identity matrix with −α added in position
(i, m). It is easy to see that the inverse of κim(α) is κim(−α). Hence:
B = P (B′κim(α))(κim(−α)U)
and element B′im has been replaced by B′im mod B′ii.
7. If the only non-zero element of the first row of D is the first element, then i can be
increased by 1.
These transformations must be applied until no further progress is possible (i.e. when in case
2). Matrix B′ is in the required form, and since all the elementary matrices π, σ and κ are
unitary, the resulting P and U are unitary. In fact, P is even a permutation matrix.
References
[1] Pierre Boulet. Array-ol revisited, multidemensional intensive signal processing specifica-
tion. Technical Report 6113, INRIA, January 2007.
[2] Eric Lenormand, 2007. Private communication.
[3] A. Schrijver. Theory of linear and integer programming. Wiley, NewYork, 1986.
[4] Rachid Seghir. Méthodes de dénombrement de points entiers de polyèdres et applications
à l’optimisation de programmes. PhD thesis, ULP, December 2006.


INRIA a CCSD electronic archive server

https://hal.inria.fr/inria-00146323/document

Elementary transformation analysis for Array-OL

Abstract

Similar works

Full text

Available Versions

Hal-Diderot

HAL Descartes

HAL-ENS-LYON

HAL Descartes

Archive Ouverte en Sciences de l'Information et de la Communication

INRIA a CCSD electronic archive server

Hal-Diderot