Data parallel language and compiler support for data intensive applications.
Processing and analyzing large volumes of data plays an increasingly important role in many domains of scientific research. High-level language and compiler support for developing applications that analyze and process such datasets has, however, been lacking so far.In this paper, we present a set of language extensions and a prototype compiler for supporting high-level object-oriented programming of data intensive reduction operations over multidimensional data. We have chosen a dialect of Java with data-parallel extensions for specifying a collection of objects, a parallel for loop, and reduction variables as our source high-level language. Our compiler analyzes parallel loops and optimizes the processing of datasets through the use of an existing run-time system, called active data repository (ADR). We show how loop fission followed by interprocedural static program slicing can be used by the compiler to extract required information for the run-time system. We present the design of a compiler/run-time interface which allows the compiler to effectively utilize the existing run-time system.A prototype compiler incorporating these techniques has been developed using the Titanium front-end from Berkeley. We have evaluated this compiler by comparing the performance of compiler generated code with hand customized ADR code for three templates, from the areas of digital microscopy and scientific simulations. Our experimental results show that the performance of compiler generated versions is, on the average 21% lower, and in all cases within a factor of two, of the performance of hand coded versions.
Introduction
Analysis and processing of very large multi-dimensional scientic datasets (i.e.
where data items are associated with points in a multidimensional attribute
space) is an important component of science and engineering. Examples of
these datasets include raw and processed sensor data from satellites [27], output
from hydrodynamics and chemical transport simulations [23], and archives
of medical images[1]. These datasets are also very large, for example, in medical
imaging, the size of a single digitized composite slide image at high power
from a light microscope is over 7GB (uncompressed), and a single large hospital
can process thousands of slides per day.
Applications that make use of multidimensional datasets are becoming increasingly
important and share several important characteristics. Both the input
and the output are often disk-resident. Applications may use only a subset
of all the data available in the datasets. Access to data items is described
by a range query, namely a multidimensional bounding box in the underlying
multidimensional attribute space of the dataset. Only the data items whose
associated coordinates fall within the multidimensional box are retrieved. The
processing structures of these applications also share common characteristics.
However, no high-level language support currently exists for developing applications
that process such datasets.
In this paper, we present our solution towards allowing high-level, yet ecient
programming of data intensive reduction operations on multidimensional
datasets. Our approach is to use a data parallel language to specify computations
that are to be applied to a portion of disk-resident datasets. Our solution
is based upon designing a prototype compiler using the titanium infrastructure
which incorporates loop ssion and slicing based techniques, and utilizing
an existing run-time system called Active Data Repository [8, 9, 10].
We have chosen a dialect of Java for expressing this class of computations. We
have chosen Java because the computations we target can be easily expressed
using the notion of objects and methods on objects, and a number of projects
are already underway for expressing parallel computations in Java and obtaining
good performance on scientic applications [4, 25, 36]. Our chosen dialect
of Java includes data-parallel extensions for specifying collection of objects,
a parallel for loop, and reduction variables. However, the approach and the
techniques developed are not intended to be language specic. Our overall
thesis is that a data-parallel framework will provide a convenient interface to
large multidimensional datasets resident on a persistent storage.
This research was supported by NSF Grant ACR-9982087, NSF Grant CCR-
9808522, and NSF CAREER award ACI-9733520.
Conceptually, our compiler design has two major new ideas. First, we have
shown how loop ssion followed by interprocedural program slicing can be
used for extracting important information from general object-oriented data-parallel
loops. This technique can be used by other compilers that use a run-time
system to optimize for locality or communication. Second, we have shown
how the compiler and the run-time system can use such information to eciently
execute data intensive reduction computations.
Our compiler extensively uses the existing run-time system ADR for optimizing
the resource usage during execution of data intensive applications. ADR
integrates storage, retrieval and processing of multidimensional datasets on a
parallel machine. While a number of applications have been developed using
ADR's low-level API and high performance has been demonstrated [9], developing
applications in this style requires detailed knowledge of the design of
ADR and is not suitable for application programmers. In comparison, our proposed
data-parallel extensions to Java enable programming of data intensive
applications at a much higher level. It is now the responsibility of the compiler
to utilize the services of ADR for memory management, data retrieval
and scheduling of processes.
Our prototype compiler has been implemented using the titanium infrastructure
from Berkeley [36]. We have performed experiments using three dierent
data intensive application templates, two of which are based upon the Virtual
Microscope application [16] and the third is based on water contamination
studies [23]. For each of these templates, we have compared the performance of
compiler generated versions with hand customized versions. Our experiments
show that the performance of compiler generated versions is, on average, 21%
lower and in all cases within a factor of two of the performance of hand coded
versions. We present an analysis of the factors behind the lower performance
of the current compiler and suggest optimizations that can be performed by
our compiler in the future.
The rest of the paper is organized as follows. In Section 2, we further describe
the charactestics of the class of data intensive applications we target.
Background information on the run-time system is provided in Section 3.
Our chosen language extensions are described in Section 4. We present our
compiler processing of the loops and slicing based analysis in Section 5. The
combined compiler and run-time processing for execution of loops is presented
in Section 6. Experimental results from our current prototype are presented
in Section 7. We compare our work with existing related research eorts in
Section 8 and conclude in Section 9.
2 Data Intensive Applications
In this section, we rst describe some of the scientic domains which involve
applications that process large datasets. Then, we describe some of the common
characteristics of the applications we target.
Data intensive applications from three scientic areas are being studied currently
as part of our project.
Analysis of Microscopy Data: The Virtual Microscope [16] is an application
to support the need to interactively view and process digitized data
arising from tissue specimens. The Virtual Microscope provides a realistic
digital emulation of a high power light microscope. The raw data for such a
system can be captured by digitally scanning collections of full microscope
slides under high power. At the basic level, it can emulate the usual behavior
of a physical microscope including continuously moving the stage and changing
magnication and focus. Used in this manner, the Virtual Microscope can
support completely digital dynamic telepathology.
Water contamination studies: Environmental scientists study the water
quality of bays and estuaries using long running hydrodynamics and chemical
transport simulations [23]. The chemical transport simulation models reactions
and transport of contaminants, using the
uid velocity data generated by the
hydrodynamics simulation. The chemical transport simulation is performed
on a dierent spatial grid than the hydrodynamics simulation, and also often
uses signicantly coarser time steps. To facilitate coupling between these two
simulation, there is a need for mapping the
uid velocity information from
the hydrodynamics grid, averaged over multiple ne-grain time steps, to the
chemical transport grid and computing smoothed
uid velocities for the points
in the chemical transport grid.
Satellite data processing: Earth scientists study the earth by processing
remotely-sensed data continuously acquired from satellite-based sensors, since
a signicant amount of Earth science research is devoted to developing correlations
between sensor radiometry and various properties of the surface of the
Earth [9]. A typical analysis processes satellite data for ten days to a year and
generates one or more composite images of the area under study. Generating a
composite image requires projection of the globe onto a two dimensional grid;
each pixel in the composite image is computed by selecting the \best" sensor
value that maps to the associated grid point.
Data intensive applications in these and related scientic areas share many
common characteristics. Access to data items is described by a range query,
namely a multidimensional bounding box in the underlying multidimensional
space of the dataset. Only the data items whose associated coordinates fall
within the multidimensional box are retrieved. The basic computation consists
of (1) mapping the coordinates of the retrieved input items to the corresponding
output items, and (2) aggregating, in some way, all the retrieved input
items mapped to the same output data items. The computation of a particular
output element is a reduction operation, i.e. the correctness of the output
usually does not depend on the order in which the input data items are aggregated

Another common characteristic of these applications is their extremely high
storage and computational requirements. For example, ten years of global
coverage satellite data at a resolution of four kilometers for our satellite data
processing application Titan consists of over 1.4TB of data [9]. For our Virtual
Microscope application, one focal plane of a single slide requires over 7GB
(uncompressed) at high power, and a hospital such as Johns Hopkins produces
hundreds of thousands of slides per year. Similarly, the computation for one
ten day composite Titan query for the entire world takes about 100 seconds
per processor on the Maryland sixteen node IBM SP2. The application scientists
typically demand real-time responses to such queries, therefore, e-cient
execution is extremely important.
3 Overview of the Run-time System
Our compiler eort targets an existing run-time infrastructure, called the Active
Data Repository (ADR) [9] that integrates storage, retrieval and processing
of multidimensional datasets on a parallel machine. We give a brief
overview of this run-time system in this section.
Processing of a data intensive data-parallel loop is carried out by ADR in two
phases: loop planning and loop execution. The objective of loop planning is to
determine a schedule to e-ciently process a range query based on the amount
of available resources in the parallel machine. A loop plan species how parts
of the nal output are computed. The loop execution service manages all the
resources in the system and carries out the loop plan generated by the loop
planning service. The primary feature of the loop execution service is its ability
to integrate data retrieval and processing for a wide variety of applications.
This is achieved by pushing processing operations into the storage manager
and allowing processing operations to access the buer used to hold data
arriving from the disk. As a result, the system avoids one or more levels of
copying that would be needed in a layered architecture where the storage
manager and the processing belong in dierent layers.
A dataset in ADR is partitioned into a set of (logical) disk blocks to achieve
high bandwidth data retrieval. The size of a logical disk block is a multiple
of the size of a physical disk block on the system and is chosen as a trade-o
between reducing disk seek time and minimizing unnecessary data transfers.
A disk block consists of one or more objects, and is the unit of I/O and
communication. The processing of a loop on a processor progresses through
the following three phases: (1) Initialization { output disk blocks (possibly
replicated on all processors) are allocated space in memory and initialized,
(2) Local Reduction { input disk blocks on the local disks of each processor are
retrieved and aggregated into the output disk blocks, and (3) Global Combine
{ if necessary, results computed in each processor in phase 2 are combined
across all processors to compute nal results for the output disk blocks.
ADR run-time support has been developed as a set of modular services implemented
in C++. ADR allows customization for application specic processing
(i.e., mapping and aggregation functions), while leveraging the commonalities
between the applications to provide support for common operations such as
memory management, data retrieval, and scheduling of processing across a
parallel machine. Customization in ADR is currently achieved through class
inheritance. That is, for each of the customizable services, ADR provides a
base class with virtual functions that are expected to be implemented by derived
classes. Adding an application-specic entry into a modular service requires
the denition of a class derived from an ADR base class for that service
and providing the appropriate implementations of the virtual functions. Current
examples of data intensive applications implemented with ADR include
Titan [9], for satellite data processing, the Virtual Microscope [16], for visualization
and analysis of microscopy data, and coupling of multiple simulations
for water contamination studies [23].
4 Java Extensions for Data Intensive Computing
In this section, we describe a dialect of Java that we have chosen for expressing
data intensive computations. Though we propose to use a dialect of Java as the
source language for the compiler, the techniques we will be developing will be
largely independent of Java and will also be applicable to suitable extensions
of other languages, such as C, C++, or Fortran 90.
4.1 Data-Parallel Constructs
We borrow two concepts from object-oriented parallel systems like Titanium [36],
HPC++ [5], and Concurrent Aggregates [11].
Interface Reducinterface f
*Any object of any class implementing
this interface is a reduction variable*
public class VMPixel f
char colors[3];
void Initialize() f
*Aggregation Function*
void Accum(VMPixel Apixel, int avgf) f
public class VMPixelOut extends VMPixel
implements Reducinterface;
public class VMScope f
static int
static int Ydimen = . ;
static
*Data Declarations*
static
static
static new VMPixel[VMSlide];
public static void main(String[] args) f
int
lowend)/subsamp];
foreach(p in Outputdomain) f
*Main Computational Loop*
foreach(p in querybox) f
Fig. 1. Example Code
Domains and Rectdomains are collections of objects of the same type. Rect-
domains have a stricter denition, in the sense that each object belonging
to such a collection has a coordinate associated with it that belongs to a
pre-specied rectilinear section of the domain.
The foreach loop, which iterates over objects in a domain or rectdomain, and
has the property that the order of iterations does not in
uence the result
of the associated computations. Further, the iterations can be performed in
parallel. We also extend the semantics of foreach to include the possibility
of updates to reduction variables, as we explain later.
We introduce a Java interface called Reducinterface. Any object of any class
implementing this interface acts as a reduction variable [18]. The semantics of a
reduction variable is analogous to that used in version 2.0 of High Performance
Fortran (HPF-2) [18] and in HPC++ [5]. A reduction variable has the property
that it can only be updated inside a foreach loop by a series of operations that
are associative and commutative. Furthermore, the intermediate value of the
reduction variable may not be used within the loop, except for self-updates.
4.2 Example Code

Figure

1 outlines an example code with our chosen extensions. This code
shows the essential computation in the virtual microscope application [16].
A large digital image is stored in disks. This image can be thought of as a
two dimensional array or collection of objects. Each element in this collection
denotes a pixel in the image. Each pixel is comprised of three characters, which
denote the color at that point in the image. The interactive user supplies two
important pieces of information. The rst is a bounding box within this two
dimensional box, which implies the area within the original image that the user
is interested in scanning. We assume that the bounding box is rectangular,
and can be specied by providing the x and y coordinates of two points. The
rst 4 arguments provided by the user are integers and together, they specify
the points lowend and hiend. The second information provided by the user
is the subsampling factor, an integer denoted by subsamp. The subsampling
factor tells the granularity at which the user is interested in viewing the image.
A subsampling factor of 1 means that all pixels of the original image must be
displayed. A subsampling factor of n means that n 2 pixels are averaged to
compute each output pixel.
The computation in this kernel is very simple. First, a querybox is created
using specied points lowend and hiend. Each pixel in the original image
which falls within the querybox is read and then used to increment the value
of the corresponding output pixel.
There are several advantages associated with specifying the analysis and processing
over multidimensional datasets in this fashion. The programs can specify
the computations assuming a single processor and
at memory. It also assumes
that the data is available in arrays of object references, and is not in
persistent storage. It is the responsibility of the compiler and run-time system
to locate individual elements of the arrays from disks. Also, it is the responsibility
of the compiler to invoke the run-time system for optimizing resource
usage.
4.3 Restrictions on the Loops
The primary goal of our compiler will be to analyze and optimize (by performing
both compile-time transformations and generating code for ADR run-time
system) foreach loops that satisfy certain properties. We assume standard semantics
of parallel for loops and reductions in languages like High Performance
Fortran (HPF) [18] and HPC++ [5]. Further, we require that no Java
threads be spawned within such loop nests, and no memory locations read or
written to inside the loop nests may be touched by another concurrent thread.
Our compiler will also assume that no Java exceptions are raised in the loop
nests and the iterations of the loop can be reordered without changing any
of the language semantics. One potential way of enabling this can be to use
bound checking optimizations [25].
5 Compiler Analysis
In this section, we rst describe how the compiler processes the given data-parallel
data intensive loop to a canonical form. We then describe how inter-procedural
program slicing can be used for extracting a number of functions
which are passed to the run-time system.
5.1 Initial Processing of the Loop
Consider any data-parallel loop in our dialect of Java, as presented in Section 4.
The memory locations modied in this loop are only the elements of collection
of objects, or temporary variables whose values are not used in other iterations
of the loop or any subsequent computations. The memory locations accessed
in this loop are either elements of collections or values which may be replicated
on all processors before the start of the execution of the loop.
For the purpose of our discussion, collections of objects whose elements are
modied in the loop are referred to as left hand side or lhs collections, and
the collections whose elements are only read in the loop are considered as right
hand side or rhs collections.
The functions used to access elements of collections of objects in the loop are
referred to as subscript functions.
Denition 1 Consider any two lhs collections or any two rhs collections.
These two collections are called congruent i
The subscript functions used to access these two collections in the loop are
identical.
The layout and partitioning of these two collections are identical. By identical
layout we mean that elements with the identical indices are put together
in the same disk block for both the collections. By identical partitioning we
mean that the disk blocks containing elements with identical indices from
these collections reside on the same processor.
Consider any loop. If multiple distinct subscript functions are used to access
rhs collections and lhs collections and these subscript functions are not
known at compile-time, tiling output and managing disk accesses while maintaining
high reuse and locality is going to be a very di-cult task for the
run-time system. In particular, the current implementation of ADR does not
support such cases. Therefore, we perform loop ssion to divide the original
loop into a set of loops, such that all lhs collections in any new loop are
congruent and all rhs collections are congruent. We now describe how such
loop ssion is performed.
Initially, we focus on lhs collections which are updated in dierent statements
of the same loop. We perform loop ssion, so that all lhs collections accessed
in any new loop are congruent. Since we are focusing on loops with no loop-carried
dependencies, performing loop ssion is straight-forward. An example
of such transformation is shown in Figure 2, part (a).
We now focus on such a new loop in which all lhs collections are congruent,
but not all rhs collections may be congruent. For any two rhs accesses in a
loop that are not congruent, there are three possibilities:
1. These two collections are used for calculating values of elements of dierent
lhs collections. In this case, loop ssion can be performed trivially.
2. These two collections Y and Z are used for calculating values of elements
of the same lhs collection. Such lhs collection X is, however, computed as
follows:
such that, op i  op j . In such a case, loop ssion can be performed, so that
the element X(f(i)) is updated using the operation op i with the values of
Y (g(i)) and Z(h(i)) in dierent loops. An example of such transformation is
shown in Figure 2, part (b).
3. These two collections Y and Z are used for calculating values of the elements
of the same lhs collection and unlike the case above, the operations used are
not identical. An example of such a case is
In this case, we need to introduce temporary collection of objects to copy
the collection Z. Then, the collection Y and the temporary collection can be
accessed using the same subscript function. An example of such transformation
is shown in Figure 2, part (c).
After such a series of loop ssion transformations, the original loop is replaced
by a series of loops. The property of each loop is that all lhs collections are
accessed with the same subscript function and all rhs collections are also
accessed with the same subscript function. However, the subscript function
for accessing the lhs collections may be dierent from the one used to access
rhs collections.
(a) foreach (p in box) f
foreach (p in box) f
foreach (p in box) f
(b) foreach (p in box) f
foreach (p in box) f
foreach (p in box) f
(c) foreach (p in box) f
foreach (p in box) f
foreach (p in box) f
Fig. 2. Examples of Loop Fission Transformations
O 1 [S L (r)] op
Om [S L (r)] op
Fig. 3. A Loop In Canonical Form
5.1.1 Discussion
Our strategy of performing loop ssion so that all lhs collections are accessed
with the same subscript function and all rhs collections are accessed with
the same subscript function is clearly not the best suited for all classes of
applications. Particularly, for stencil computations, it may result in several
accesses to each disk block. However, for the class of data intensive reductions
we have focused on, this strategy works extremely well and simplies the later
loop execution. In the future, we will incorporate some of the techniques from
parallel database join operations for loop execution, which will alleviate the
need for performing loop ssion in all cases.
5.1.2 Terminology
After loop ssion, we focus on one individual loop at a time. We introduce
some notation about this loop which is used for presenting our solution. The
terminology presented here is illustrated by the example loop in Figure 3.
The domain over which the iterator iterates is denoted by R. Let there be n
rhs collection of objects read in this loop, which are denoted by I
Similarly, let the lhs collections written in the loop be denoted by O
Further, we denote the subscript function used for accessing right hand side
collections by SR and the subscript function used for accessing left hand side
collections by S L .
int
Fig. 4. Slice for Subscript Function (left) and for Aggregation Function (right)
Given a point r in the range for the loop, elements S L (r) of the output collections
are updated using one or more of the values I 1
and other scalar values in the program. We denote by A i the function used for
creating the value which is used later for updating the element of the output
collection O i . The operator used for performing this update is op i .
5.2 Slicing Based Interprocedural Analysis
We are primarily concerned with extracting three sets of functions, the range
function R, the subscript functions SR and S L , and the aggregation functions
Similar information is often extracted by various data-parallel
Fortran compilers. One important dierence is that we are working with an
object-oriented language (Java), which is signicantly more di-cult to ana-
lyze. This is mainly because the object-oriented programming methodology
frequently leads to small procedures and frequent procedure calls. As a result,
analysis across multiple procedures may be required in order to extract range,
subscript and aggregation functions.
We use the technique of interprocedural program slicing for extracting these
three sets of functions. Initially, we give background information on program
slicing and give references to show that program slicing can be performed
across procedure boundaries, and in the presence of language features like
polymorphism, aliases, and exceptions.
5.2.1 Background: Program Slicing
The basic denition of a program slice is as follows. Given a slicing criterion
(s; x), where s is a program point in the program and x is a variable,
the program slice is a subset of statements in the programs such that these
statements, when executed on any input, will produce the same value of the
variable x at the program point s as the original program.
The basic idea behind any algorithm for computing program slices is as follows.
Starting from the statement p in the program, we trace any statements on
which p is data or control dependent and add them to the slice. The same is
repeated for any statement which has already been included in the slice, until
no more statements can be added in the slice.
ADR Pt outpoint(2);
ADR Pt lowend(2);
int
return outpoint ;
void Accumulate(ADR Box current block, ADR Box
current tile, ADR Box querybox) f
current block.intersect(querybox);
ADR Pt inputpt(2);
ADR Pt outputpt(2);
int
for
for
if (project(inputpt, outputpt, current tile)) f
Output[outputpt].Accum(VScope[inputpt],
Fig. 5. Compiler Generated Subscript and Aggregation Functions
Slicing has been very frequently used in software development environments,
for debugging, program analysis, merging two dierent versions of the code,
and software maintenance and testing. A number of techniques have been
presented for accurate program slicing across procedure boundaries [32]. Since
object-oriented languages have been frequently used for developing production
level software, signicant attention has been paid towards developing slicing
techniques in the presence of object-oriented features like object references,
polymorphism, and more recently, Java features like threads and exceptions.
Harrold et al. and Tonnela et al. have particularly focused on slicing in the
presence of polymorphism, object references, and exceptions [17, 28]. Slicing
in the presence of aliases and reference types has also been addressed [3].
5.2.2 Extracting Range Function
We need to determine the rhs and lhs collection of objects for this loop. We
also need to provide the range function R.
The rhs and lhs collection of objects can be computed easily by inspecting
the assignment statements inside the loop and in any functions called inside
the loop. Any collection which is modied in the loop is considered a lhs
collection, and any other collection touched in the loop is considered a rhs
collection.
For computing the domain, we inspect the foreach loop and look at the domain
over which the loop iterates. Then, we compute a slice of the program using
the entry of the loop as the program point and the domain as the variable.
5.2.3 Extracting Subscript Functions
The subscript functions SR and S L are particularly important for the run-time
system, as it determines the size of lhs collections written in the loop and the
rhs disk blocks from each collection that contributes to the lhs collections.
The function S L can be extracted using slicing as follows. Consider any statement
in the loop which modies any lhs collection. We focus on the variable
or expression used to access an element in the lhs collection. The slicing criterion
we choose is the value of this variable or expression at the beginning of
the statement where the lhs collection is modied.
The function SR can be extracted similarly. Consider any statement in the
loop which reads from any rhs collection. The slicing criterion we use is the
value of the expression used to access the collection at the beginning of such
a statement.
Typically, the value of the iterator will be included in such slices. Suppose the
iterator is p. After rst encountering p in the slice, we do not follow data-dependencies
for p any further. Instead, the functions returned by the slice
use such iterator as the input parameter.
For the virtual microscope template presented in Figure 1, the slice computed
for the subscript function S L is shown at the left hand side of Figure 4 and the
code generated by the compiler is shown on the left hand side of Figure 5. In
the original source code, the rhs collection is accessed with just the iterator
therefore, the subscript function SR is the identity function. The function
receives the coordinates of an element in the rhs collection as parameter
(iterpt) from the run-time system and returns the coordinates of the corresponding
lhs element. Titanium multidimensional points are supported by
ADR as a class named ADR Pt. Also, in practice, the command line parameters
passed to the program are extracted and stored in a data-structure, so that
the run-time system does not need to explicitly read args array.
5.2.4 Extracting Aggregation Functions
For extracting the aggregation function A i , we look at the statement in the
loop where the lhs collection O i is modied. The slicing criterion we choose is
the value of the element from the collection which is modied in this statement,
at the beginning of this statement.
The virtual microscope template presented in Figure 1 has only one aggregation
function. The slice for this aggregation function is shown in Figure 4 and
the actual code generated by the compiler is shown in Figure 5. The function
Accum accessed in this code is obviously part of the slice, but is not shown here.
The generated function iterates over the elements of a disk block and applies
aggregation functions on each element, if that element intersects with the range
of the loop and the current tile. The function is presented as a parameter of
current block (the disk block being processed), the current tile (the portion
of lhs collection which is currently allocated in memory), and querybox
which is the iteration range for the loop. Titanium rectangular domains are
supported by the run-time as ADR Box. Further details of this aggregation
function are explained after presenting the combined compiler/run-time loop
processing.
6 Combined Compiler and Run-time Processing
In this section we explain how the compiler and run-time system can work
jointly towards performing data intensive computations.
6.1 Initial Processing of the Input
The system stores information about how each of the rhs collections of objects
I i is stored across disks. Note that after we apply loop ssion, all rhs
collections accessed in the same loop have identical layout and partitioning.
The compiler generates appropriate ADR functions to analyze the meta-data
about collections I i , the range function R, and the subscript function SR , and
compute the list of disk blocks of I i that are accessed in the loop. The domain
of each rhs collection accessed in the loop is SR  R. Note that if a disk block
is included in this list, it is not necessary that all elements in this disk block
are accessed during the loop. However, for the initial planning phase, we focus
on the list of disk blocks.
We assume a model of parallel data intensive computation in which a set of
disks is associated with each node of the parallel machine. This is consistent
with systems like IBM SP and cluster of workstations. Let the set
denote the list of processors in the system. Then, the information computed
by the run-time system after analyzing the range function, the input subscript
function and the meta-data about each of the collections of objects I i is the
. For a given input collection I i and a processor is the set of disk
blocks b that contain data for collection I i , is resident on a disk connected to
processor intersects with SR  R.
Further, for each disk block b ijk belonging to the set B ij , we compute the
information D(b ijk ), which denotes the subset of the domain SR  R which is
resident on the disk block b. Clearly the union of the domains covered by all
selected disk blocks will cover the entire area of interest, or in formal terms,
6.2 Work Partitioning
One of the issues in processing any loop in parallel is work or iteration parti-
tioning, i.e., deciding which iterations are executed on which processor.
The work distribution policy we use is that each iteration is performed on the
owner of the element read in that iteration. This policy is opposite to the owner
computes policy [19] which has been commonly used in distributed memory
compilers, in which the owner of the lhs element works on the iteration. The
rationale behind the approach is that the system will not have to communicate
blocks of the rhs collections. Instead, only replicated elements of the lhs
collections need to be communicated to complete the computation. Note that
the assumptions on the nature of loops we have placed requires replacing an
initial loop by a sequence of canonical loops, which may also increase net
communication between processors. However, we do not consider it to be a
problem for the set of applications we target.
6.3 Allocating Output Buers and Strip Mining
The distribution of rhs collections is decided before performing the processing,
and we have decided to partition the iterations accordingly. We now need to
allocate buers to accumulate the local contributions to the nal lhs objects.
We use run-time analysis to determine the elements of output collections which
are updated by the iterations performed on a given processor. This run-time
analysis is similar to the one performed by run-time libraries for executing
irregular applications on distributed memory machines. Any element which is
updated by more than one processor is initially replicated on all processors
by which it is updated. Several dierent strategies for allocation of buers
have been studies in the context of the run-time system [9]. Selecting among
these dierent strategies for the compiler generated code is a topic for future
research.
The memory requirements of the replicated output space are typically higher
than the available memory on each processor. Therefore, we need to divide
the replicated output buer into chunks that can be allocated on the main
memory of each processor. This is the same issue as strip mining or tiling used
for improving cache performance. We have so far used only a very simple strip
mining strategy. We query the run-time system to determine the available
memory that can be allocated on a given processor. Then, we divide the lhs
space into blocks of that size. Formally, we divide the lhs domain S L  R
into a set of smaller domains (called strips) fS 1 g. Since each of the
lhs collection of objects in the loop is accessed through the same subscript
function, the same strip mining is done for each of them.
In performing the computations on each processor, we will iterate over the set
of strips, allocate that strip for each of the n output collections, and compute
local contributions to each strip, before allocating the next strip. To facilitate
this, we compute the set of rhs disk blocks that will contribute to each strip
of the lhs.
6.4 Mapping Input to the Output
We use subscript functions SR and S L for computing the set of rhs disk blocks
that will contribute to each strip of the lhs as indicated above. To do this
we apply the function S L  (S 1
R ) to each D(b ijk ) to obtain the corresponding
domain in the lhs region. These output domains that each disk block can
contribute towards are denoted as OD(b ijk ). If D(b ijk ) is a rectangular domain
and if the subscript functions are monotonic, OD(b ijk ) will be a rectangular
domain and can easily be computed by applying the subscript function to the
two extreme corners. If this is not the case, the subscript function needs to be
applied on each element of D(b ijk ) and the resulting OD(b ijk ) will just be a
domain and not a rectangular domain. Formally, we compute the sets L jl , for
each processor j and each output strip l, such that
6.5 Actual Execution
The computation of sets L il marks the end of the loop planning phase of the
run-time system. Using this information, now the actual computations are
performed on each processor. The structure of the computation is shown in

Figure

6. In practice, the computation associated with each rhs disk block and
retrieval of disk blocks is overlapped, by using asynchronous I/O operations.
We now explain the aggregation function generated by the compiler for the
virtual microscope template presented in Figure 1, shown in Figure 5 on the
right hand side. The accumulation function output by the compiler captures
the Foreach element part of the loop execution model shown in Figure 6. The
For each output strip S l :
Execute on each Processor
Allocate and initialize strip S l for O
Foreach
Read blocks b ijk disks
Foreach element e of D(b ijk )
If the output pt. intersects with S l
Evaluate functions A
Global reduction to nalize the values for S l
Fig. 6. Loop Execution on Each Processor
run-time system computes the sets L jl as explained previously and invokes the
aggregation function in a loop that iterates over each disk block in such a set.
The current compiler generated code computes the rectangular domain D(b ijk )
in each invocation of the aggregation function, by doing an intersection of the
current block and query block. The resulting rectangular domain is denoted
by box.
The aggregation function iterates over the elements of box. The conditional if
project() achieves two goals. First, it applies subscript functions to determine
the lhs element outputpt corresponding to the rhs element inputpt. Second,
it checks if outputpt belongs to current tile. The actual aggregation is
applied only if outputpt belongs to current tile. This test represents a
signicant source of ine-ciency in the compiler generated code. If the tile
or strip being currently processed represents a rectangular rhs region and
the subscript functions are monotonic, then the intersection of the box and
the tile can be performed before the inner loop. This is in fact done in the
hand customization of ADR for virtual microscope [16]. Performing such an
optimization automatically is a topic for future research and beyond the scope
of our current work.
7 Current Implementation and Experimental Results
In this section, we describe some of the features of the current compiler and
then present experimental results comparing the performance of compiler generated
customization for three codes with hand customized versions.
Our compiler has been implemented using the publicly available Titanium
infrastructure from Berkeley [36]. Our current compiler and run-time system
only implement a subset of the analysis and processing techniques described
in this paper. Two key limitations are as follows. We can only process codes
in which the rhs subscript function is the identity function. It also requires
that the domain over which the loop iterates is a rectangular domain and
all subscript functions are monotonic. Titanium language is an explicitly parallel
dialect of Java for numerical computing. We could use the Titanium
front-end without performing any modications. Titanium language includes
Point, RectDomain, and foreach loop which we required for our purposes.
The concept of reducinterface is not part of Titanium language, but no
modications to the parser were required for this purpose. Titanium also includes
a large number of additional directives which we did not require, and
has signicantly dierent semantics for foreach loops.
We have used three templates for our experiments.
VMScope1 is identical to the code presented in Figure 1. It models a virtual
microscope, which provides a realistic digital emulation of a microscope,
by allowing the users to specify a rectangular window and a subsampling
factor. The version VMScope1 averages the colors of neighboring pixels to
create a pixel in the output image.
VMScope2 is similar to VMScope1, except for one important dierence. Instead
of taking the average of the pixels, it only picks every subsamp th
element along each dimension to create the nal output image. Thus, only
memory accesses and copying is involved in this template, no computations
are performed.
Bess models computations associated with water contamination studies over
bays and estuaries. The computation performed in this application determines
the transport of contaminants, and accesses large
uid velocity data-sets
generated by a previous simulation.
Virtual Microscope (averaging)103050701 2 4 8
# of processors
Execution
Time
Compiled
Original
Fig. 7. Comparison of Compiler and Hand Generated Versions for VMScope1
These three templates represent data intensive applications from two important
domains, digital microscopy and scientic simulations. The computations
and data accesses associated with these computations are relatively simple and
can be handled by our current prototype compiler. Moreover, we had access to
hand coded ADR customization for each of these three templates. This allowed
us to compare the performance of compiler generated versions against hand
coded versions whose performance had been reported in previously published
work [16, 23].
Our experiments were performed using the ADR run-time system ported on a
cluster of dual processor 400 MHz Intel Pentium nodes connected by gigabit
ethernet. Each node has 256MB main memory and GB of internal disk.
Experiments were performed using 1, 2, 4 and 8 nodes of the cluster. ADR run-time
system's current design assumes a shared nothing architecture and does
not exploit multiple CPUs on the same node. Therefore, only one processor
on each node was used for our experiments.
The results comparing the performance of compiler generated and hand customized
VMScope1 are shown in Figure 7. A microscope image of 19; 760
15; 360 pixels was used. Since each pixel in this application takes 3 bytes, a
total of 910 MB are required for storage of such an image. A query with a
bounding box of size 10; 00010; 000 with a subsampling factor of 8 was used.
The time taken by the compiler generated version ranged from 73 seconds on
1 processor to 13 seconds on 8 processors. The speedup on 2, 4, and 8 processors
was 1.86, 3.32, and 5.46, respectively. The time taken by the hand coded
version ranged from 68 seconds on 1 processor to 8.3 seconds on 8 processors.
The speedup on 2, 4, and 8 processors was 2.03, 4.09, and 8.2, respectively.
Since the code is largely I/O and memory bound, slightly better than linear
speedup is possible. The performance of compiler generated code was lower
by 7%, 10%, 25%, and 38% on 1, 2, 4, and 8 processors respectively.
From this data, we see that the performance of compiler generated code is
very close to the hand coded one on the 1 processor case, but is substantially
lower on the 8 processor case. We carefully compared the compiler generated
and hand coded versions to understand these performance factors. The two
codes use dierent tiling strategies of the lhs collections. In the hand coded
version, an irregular strategy is used which ensures that each input disk block
maps entirely into a single tile. In the compiler version, a simple regular tiling
is used, in which each input disk block can map to multiple tiles. As shown
in

Figure

5, the compiler generated code performs an additional check in each
iteration, to determine if the lhs element intersects with the tile currently
being processed. In comparison, the tiling strategy used for the hand coded
version ensures that this check always returns true, and therefore does not
need to be inserted in the code. But, because of the irregular tiling strategy,
an irregular mapping is required between the bounding box associated with
each disk block and the actual coordinates on the allocated output tile. This
mapping needs to be carried out after each rhs disk block is read from the
memory. The time required for performing such mapping is proportional to the
number of rhs disk blocks processed by each processor for each tile. Since the
output dataset is actually quite small in our experiments, the number of rhs
disk blocks processed by each processor per tile decreases as we go to larger
congurations. As a result, the time required for this extra processing reduces.
In comparison, the percentage overhead associated with extra checks in each
iteration performed by the compiler generated version remains unchanged.
This dierence explains why the compiler generated code is slower than the
hand coded one, and why the dierence in performance increases as we go to
larger number of processors.
Virtual Microscope (subsampling)103050
# of processors
Execution
Time
Compiled
Original
Fig. 8. Comparison of Compiler and Hand Generated Versions for VMScope2
The results comparing the performance of compiler generated and hand coded
VMScope2 are shown in Figure 8. This code was executed on the same dataset
and query as used for VMScope1. The time taken by the compiler generated
version ranged from 44 seconds on 1 processor to 9 seconds on 1 processor.
The hand coded version took 47 seconds on 1 processor and nearly 5 seconds
on 8 processors. The speedup of the compiler generated version was 2.03, 3.31,
and 4.88 on 2, 4, and 8 processors respectively. The speedup of the hand coded
version was 2.38, 4.98, 10.0 on 2, 4, and 8 processors respectively.
A number of important observations can be made. First, though the same
query is executed for VMScope1 and VMScope2 templates, the execution times
are lower for VMScope2. This is clearly because no computations are performed
in VMScope2. However, a factor of less than two dierence in execution times
shows that both the codes are memory and I/O bound and even in VMScope1,
the computation time does not dominate. The speedups for hand coded version
of VMScope2 are higher. This again is clearly because this code is I/O and
memory bound.
The performance of the compiler generated version was better by 6% on 1
processor, and was lower by 10% on 2 processors, 29% on 4 processors, and
48% on 8 processors. This dierence in performance is again because of the
dierence in tiling strategies used, as explained previously. Since this template
does not perform any computations, the dierence in the conditionals and extra
processing for each disk block has more signicant eect on the overall
performance. In the 1 processor case, the additional processing required for
each disk block becomes so high that the compiler generated version is slightly
faster. Note that the hand coded version was developed for optimized execution
on parallel systems and therefore is not highly tuned for sequential case.
For the 8 processor case, the extra cost of conditional in each iteration becomes
dominant for the compiler generated version. Therefore, the compiler
generated version is almost a factor of 2 slower than the hand coded one.
Bays and Estuaries Simulation System2060100140
# of processors
Execution
Time
Compiled
Original
Fig. 9. Comparison of Compiler and Hand Generated Versions for Bess
The results comparing performance of compiler generated and hand coded
version for Bess are shown in Figure 9. The dataset comprises of a grid with
2113 elements and 46,080 time-steps. Each time-step has 4 4-byte
oating
point numbers per grid element, denoting simulated hydrodynamics parameters
previously computed. Therefore, the memory requirements of the dataset
are 1.6 GB. The Bess template we experimented with performed weighted
averaging of each of the 4 values for each column, over a specied number of
time-steps. The number of time-steps used for our experiments was 30,000.
The execution times for the compiler generated version ranged from 131 seconds
on 1 processor to 11 seconds on 8 processors. The speedup on 2, 4, and
8 processors was 1.98, 5.53, and 11.75 respectively. The execution times for
the hand coded version ranged from 97 seconds on 1 processor to 9 seconds
on 8 processors. The speedup on 2, 4, and 8 processors was 1.8, 5.4, and 10.7
respectively. The compiler generated version was slower by a factor of 25%,
19%, 24%, and 19% on 1, 2, 4, and 8 processors respectively.
We now discuss the factors behind the dierence in performance of compiler
generated and hand coded Bess versions. As with both the VMScope versions,
the compiler generated Bess performs checks for intersecting with the tile for
each pixel. The output for this application is very small, and as a result, the
hand coded version explicitly assumes a single output tile. The compiler generated
version cannot make this assumption and still inserts the conditionals.
The amount of computation associated with each iteration is much higher for
this application. Therefore, the percentage overhead of the extra test is not
as large as the VMScope templates. The second important dierence between
the compiler generated and hand coded versions is how averaging is done. In
the compiler generated code, each value to be added is rst divided by the
total number of values which are being added. In comparison, the hand coded
version performs the summation of all values rst, and then performs a single
division. The percentage overhead of this is independent of the number of processors
used. We believe that the second factor is the dominant reason for the
dierence in performance of two versions. This also explains why the percentage
dierence in performance remains unchanged as the number of processors
is increased. The performance of compiler generated code can be improved by
performing the standard strength reduction optimization. However, the compiler
needs to perform this optimization interprocedurally, which is a topic for
future work.
As an average over these three templates and 1, 2, 4, and 8 processor congu-
rations, the compiler generated versions are 21% slower than hand coded ones.
Considering the high programming eort involved in managing and optimizing
disk accesses and computations on a parallel machine, we believe that a 21%
slow-down from automatically generated code will be more than acceptable
to the application developers. It should also be noted that the previous work
in the area of out-of-core and data intensive compilation has focused only
on evaluating the eectiveness of optimizations, and not on any comparisons
against hand coded versions.
Our analysis of performance dierences between compiler generated and hand
coded versions has pointed us to a number of directions for future research.
First, we need to consider more sophisticated tiling strategies to avoid large
performance penalties associated with performing extra tests during loop ex-
ecution. Second, we need to consider more advanced optimizations like inter-procedural
code motion and interprocedural strength reduction to improve the
performance of compiler generated code.
8 Related Work
Our work on providing high-level support for data intensive computing can
be considered as developing an out-of-core Java compiler. Compiler optimizations
for improving I/O accesses have been considered by several projects. The
PASSION project at Northwestern University has considered several dierent
optimizations for improving locality in out-of-core applications [6, 20]. Some
of these optimizations have also been implemented as part of the Fortran D
compilation system's support for out-of-core applications [29]. Mowry et al.
have shown how a compiler can generate prefetching hints for improving the
performance of a virtual memory system [26]. These projects have concentrated
on relatively simple stencil computations written in Fortran. Besides
the use of an object-oriented language, our work is signicantly dierent in
the class of applications we focus on. Our techniques for executions of loops
are particularly targeted towards reduction operations, whereas previous work
has concentrated on stencil computations. Our slicing based information extraction
for the runtime system allows us to handle applications which require
complex data distributions across processors and disks and for which only
limited information about access patterns may be known at compile-time.
Many researchers have developed aggressive optimization techniques for Java,
targeted at parallel and scientic computations. javar and javab are compilation
systems targeting parallel computing using Java [4]. Data-parallel
extensions to Java have been considered by at least two other projects: Titanium
[36] and HP Java [7]. Loop transformations and techniques for removing
redundant array bounds checking have been developed [12, 25]. Our eort is
also unique in considering persistent storage, complex distributions of data on
processors and disks, and the use of a sophisticated runtime system for optimizing
resources. Other object-oriented data-parallel compilation projects
have also not considered data residing on persistent storage [5, 11, 30].
Program slicing has been actively used for many software engineering applications
like program based testing, regression testing, debugging and software
maintenance over the last two decades [34]. In the area of parallel compi-
lation, slicing has been used for communication optimizations by Pugh and
Rosser [31] and for transforming multiple levels of indirection by Das and
Saltz [15]. We are not aware of any previous work on using program slicing
for extracting information for the runtime system.
Several research projects have focused on parallelizing irregular applications,
such as computational
uid dynamics codes on irregular meshes and sparse
matrix computations. This research has demonstrated that by developing run-time
libraries and through compiler analysis that can place these runtime calls,
such irregular codes can be compiled for e-cient execution [2, 21, 24, 35]. Our
project is related to these eorts in the sense that our compiler also heavily
uses a runtime system. However, our project is also signicantly dierent. The
language we need to handle can have aliases and object references, the applications
involve disks accesses and persistent storage, and the runtime system
we need to interface to works very dierently.
Several runtime support libraries and le systems have been developed to support
e-cient I/O in a parallel environment [13, 14, 22, 33]. They also usually
provide a collective I/O interface, in which all processing nodes cooperate to
make a single large I/O request. Our work is dierent in two important ways.
First, we are supporting a much higher level of programming by involving a
compiler. Second, our target runtime system, ADR, also diers from these
systems in several ways. The computation is an integral part of the ADR
framework. With the collective I/O interfaces provided by many parallel I/O
systems, data processing usually cannot begin until the entire collective I/O
operation completes. Also, data placement algorithms optimized for range
queries are also integrated as part of the ADR framework.
9 Conclusions
In this paper we have addressed the problem of expressing data intensive computations
in a high-level language and then compiling such codes to e-ciently
manage data retrieval and processing on a parallel machine. We have developed
data-parallel extensions to Java for expressing this important class of
applications. Using our extensions, the programmers can specify the computations
assuming that there is a single processor and
at memory.
Conceptually, our compiler design has two major new ideas. First, we have
shown how loop ssion followed by interprocedural program slicing can be
used for extracting important information from general object-oriented data-parallel
loops. This technique can be used by other compilers that use a run-time
system to optimize for locality or communication. Second, we have shown
how the compiler and run-time system can use such information to e-ciently
execute data intensive reduction computations. This technique for processing
such loops is independent of the source language.
These techniques have been implemented in a prototype compiler built using
the Titanium front-end. We have used three templates, from the areas of digital
microscopy and scientic simulations, for evaluating the performance of
this compiler. We have compared the performance of compiler generated code
with the performance of codes developed by customizing the run-time system
ADR manually. Our experiments have shown that the performance of compiler
generated codes is, on the average, 21% slower than the hand coded ones, and
in all cases within a factor of 2. We believe that these results establish that
our approach can be very eective. Considering the high programming eort
involved in managing and optimizing disk accesses and computation on a parallel
machine, we believe that a 21% slow-down from automatically generated
code will be more than acceptable to the application developers. It should also
be noted that the previous work in the area of out-of-core and data intensive
compilation has focused only on evaluating the eect of optimizations, and not
on any comparisons against hand coded versions. Further, we believe that by
considering more sophisticated tiling strategies and other optimizations like
interprocedural code motion and strength reduction, the performance of the
compiler generated codes can be further improved.

Acknowledgments

We are grateful to Chialin Chang, Anurag Acharya, Tahsin Kurc, Alan Sussman
and other members of the ADR team for developing the run-time system,
developing hand customized versions of applications, helping us with the ex-
periments, and for many fruitful discussions we had with them during the
course of this work.



--R

Angelo De- marzo
Interprocedural compilation of irregular applications for distributed memory machines.

A prototype Java restructing compiler.
Distributed pC
A model and compilation strategy for out-of-core data parallel programs

A customizable parallel database for multi-dimensional data
Infrastructure for building parallel database systems for multi-dimensional data
Alan Suss- man
Concurrent aggregates (CA).

The Vesta parallel
Input/Output characteristics of Scalable Parallel Applications.
Paul Havlak
The Virtual Microscope.

High Performance Fortran Forum.
Compiling Fortran D for MIMD distributed-memory machines
Improving the performance of out-of-core computations
Compiling global name-space parallel loops for distributed execution

Coupling multiple simulations via a high performance customizable database system.
Exploiting spatial regularity in irregular iterative applications.

Automatic compiler-inserted i/o prefetching for out-of-core applications
NASA Goddard Distributed Active Archive Center (DAAC).

Compiler support for out-of- core arrays on parallel machines
object-oriented languages
Iteration space slicing and its application to communication optimization.
Speeding up slicing.

A survey of program slicing techniques.

--TR
Dynamic slicing in the presence of unconstrained pointers
Compiling Fortran D for MIMD distributed-memory machines
object-oriented languages
Speeding up slicing
A model and compilation strategy for out-of-core data parallel programs
Interprocedural compilation of irregular applications for distributed memory machines
Input/output characteristics of scalable parallel applications
Index array flattening through program transformation
The Vesta parallel file system
Automatic compiler-inserted I/O prefetching for out-of-core applications
Flow insensitive C++ pointers and polymorphism analysis and its application to slicing
Iteration space slicing and its application to communication optimization
Reuse-driven interprocedural slicing
Concurrent aggregates (CA)
Passion
Distributed Memory Compiler Design For Sparse Problems
Compiling Global Name-Space Parallel Loops for Distributed Execution
Titan
Improving the Performance of Out-of-Core Computations
Infrastructure for Building Parallel Database Systems for Multi-Dimensional Data
Exploiting spatial regularity in irregular iterative applications
Compiler support for out-of-core arrays on parallel machines
data intensive applications;data parallel language;run-time support;compiler techniques