An Efficient Solution to the Cache Thrashing Problem Caused by True Data Sharing.
AbstractWhen parallel programs are executed on multiprocessors with private caches, a set of data may be repeatedly used and modified by different threads. Such data sharing can often result in cache thrashing, which degrades memory performance. This paper presents and evaluates a loop restructuring method to reduce or even eliminate cache thrashing caused by true data sharing in nested parallel loops. This method uses a compiler analysis which applies linear algebra and the theory of numbers to the subscript expressions of array references. Due to this method's simplicity, it can be efficiently implemented in any parallel compiler. Experimental results show quite significant performance improvements over existing static and dynamic scheduling methods.
Introduction
Parallel processing systems with memory hierarchies have become quite common today. Commonly,
most multiprocessor systems have a local cache in each processor to bridge the speed gap between the
processor and the main memory. Some systems use multi-level caches [5, 14]. Very often, a copy-back
snoopy cache protocol is employed to maintain cache coherence in these multiprocessor systems. Certain
supercomputers also use a local memory which can be viewed as a program-controlled cache. When
programs with nested parallel loops are executed on such parallel processing systems, it is important to
assign parallel loop iterations to the processors in such a way that unnecessary data movement between
different caches is minimized. For convenience, in this paper, we call each iteration of a parallel loop a
thread. The following loop nest is an example.
DO I=1,100
DO K=1,100
ENDDO
ENDDO
In this example, loop J is a parallel loop because, with the value of I fixed, statement S has no
loop-carried dependences, while J varies. Loop J is executed 100 times in the loop nest, creating 10,000
iterations, or 10,000 threads, in total. Each thread addresses 100 elements of array A. Many array
elements are repeatedly accessed by these threads as shown in Table 1 and Figure 1, where T i;j denotes
the thread corresponding to loop index values I = i and denotes the set of threads
created by the index value I = i. As shown in Figure 2, there exist lists of threads: (T 1;1 ), (T 1;2 , T 2;1 ),
that each thread modifies and reuses most of the
array elements accessed by the neighboring threads in the same list. If the threads in the same list are
assigned to different processors, the data of array A will unnecessarily move back and forth between
different caches in the system, causing a cache thrashing problem due to true data sharing [12].
The nested loop construct shown in the above example is quite common in parallel code used for
scientific computation. J. Fang and M. Lu studied a large number of programs including the LINPACK
benchmarks, the PERFECT Club benchmarks, and programs for mechanical CAE, computational
chemistry, image and signal processing, and petroleum applications [11]. They reported that almost
all of the most time-consuming loop nests contain at least three loop levels, out of which 60% contain
at least one parallel loop. Even after using loop interchange to move parallel loops outwards when
it was legal, they still found 94% of the parallel loops enclosed by sequential loops. Such loop nests
include the cases in which a parallel loop appears in the outermost loop level in a subroutine, but the
subroutine is called by a call-statement which is contained by a sequential loop. Most of these loop
nests are not perfectly nested, i.e. there exist statements right before or after an inner loop. Fang and
Lu proposed a thread alignment algorithm to solve the cache thrashing problem which may be created
by such multi-nested loops. Their algorithm, however, assigns threads to processors either by solving
linear equations at run time or by storing the precomputed numerical solutions to the equations in the
processor's memory. Since storing all the numerical solutions requires quite a large memory space, and
since the exact number of threads often cannot be determined statically due to unknown loop bounds,
they favor the on-line computation approach. In this paper, we present a method to reduce the run-time
overhead in Fang and Lu's algorithm by using a thorough compiler analysis of array references
to derive a closed-form formula that solves the thread assignment equations. The thread assignment
then becomes highly efficient at run time. Previously, we presented preliminary algorithms [22, 23] to
deal with a simple case in which data-dependent array references use the same linear function in the
subscripts. No experimental data were given. In this paper, we extend the work by covering multiple
linear functions and by clarifying the underlying theory. We report experimental results using a Silicon
Graphics (SGI) multiprocessor.

Table

1: The elements of A accessed by each thread.
With our method, the compiler analyzes the data dependences between the threads and uses that
information to restructure the nested loop, perfectly nested or otherwise, in order to reduce or even
eliminate true data sharing, which causes cache thrashing. Our method can be efficiently implemented
in any parallel compiler, and our experimental results show quite significant improvement over existing
static and dynamic scheduling methods.
In what follows, we first address related work. We then introduce basic concepts and assumptions.
After that, we present solutions to the cache thrashing problem due to true data sharing, and lastly we
show the experimental results conducted on an SGI multiprocessor system.
Related Work
Extensive research regarding efficient memory hierarchies has been reported in the literature.
Abu-Sufah, Kuck and Lawrie use loop blocking to improve paging performance by improving the
locality of references [2]. Wolfe proposes iteration space tiling as a way to improve data reuse in a cache
or a local memory [35]. Gallivan, Jalby and Gannon define a reference window for a dependence as the
variables referenced by both the source and the sink of the dependence [15, 16]. After executing the
source of the dependence, they save the associated reference window in the cache until the sink has also
T1,3

Figure

1: The elements of A accessed by different threads.
been executed, which may increase the number of cache hits. Carr, Callahan and Kennedy [7, 8] discuss
options for compiler control of a uniprocessor's memory hierarchy. Wolf and Lam develop an algorithm
that estimates all temporal and spatial reuse of a given loop permutation [34]. These optimizations all
attempt to maximize the reuse of cached data on a single processor. They also have a secondary effect
of improving multiprocessor performance by reducing the bandwidth requirement of each processor,
thereby reducing contention in the memory system. In contrast, our work considers a multiprocessor
environment where each processor has its own local cache or its own local memory, and where different
processors may share data.
The work by Peir and Cytron [28], Shang and Fortes [30], and by D'Hollander [9] share the common
goal of partitioning an index set into independent execution subsets such that the corresponding loop
iterations can execute on different processors without interprocessor communication. Their methods
apply to a specific type of loop nest called a uniform recurrence or a uniform dependence algorithm, in
which the loops are perfectly nested, the loop bounds are constant, the loop-carried dependences have
constant distances, and the array subscripts are of the form is a loop index and c an
integer constant. Hudak and Abraham [1, 18] develop a static partitioning approach called adaptive data
partitioning (ADP) to reduce interprocessor communication for iterative data-parallel loops. They also
assume perfectly nested loops. The loop body is restricted to update a single data point A(i;
a two-dimensional global matrix A. The subscript expressions of right-hand side array references are
restricted to be the sum of a parallel loop index and a small constant, while the subscript expressions of
left-hand array references are restricted to contain the parallel loop indices only. Tomko and Abraham
T1,3
T100,3

Figure

2: Lists of threads accessing similar array elements.
[32] develop iteration partitioning techniques for data-parallel application programs. They assume that
there is only one pair of data access functions and that each loop index variable can appear in only one
dimension of each array subscript expression. Agarwal, Kranz, and Natarajan [3] propose a framework
for automatically partitioning parallel loops to minimize cache coherence traffic on shared-memory
multiprocessors. They restrict their discussion to perfectly nested doall loops, and assume rectangular
iteration spaces. Unlike these previous works, our work considers nested loops which are not necessarily
perfectly nested. Loop bounds can be any variables, and array subscript expressions are much more
general. Many researchers have studied the cache false sharing problem [10, 17, 19, 33] in which cache
thrashing occurs when different processors share the same cache line of multiple words, although the
processors do not share the same word. Many algorithms have been proposed to reduce false sharing by
better memory allocation, better thread scheduling, or by program transformations. Our work considers
cache thrashing which is due to the true sharing of data words.
Our work is most closely related to the research done by Fang and Lu [11, 12, 26, 13]. In their work,
the iteration space is partitioned into a set of equivalence classes, and each processor uses a formula to
determine which iterations belong to the same equivalence class at execution time. Each processor then
executes the corresponding iterations so as to reduce or eliminate cache thrashing. These iterations are
the solution vectors of a linear integer system. In Fang and Lu's work, these vectors may either be
computed at run time or may be precomputed and later retrieved at run time when loop bounds are
known before execution. Both approaches require additional execution time when a processor fetches
the next iteration. Unlike Fang and Lu's approaches, we solve the thrashing problem at compile time to
reduce run-time overhead, while we achieve the same effect of reducing cache thrashing. Our new method
restructures the loops at compile time and is based on a thorough analysis of the relationship between
the array element accesses and the loop indices in the nested loop. We have performed experiments on
a commercial multiprocessor, namely a Silicon Graphics Challenge Cluster, thereby obtaining real data
regarding cache thrashing and its reduction. In contrast, previous data were mainly from simulations
[11, 12, 26].
3 Basic Concepts and Assumptions
Data dependences between statements are defined in [6, 24, 4, 25]. If a statement S 1 uses the result
of another statement S 2 , then S 1 is flow-dependent on S 2 . If S 1 can safely store its result only after
fetches the old data stored in that location, then S 1 is anti-dependent on S 2 . If S 1 overwrites the
result of S 2 , then S 1 is output-dependent on S 2 . A dependence within an iteration of a loop is called
a loop-independent dependence. A dependence across the iterations of a loop is called a loop-carried
dependence.
There can be no cache thrashing due to true data sharing if the outermost loop is parallel, because
no data dependences will cross the threads. Therefore, in this paper, we consider only loop nests whose
outermost loops are sequential. To simplify our discussion, we make the following assumptions about
the program pattern: 1) All functions representing array subscript expressions are linear. 2) The loop
construct considered here consists of a sequential loop which embraces one or several single-level parallel
loops. If there exist multilevel parallel loops, only one level is parallelized, as on most commercial
shared-memory multiprocessor systems. Hence, as shown below, a loop nest in our model has three
levels: a parallel loop, its immediately enclosed sequential loop, and its immediately enclosing sequential
loop:
ENDDO
ENDDO

Figure

3: The loop nest model.
are linear mappings from the iteration space N 1 \Theta N 2 \Theta N 3 to the domain space M 1 \Theta M 2
of A:
and they can be expressed as:
array k). The loops in the above example are not necessarily perfectly
nested. Our restructuring techniques, to be presented later, assume arbitrary loop bounds, although
we are showing lower bounds of 1 here for simplicity of notation. Multiple array variables and multiple
linear subscript functions may exist in the nested loop.
Since we are considering cache thrashing due to true data sharing, i.e. due to data dependences
between threads, we can also write the loop nest in Figure 3 as:
ENDDO
ENDDO
where A fl m) is an array name appearing in the loop body, ~ h
m) are linear
mappings from iteration space N 1 \Theta N 2 \Theta N 3 to domain space M 1 fl
\Theta M 2 fl
of A fl ,
A
m) are potentially dependent reference pairs, and m is the number of such
pairs. Fang and Lu [11] reported that arrays involved in nested loops are usually two-dimensional or
three-dimensional with a small-sized third dimension. The latter can be treated as a small number of
two-dimensional arrays. Nested loops with the parallel loop at the innermost level are degenerate cases
of the loop nest in Figure 3. Therefore, our loop nest model seems quite general. Our method can also
be applied to a loop nest which contains several separate parallel loops at the middle level. Each of
these parallel loops may be restructured according to its own reference patterns, such that the threads in
different instances of the same parallel loop are aligned. We currently do not align the threads created
by different parallel inner loops. For programs with more complicated loop nests, pattern-matching
techniques can be used to identify a loop subnest that matches the nest shown in Figure 3. Other
outer- or inner- loops that are not a part of the subnest can be ignored, as long as their loop indices do
not appear in the array subscripts.
The compiler analysis is based on a simple multiprocessor model in which the cache memory has the
following characteristics: 1) it is local to a processor; 2) it uses a copy-back snoopy coherence strategy;
and its line size is one word. The transformed code, however, will execute correctly on machines
which have a multiword cache line and multilevel caches. Furthermore, as the experimental results will
show, the performance of the transformed code is quite good on such realistic machines. Our analysis
can also be extended to incorporate more machine parameters such as the cache line size.
Solutions
In this section, we analyze the relationship between linear functions in the array subscripts. Based
on this analysis, we restructure a given loop nest to reduce or eliminate the cache thrashing due to
true data sharing. We consider nested loops which are not necessarily perfectly nested and which may
have variable loop bounds. For clarity of presentation, in Section 4.1 we first discuss how to deal with
dependent reference pairs such that the same subscript function is used in both references in each pair
(Different pairs may use different subscript functions). Later in Section 4.2, we will discuss how to deal
with more general cases by using simple affine transforms to fit to this model.
4.1 The Basic Model
In this subsection, we assume that for each pair of dependent references, the same subscript function is
used in both references. Under this assumption, if we extract the subscript function ~ h fl (I; J; K) from
each pair of dependent references, then a model for a nested loop which has m pairs of dependent
references can be illustrated by the following code segment.
ENDDO
ENDDO
Without loss of generality, suppose that all m linear subscript functions above are different. We assume
that each function ~ h m, is of rank 2 and is in the form of ~ h fl (i; j;
where
Take the following example.
ENDDO
ENDDO

Figure

4: A nested loop with multiple linear subscript functions.
In this example, no data dependences exist within loop J . However, two data dependences exist in
the whole loop nest, one between the references to A, and the other between those to B. We have two
linear functions to consider, one for each dependence:
The iteration subspace N 1 \Theta N 2 is called the reduced iteration space because it omits the K loop. In
order to find the iterations in the reduced iteration space which may access common memory locations
within the corresponding threads, we define a set of elements of array A fl which are accessed within
thread using subscript function ~ h fl .
Definition 1: Given iteration (i in the reduced iteration space, the elements A fl (f fl (i
are accessed within thread
. They are
denoted by A i 0 ;j 0
g.
Definition 2: If we suppose that T i;j and T i 0 ;j 0 are two threads corresponding to (I; J)=(i;
) in the reduced iteration space of the given loop nest such that A i;j
m); we say T has a dependence because of ~ h fl , denoted by T i;j
Definition 3: If there exists fl, 1 - fl - m, such that T i;j
Since both f fl and g fl are linear in terms of i; j and k, the following lemma is obvious.
Lemma 1: In the program pattern described above, if there exist fl, 1 - fl - m, and two iterations
in the iteration space N 1 \Theta N 2 \Theta N 3 , (i, j,
, such that
then for any constant n 0 , we have a series of iterations in the space, (i;
that satisfy the following equations:
The following lemma and two theorems establish the relationship between the loop indexes corresponding
to two inter-dependent threads. We will use this index relationship to stagger the loop
iteration space such that inter-dependent threads can be assigned to the same processors.
Lemma 2: Let T i;j , T i 0 ;j 0 be two threads,
. T i;j exist k, k 0
such that
a
a
Theorem 1: Let b fl;1 a exist k, k 0
b fl;1 a fl;2 \Gamma a fl;1 b fl;2
a fl;2 c fl;1 \Gamma a fl;1 c fl;2
b fl;1 a fl;2 \Gamma a fl;1 b fl;2
The proofs of Lemma 2 and Theorem 1 are obvious from Definition 3. We now consider the case
of b fl;1 a fl;2 \Gamma a fl;1 b assuming that the loop bounds, N 2 and N 3 , are large enough to satisfy the
c fl;1
These assumptions are almost always true in practice [31]. When they are not true, the parallel loops
will be too small to be important. With these assumptions, we have the following theorem.
Theorem 2 [21]: Let b fl;1 a fl;2 \Gamma a fl;1 b The fact that the J loop at the middle level is a
loop guarantees that T i;j
(1) a fl;1 (i 0
(2) a fl;2 (i 0
In order to find the threads which have data dependence relations with thread T i;j , we make the
following definition.
Definition 4: Given iteration (i; j) in the reduced iteration space, we let S i;j denote the following
set of iterations in the space:
r
r
where L fl;1 (L fl;1 6= 0) and L fl;2 (1 - fl - m) are defined as:
and
with GCD fl equal to G:C:D:(b fl;1 c fl;2 \Gammac fl;1 b fl;2 ; a fl;2 c fl;1 \Gammaa fl;1 c fl;2 ; b fl;1 a fl;2 \Gammaa fl;1 b fl;2 ) or equal to \GammaG:C:D:
(b fl;1 c fl;2 - c fl;1 b fl;2 , a fl;2 c fl;1 -a fl;1 c fl;2 , b fl;1 a fl;2 -a fl;1 b fl;2 ) to guarantee L fl;1 ? 0;
(2)
and L
with GCD fl equal to G:C:D:(a fl;1 ; b fl;1 ) or equal to \GammaG:C:D:(a fl;1 ; b fl;1 ) to guarantee L fl;1 ? 0;
and L
with GCD fl equal to G:C:D:(a fl;2 ; b fl;2 ) or equal to \GammaG:C:D:(a fl;2 ; b fl;2 ) to guarantee L fl;1 ? 0;
called the staggering parameter corresponding to linear function ~ h fl . If there exist no
data dependences between the given pair of references, we define the staggering parameter (L fl;1 ; L fl;2 )
as (0; 0).
The staggering parameters for the example in Figure 4 can be calculated to be: (L 1;1 ; L 1;2
and (L 2;1 ; L 2;2 according to Definition 4(1).
The following theorem, derived from Theorems 1 and 2 and Definition 4, states that we can use the
staggering parameters to uniquely partition the threads into independent sets.
Theorem 3: S i;j as defined above satisfies:
(1) if (i;
The theorem above indicates that S i;j includes all the iterations whose corresponding threads have
a data dependence relation with T i;j . We call S i;j an equivalence class of the reduced iteration space.
In order to eliminate true data sharing, threads in the same equivalence class should be assigned to
the same processor. We want to restructure the reduced iteration space such that threads in the same
equivalence class will appear in the same column. Each staggering parameter (L computed for a
dependent reference pair tells us that if we stagger the (i row in the reduced iteration space
by columns to the right if L 2 ! 0, or to the left if relative to the i-th, then the threads
involved in the dependence pair will be aligned in the same column. Different staggering parameters
may require staggering the iteration space in different ways. However, if these staggering parameters
are in proportion, then staggering by the unified staggering parameter defined below will satisfy all the
requirements simultaneously.
Definition 5: Given staggering parameters (L
, and then we call (g; L1;2
the unified staggering parameter.
Lemma 3 [21]: If the condition L k;1
m) in Definition 5 is true, then (a) the
iterations (i; belong to two different equivalent classes; and (b) the iterations
belong to two different equivalence classes.
Theorem 4 [21]: If the condition L k;1
m) in Definition 5 is true, then the
reduced iteration space must be staggered according to the unified staggering parameter (g; L1;2
L1;1 g) in
order to reduce or eliminate data sharing among the threads, i.e. the (i g)-th row in the reduced
iteration space must be staggered by j L1;2
L1;1 gj columns to the right if L 1;2 ! 0, or to the left if L 1;2 ? 0,
relative to the i-th row.
If a given loop nest satisfies the condition L k;1
m) in Definition 5, then, according
to Theorem 4 above, the reduced iteration space can be transformed into a staggered and reduced
iteration space (SRIS) by leaving the first g rows unchanged,
staggering each of the remaining rows using the unified staggering parameter. There will be no data
dependences between different columns in the SRIS.
However, if the staggering parameters are not in proportion, i.e, if there exist (j; k) such that
, then we can no longer obtain a unique unified staggering parameter.
Moreover, staggering alone is no longer sufficient for eliminating data dependences between the different
columns in the restructured iteration space. This is because some threads in the same equivalence class
are still in different columns. We perform a procedure called compacting which stacks these columns
onto each other. We will discuss staggering first.
Definition Given staggering parameters (L
(L 1;1 , L 2;1 , ., L m;1 ), suppose there exists (j; k) such that 1 -
. According to
the theory of numbers [27], there exist integers a 1 , a 2 , ., am that satisfy
a
a fl L fl;2 . We call (g; g 0
unified staggering parameter.
Note that since the m-tuple (a 1 , a 2 , ., am ) is not necessarily unique, the (g; g 0
may not be
unique either. With Definition 6, a unified staggering parameter (g; g 0
) of the example in Figure 4 is
found to be
After staggering by using any unified staggering parameter (g; g 0
), the resulting SRIS has four
possible shapes, as shown in Figure 5(b \Gamma e). Figure 5(a) shows details of one of these shapes. Figure
6(a) and 6(b) show the reduced iteration space for the example in Figure 4 before and after staggering
with (3,-3) as the unified staggering parameter.
Next, we compute the compacting parameter d using Algorithm 1 and 2, to be presented shortly.
We then partition the SRIS into n chunks, where
d
, which is the total number of
columns in the SRIS devided by the compacting parameter d (Figure 5(d \Gamma e)). These d-wide chunks
are stacked onto each other to form a compacted iteration space of width d, as shown in Figure 7. As we
will explain later, the threads in different columns after compacting the SRIS with d are independent.
Moreover, the product of d and g equals the number of equivalence classes. The SRIS shown in Figure
6(b) for our example is transformed by being compacted with into the form shown in Figure 8.
The following algorithm computes the compacting parameter d.
Algorithm 1:
Input: A set of staggering parameters (L
Output: The compacting parameter d.
Step 1: For each 2-element subset, fL i;1 ; L j;1 g, of fL 1;1 ; L 2;1 ; :::; L m;1 g, compute
of all such d 2 hL i;1 ; L j;1 i.
Step 2: For each j-element subset, fL pick any element, say
a) g'<0, g>1
d
d
d) g'<0, g>1 e) g'>0, g>1

Figure

5: SRIS and outlines.
Using the Euclidean Algorithm, compute integers b
Apply Algorithm 2 below to find nonzero integers r 2 ; :::; r j such that
r
a) Original reduced iteration space
(1,1) .

Figure

The reduced iteration space before and after rearrangement.
Let
r
Step 3: For j from 3 to m, compute
Step 4:
As will be established later, d is unique regardless of the choice of L i 1 ;1 in Step 2.
To calculate the compacting parameter d, non-zero integers r need to be found in Algorithm 1
from the integer coefficients b computed by the Euclidean Algorithm. Algorithm 2 is therefore
invoked to derive a group of non-zero integer coefficients from a group of any integer coefficients of a
linear expression.
Algorithm 2:
Input: Non-zero positive integers such that
Output: non-zero integers such that
1: If there are an even number of zero coefficients a
(0 - 2k - p) among
d)
d)
d)
d)
d)
d)
d)
d)
d)
d)
d)
a

Figure

7: Compacted SRIS.
Step 2: If there are an odd number of zero coefficients a i 1
(0
Obviously, the non-zero integers computed by Algorithm 2 satisfy
For the example in Figure 4, since there are only two linear functions in the loop nest, only Step 1
and Step 4 of Algorithm 1 are used to calculate the compacting parameter d, that is,
Next, we need to establish two important facts. First, after compacting with d, the threads in different
columns are independent. Second, the compacting parameter d computed by Algorithm 1 is the

Figure

8: The reduced iteration space after compacting.
largest number of independent columns possible as the result of compacting the SRIS with a constant
value. The first fact is established by Theorem 5, Theorem 6, and Corollary 1. To do so, we introduce
the following definition.
Definition 7: Given an iteration (i; j) in the reduced iteration space, staggering parameters
and the unified staggering parameter (g; g 0
are integers that satisfy
a
a set of iterations S 0
i;j is constructed as follows:
(1) For any integer r, iteration (i
) in the space belongs to S 0
a
(2) If there exist integers r not all zero, integer r, and iterations (i 0
the space, such that
r
r
and
a
The following three lemmas and Theorem 5 show that S 0
i;j is the same as the equivalence class S i;j .
From the process of constructing S 0
i;j , we immediately have the following lemma.
Lemma 4: Given iterations (i;
) in the reduced iteration space, and a unified staggering
parameter (g; g 0
), if there exists integer r such that
Lemma 5 [21]: Given iterations (i 0
) in the reduced iteration space, if there
exist integers r all zero, such that
r fl L fl;2
r
Lemma 6 [21]: Given iterations (i;
) in the reduced iteration space, if (i 0
Theorem 5 [21]: Given staggering parameters (L
for any iteration (i; j) in the reduced iteration space.
Next, we establish that S 0
i;j is the result of staggering with (g; g 0
followed by compacting with d.
This is stated by Corollary 1 below.
Lemma 7 [21]: Given staggering parameters (L
that are integers that satisfy
(1)
r
(2) if there exist integers r 0
satisfying
1 .
For any integers r 00
satisfying
r 00
there exists an integer k - 1 such that r 00
Theorem 6 [21]: If d is the compacting parameter determined by Algorithm 1, and d 0
r fl L fl;2 ,
are integers, not all zero, which satisfy
r
then there exists an integer k such that d 0
Corollary 1: The set S 0
i;j in Definition 7 satisfies
where k and r are integers, (g; g 0
) is the unified staggering parameter in Definition 6, and d is the
compacting parameter computed by Algorithm 1.
From the above result, the threads in different columns after compacting the SRIS with d are inde-
pendent. Next, we establish with Theorem 7 that any two columns which are d columns apart, where d
is computed by Algorithm 1, should be dependent and that, therefore, d is the largest possible number
of independent columns as the result of compacting the SRIS with a constant number.
Theorem 7: Given (i; j), we have
i;j .
Proof: According to how d is computed in Algorithm 1, there exist integers r such that
By the definition of S 0
To further simplify the process of the staggering and the compacting of the reduced iteration space,
the following theorem can be used to replace multiple staggering parameters, which are in proportion,
with a single staggering parameter.
Theorem 8 [21]: Given staggering parameters (L
m) are the staggering parameters satisfying
there exists an integer r satisfying
We now estimate the time needed by the compiler to compute the staggering parameters, a unified
staggering parameter, and the compacting parameter. Suppose there are m reference pairs. The
complexity of determining all the staggering parameters is O(m). A unified staggering parameter
of these staggering parameters can be determined in O(m) with the Euclidean Algorithm. Let
be the number of groups of staggering parameters such that all parameters in the same group are
in proportion. m 0
is very small in practice. According to Theorem 8, we only need to consider one
representative from each group. The complexity of Algorithm 1 and 2 for computing the compacting
parameter is C 2
Lastly, we show the result of restructuring the original loop nest based on staggering and compacting.
Note that if all staggering parameters are in proportion, then compacting is unnecessary for data
dependence elimination. However, to improve load balance, we compact the SRIS by a compacting
factor d, which equals the number of the available processors. The restructured code is parameterized
by the loop bounds and by the number of available processors, which can be obtained by a system call
at runtime. There is no need to recompile for a different number of available processors.
If the given loop nest is perfectly nested, then the resulting code, after staggering and compacting,
is shown in Code Segment 1 listed below.
Code Segment 1 (The result after restructuring the perfectly nested loop with multiple linear
ENDDO
ENDDO
ENDDO
If the given loop nest is not perfectly nested, then the resulting code has two variants, one for
and the other for g 0
=0 . We show the code for g 0
listed below (the code for
=0 is similar [21]). In this code segment, LB 1 and LB 2 are the lower bounds of I and J , UB 1 and
UB 2 are the upper bounds of these two loops, (g; g 0
and d are the unified staggering parameter and
the compacting parameter, respectively, which have been determined above. PSI and PSJ are local
variables that each processor uses to determine the first iteration J 0 of loop J to be executed on it.
Variables J 0 and OFFSET are also local to each processor. PSI and PSJ for each processor are modified
every g iterations of the loop I, according to the staggering and the compacting parameters. The
values of (PSI; PSJ) are initialized for the d different processors to (LB
respectively. We define a function mod* such that x mod*
Code Segment 2 (The restructured code for the case of g 0 6= 0):
endif
endif
ENDDO
ENDDO
ENDDO
4.2 An Extended Model
The theory we developed in the previous subsection can be extended to more general cases in which the
subscript functions in the same pair of references are not necessary the same. Suppose the following
two linear functions
~
and
~
belong to the same pair of references. In order to determine which iterations in the reduced iteration
space are dependent due to this reference pair, we consider an affine transformation
such that the linear function ~
h 2 can be expressed as
~
a 2;1
a
a 2;2
which we denote by ~
). In order to use the previous results from Section 4.1, we let ~
be identical to ~
which implies
a 1;1
a 2;1
a 1;2
a 2;2
c 1;1
c 2;1
c 2;2
and
a 2;1
a 2;2
We can now apply the algorithms in Section 4.1 to ~
2 and ~
which yield a staggering parameter,
say For a given iteration (i 0
). The iteration (i; must
have a dependence with (i 00
before the affine transformation if and only if the iteration (i 0
a dependence with (i 00
after the transformation. We denote the distance between (i;
as (L 0
2 ), which can be calculated as:
or
such that L 0
not be constant, meaning that the iterations cannot be
aligned with a constant staggering parameter. In common practice, since loop J is DOALL in our loop
nest model, the two linear functions ~
will have the same cofficients for loop index variables I
and J , which implies that ff 1. In this paper, we will consider the case of ff
We now have
or
We define (L 0
2 ), which are two constants given staggering parameter in this case.
If Equations (1) and (2) have a unique solution for we have a unique staggering
parameter (L 0
On the other hand, if there exist multiple solutions for then the following theorem shows
that under certain conditions, (L 0
determined by different should be in proportion.
Theorem 9: Assume ff 1. If the staggering parameter (L of the subscript function
~
after the affine transformation is a solution for Equations (1) and (2), then (L 0
is equal to (fi or to (L proportion
with solution (fi Equations (1) and (2).
proportion with prove that
in proportion to (L supposing that Every solution to Equations (1)
and (2) can be written as 0
solution of the homogeneous system associated with Equations (1) and (2), that
is,
a
a
So, if a 1;2 b 1;1 \Gamma b 1;2 a 1;1 6= 0, we have
a 2;2 b
a 1;2 b 1;1 \Gamma b 1;2 a 1;1
a 2;2 b
c 1;2 a
a 1;2 b 1;1 \Gamma b 1;2 a 1;1
Suppose
then we have c 1;2 a
according to Definition 4, we have
c 1;2 a
For the case of a 1;2 b as in Theorem 2, we have: (1) a 1;1
(2) a 1;2 Therefore, according to Definition 4,
we also have (- proportion with are in
proportion with
2 ) is in proportion with
If the condition in Theorem 9 is met, we choose (L 0
as the staggering
parameter for the reference pair ~

Table

shows examples of staggering parameters for different subscript functions appearing in the
dependent reference pair, where the loop index variables are listed in the order from the outermost loop
level to the innermost. If we simultaneously consider two reference pairs: A(I; J) with
and B(I; J) with then the thread T i;j will share the same array element A(i;
thread T i+3;j+1 and the same array element B(i; j) with thread T i+1;j+3 . Using Theorem 9, the staggering
parameters (L 0
these two pairs are (3,1) and (1,3) respectively. A unified staggering
parameter and compacting parameter can be calculated as (g; g 0
8.

Table

2: Examples of different functions in the same dependent reference pair.
Loop nest Dependent reference pair ( ~
5 Experimental Results
The thread alignment techniques described in this paper have been implemented as backend optimizations
in KD-PARPRO [20], a knowledge-based parallelizing tool which can perform intra- and inter-procedural
data dependence analysis and a large number of parallelizing transformations, including loop
interchange, loop distribution, loop skewing, and strip mining for FORTRAN programs.
To evaluate the effect of the thread alignment techniques on the performance of multiprocessors with
memory hierarchies, we experimented with three programs from the LINPACK benchmarks on a SGI
Challenge cluster which can be configured to contain up to twenty MIPS 4400 processors. First, the
programs were parallelized and optimized using KD-PARPRO. To reduce or eliminate cache thrashing
due to true data sharing, our tool recognized the nested loops which may cause the thrashing. It
applied the techniques described in the previous section to analyze and restructure the loop nests. The
parallelized programs were then compiled using SGI's f77 compiler with the optimization option -O2.
The sequential versions of the programs were compiled on the same machine using the same optimization
option for f77. The output binary codes were then executed on various configurations with a different
number of processors during dedicated time. Each MIPS 4400 processor has a 16K-byte primary data
cache and a 4M-byte secondary cache. The cache block size is 32 bytes for the primary data cache and
128 bytes for the secondary cache. A fast and wide split transaction bus POWERpath-2 is used as its
coherent interconnect. Cache coherence is maintained with a snoopy write-invalidate strategy.
We compared the results obtained by using our algorithm to align the threads with those obtained
by using four different loop scheduling strategies provided by SGI system software, namely, simple,
interleave, dynamic, and gss. The simple method divides the iterations by the number of processors
and then assigns each chunk of consecutive iterations to one processor. The interleave scheduling
method divides the iterations into chunks of the size specified by the CHUNK option, and execution of
those chunks is statically interleaved among the processes. With dynamic scheduling, the iterations are
also divided into CHUNK-sized chunks. As each process finishes a chunk, however, it enters a critical
section to grab the next available chunk. With gss scheduling [29] , the chunk size is varied, depending
on the number of iterations remaining. None of these SGI-provided methods consider task alignment.
The speedup of a parallel execution on shared memory machines just like the SGI cluster can
be affected by many factors, including: program parallelism, data locality, scheduling overhead, and
load balance. Usually gss, dynamic and interleave schedulings with a small chunk size are supposed
to show better load balance than simple scheduling. On the other hand, they tend to incur more
scheduling overhead than simple. Furthermore, simple captures more data locality in most cases than
other schedulings do.
The programs we selected from LINPACK are SGEFA, SPODI, and SSIFA. SGEFA factors a double
precision matrix by Gaussian elimination. The main loop structure in this program consists of three
imperfectly nested loops. The innermost loop is inside subroutine SAXPY, which multiplies a vector
by a constant and then adds the result to another vector. In order to show the array access pattern
inside the loop body, we inlined the SAXPY in the code section given below. However, we kept the
subroutine call when we applied our techniques to the program.
ENDDO
ENDDO
The SRIS, after staggering the iterations in the reduced iteration space (K; J), is shown in Figure
9. For this program, we only consider the linear function (I; J). The staggering parameter is (1,0)
according to Definition 4(4). The number of processors is used to determine the compacting factor.

Figure

9: SRIS for SGEFA.
SPODI computes the determinant and inverse of a certain double precision symmetric positive-definite
matrix. There are two main loop nests in this program, as shown below. We restructured both
loop nests. The same as in SGEFA, the innermost loop is contained in the subroutine SAXPY.
/* the first loop nest */
ENDDO
ENDDO
/* the second loop nest */
ENDDO
ENDDO
The SRISs, after staggering the iterations in the reduced iteration spaces (K; J) and (J; K) for
these two loop nests, respectively, are shown in Figure 10(a)(b). The linear functions we considered are:
(I; J) and (I; K), respectively. Their staggering parameters are both (1,0), according to Definition 4(4).
(a) Loop nest 1
(b) Loop nest 2

Figure

10: SRISs for SPODI.
SSIFA factors a double precision symmetric matrix by elimination. The main loop nest in this
program is shown below. We view the backward GOTO loop as the outermost sequential loop, within
which the value of kstep may change between 1 and 2 in different iterations, based on the input matrix.
Depending on the value of kstep, one of the two parallel loop nests inside the outermost sequential loop
will be executed for each iteration of the outermost loop. The index step of the outermost loop equals
i.e. \Gamma1 or \Gamma2. The array access patterns for these two kstep values are slightly different. The
innermost loop is again inside the SAXPY subroutine.
10: CONTINUE
IF (K .EQ.
ENDDO
ENDDO
ENDDO
20: CONTINUE
The SRISs, after staggering the iterations in the reduced iteration space (K; JJ) for two different
ksteps, are shown in Figure 11. For clarity, we use c kstep to denote the value of kstep in the current
K iteration, and we use p kstep for its value in the previous K iteration. All threads will be aligned
well if we properly align the threads created in the current K iteration with those in the previous K
iteration. We need to consider four possible cases of dependences: one is between two references to
another is between two references to A(I; the third is from A(I; K \Gamma JJ)
to A(I; and the last is from A(I; K \Gamma For all of these cases, the
staggering parameter is (L
(a)

Figure

11: SRISs for SSIFA.
The problem sizes we used in our experiments are n=100 and 1000. The performance of the parallel
codes transformed by our techniques, compared with the performance achieved by the scheduling
methods provided by SGI, are shown in Figures 12-15. Both SGEFA and SSIFA may require pivoting
for non-positive-definite symmetric matrices, but not for positive-definite symmetric matrices. We show
data for SGEFA with pivoting, and we show data both with and without pivoting for SSIFA. Pivoting
may potentially destroy the task alignment.1352 4 8
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(a)
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(b) n=1000

Figure

12: Gaussian elimination (SGEFA).12345
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(a) n=10026101418
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(b) n=1000

Figure

13: Determinant and inverse of a symmetric positive matrix (SPODI).
As the figures show, our method always outperforms all of the SGI's scheduling methods, with the
exception of program SGEFA. For this program, our method's performance is almost the same as that
of simple, although our method outperforms simple by 14% on 16 processors with should
be attributed to the reduction of cache thrashing due to true data sharing, a problem that tends to
be more severe when more processors are running. The simple scheduling method tends to get better
performance than the dynamic, gss, and interleave methods, because it results in better locality and
less cache thrashing in most cases, and it also incurs less scheduling overhead. But when the programs
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(a) n=100246810
Number of processors
Interleave (chunk
Interleave (chunk
Dynamic
Gss
Simple
Our method
(b) n=1000

Figure

14: Factorization of a symmetric matrix (SSIFA).
do not exhibit a good load balance, like SPODI, SSIFA, and other programs in LINPACK, which deal
with symmetric matrices, simple's performance results degrade substantially. Our method outperforms
simple quite significantly in most cases, especially for SPODI (Figure 13), as well as for SSIFA without
pivoting (Figure 15), where our method beats simple by as much as 105%. We are not able to get
improvement over simple for program SGEFA when pivoting is much more likely to
destroy the locality we try to keep. For the rest of the programs, we attribute our performance gain
over simple both to the reduction of cache thrashing due to true data sharing and to a better load
balance, although, for SSIFA with pivoting, we believe our method benefits more from load balancing.
We note that the SGI system software cannot pick the right scheduling method automatically to fit the
particular program. On the other hand, our method seems more capable of delivering good performance
for different loop shapes. As the thrashing problem becomes more serious on parallel systems with more
processors and greater communication overhead, our method will likely be even more effective.
6 Conclusions
This paper presents a method in which the reduced iteration space is rearranged according to the staggering
and the compacting parameters. The nested loop (either perfectly nested or imperfectly nested)
is restructured to reduce or even eliminate cache thrashing due to true data sharing. This method can
be efficiently implemented in any parallel compiler. Although the analysis per se is based on a simple
machine model, the resulting code executes correctly on more complex models. Our experimental
results show that the transformed code can perform quite well on a real machine. How to extend the
techniques proposed in this paper to incorporate additional machine parameters is interesting future
work.



--R


On the performance enhancement of paging systems through program analysis and transformations.
Automatic partitioning of parallel loops and data arrays for distributed shared-memory multiprocessors
Automatic loop interchange.
Multilevel cache hierarchies: organizations
Dependence analysis for supercomputing.
Improving register allocation for subscripted variables.
Compiling scientific code for complex memory hierarchies.
Partitioning and labeling of loops by unimodular transformations.
Eliminating False Sharing.
A solution of cache ping-pong problem in RISC based parallel processing systems
Cache or local memory thrashing and compiler strategy in parallel processing systems.
An iteration partition approach for cache or local memory thrashing on parallel processing.
Performance optimizations
On the problem of optimizing data transfers for complex memory systems.
Strategies for cache and local memory management by global program transformation.
Effects of program parallelization and stripmining transformation on cache performance in a multiprocessor.
Compiler techniques for data partitioning of sequentially iterated parallel loops.
Reducing false sharing on shared memory multiprocessors through compile-time data transformations
The design and the implementation of a knowledge-based parallelizing tool
An efficient solution to the cache thrashing problem (Extended Version).
Loop restructuring techniques for the thrashing problem.
Loop staggering
The structure of computers and computations.
Dependence graphs and compiler op- timizations
A solution of the cache ping-pong problem in multiprocessor systems
An introduction to the theory of numbers.
Minimum distance: a method for partitioning recurrences for multiproces- sors
Guided self-scheduling: a practical scheduling scheme for parallel supercomputers
Time optimal linear schedules for algorithms with uniform dependencies.
An empirical study of Fortran programs for parallelizing compilers.
Iteration partitioning for resolving stride conflicts on cache-coherent multiprocessors
False sharing and spatial locality in multiprocessor caches.
A data locality optimizing algorithm.
More iteration space tiling.
--TR
multiprocessors;parallelizing compilers;parallel threads;loop transformations;cache thrashing;true data sharing