Synthesizing Efficient Out-of-Core Programs for Block Recursive Algorithms Using Block-Cyclic Data Distributions.
AbstractIn this paper, we present a framework for synthesizing I/O efficient out-of-core programs for block recursive algorithms, such as the fast Fourier transform (FFT) and block matrix transposition algorithms. Our framework uses an algebraic representation which is based on tensor products and other matrix operations. The programs are optimized for the striped Vitter and Shriver's two-level memory model in which data can be distributed using various $cyclic(B)$ distributions in contrast to the normally used physical track distribution $cyclic(B_d)$, where $B_d$ is the physical disk block size. We first introduce tensor bases to capture the semantics of block-cyclic data distributions of out-of-core data and also data access patterns to out-of-core data. We then present program generation techniques for tensor products and matrix transposition. We accurately represent the number of parallel I/O operations required for the synthesized programs for tensor products and matrix transposition as a function of tensor bases and data distributions. We introduce an algorithm to determine the data distribution which optimizes the performance of the synthesized programs. Further, we formalize the procedure of synthesizing efficient out-of-core programs for tensor product formulas with various block-cyclic distributions as a dynamic programming problem. We demonstrate the effectiveness of our approach through several examples. We show that the choice of an appropriate data distribution can reduce the number of passes to access out-of-core data by as large as eight times for a tensor product and the dynamic programming approach can largely reduce the number of passes to access out-of-core data for the overall tensor product formulas.
Introduction
Due to the rapid increase in the performance of processors and communication networks in the last two
decades, the cost of memory access has become the main bottleneck in achieving high-performance for
many applications. Modern computers, including parallel computers, use a sophisticated memory hierarchy
consisting of, for example, caches, main memory, and disk arrays, to narrow the gap between the processor
and memory system performance. However, the efficient use of this deep memory hierarchy is becoming more
and more challenging. For out-of-core applications, such as computational fluid dynamics and seismic data
processing, which involve a large volume of data, the task of efficiently using the I/O subsystem becomes
y Supported by NSF Grant NSF-IRI-91-00681, Rome Labs Contracts F30602-94-C-0037, ARPA/SISTO contracts N00014-
91-J-1985, and N00014-92-C-0182 under subcontract KI-92-01-0182.
z Surface mail: Department of Computer Science, Duke University, Box 90129, Durham, N.C. 27708-0129.
extremely important. This has spurred a large interest in various aspects of out-of-core applications, including
language support, out-of-core compilers, parallel file systems, out-of-core algorithms, and out-of-core program
synthesis [2, 18, 7, 4].
Program synthesis (or automatic program generation) has a long history in computer science [16]. In the
recent past, tensor (Kronecker) product algebra has been successfully used to synthesize programs for the
class of block recursive algorithms for various architectures such as vector, shared memory and distributed
memory machines [11, 9, 5], and for memory hierarchies such as cache and single disk systems [14, 13]. We
have recently enhanced this program synthesis framework for multiple disk systems with the fixed physical
track data distribution [10] as captured by the two-level disk model proposed by Vitter and Shriver [19].
In this paper, we present a framework of using tensor products to synthesize programs for block recursive
algorithms for the striped Vitter and Shriver's two-level memory model which permits various block-cyclic
distributions of the out-of-core data on the disk array. The framework presented in this paper generalizes
the framework presented in [10]. We use the algebraic properties of the tensor products to capture the
semantics of block-cyclic data distributions cyclic(B), where B is the logical block size, on the disk array.
We investigate the implications of various block-cyclic distributions cyclic(B) on the performance of out-of-
core block recursive algorithms, such as the fast Fourier transform (FFT) and block matrix transposition
algorithm.
Tensor product representations of block recursive algorithms may involve stride permutations. Since a
stride permutation can be interpreted as a matrix transposition, synthesizing efficient out-of-core programs
for stride permutations is important [17, 6, 12]. We present a procedure for synthesizing efficient out-of-core
programs for stride permutations using a cyclic(B) distribution of the data. An algorithm for determining
the block size B which optimizes the performance of the synthesized programs is also presented.
We then discuss program generation techniques for tensor products with various block-cyclic data dis-
tributions. We discuss several strategies, such as factor grouping and data rearrangement, to improve the
performance for the tensor product formulas. We formalize the procedure of synthesizing efficient out-of-
core programs for the tensor product formulas with various data distributions as a dynamic programming
problem. However, since data rearrangement is too expensive in our target model as discussed in Section 7,
we have not incorporated it into our dynamic programming approach. In this sense, the stride permutations
should be mainly understood as a method of program synthesis for matrix transpositions not for the tensor
product formulas. We illustrate the effectiveness of this dynamic programming approach through an example
out-of-core FFT program.
We further examine the performance issues of synthesized programs. We show that:
1. The choice of data distribution has a large influence on the performance of the synthesized programs,
2. Our simple algorithm for selecting the appropriate data distribution size is very effective, and
3. The dynamic programming approach can always reduce the number of passes to access out-of-core
data.
The paper is organized as follows. Section 2 discusses formulation of block recursive algorithms using
tensor products and other matrix operations. In Section 3, we introduce a two-level computation model and
present the semantics of data distributions and data access patterns. We also argue the advantages of using
various block-cyclic distributions. Section 4 presents an overview of our approach of out-of-core program
synthesis. A generic program for synthesizing out-of-core programs for tensor products from given tensor
bases is also discussed. Section 5 presents a framework for synthesizing programs for the various block-cyclic
data distributions for stride permutations. In Section 6, we present our approach for synthesizing out-of-
core programs for tensor products. Section 7 presents a multi-step dynamic programming algorithm for
synthesizing programs for tensor product formulas. In Section 8, we summarize the performance results and
show the effectiveness of using various block-cyclic data distributions. Section 9 discusses related research.
Conclusions are provided in Section 10. In Appendix A, we list a selected set of symbols used in this paper.


Appendix

B and Appendix C present additional details for proving some of the properties discussed in
Section 4 and Section 5, respectively.
Product Algebra
The tensor (Kronecker) product [8] of an m \Theta n matrix A m;n and a p \Theta q matrix B p;q is a block matrix
A
obtained by replacing each element a i;j of A m;n by the matrix [a i;j B p;q ].
A
a
A tensor product involving an identity matrix can be implemented as parallel operation. An identity
matrix of order n is denoted as I n . Consider the application of (I
to a vector
n. This can be interpreted as n copies of B p;q acting in parallel on
disjoint segments of X nq . However, to interpret the application of
m;n\Omega I p ) to Y np as parallel operations
we need to understand stride permutations.
A stride permutation L mn
n is an mn \Theta mn permutation matrix. The application of L mn
n to X mn results
in a vector Y mn such that:
vector consisting of the elements in the set fX i+j \Lambdan
1)g. One important property of the stride permutation is L mp
p;q\Omega A m;n )L nq
q .
Using stride permutations, the application
m;n\Omega I p ) to the vector Y np can also be interpreted as p parallel
applications of A m;n to disjoint segments of Y np by using the identity L mp
p\Omega A m;n )L np
.
In this case, however, the inputs for each application of A m;n are accessed at a stride of p and the outputs
are also stored at a stride of p. In general, (I
m\Omega A
n;p\Omega I q ) can be interpreted as mq parallel applications of
A n;p .
The properties of tensor products can be used to transform the tensor product representation of an
algorithm into another equivalent form, which can take the advantage of the parallel operations discussed
above. For example, by using the following tensor product factorizations,
A
m;n\Omega I p )(I
m;n\Omega I q ); (4)
A\Omega B can be implemented by first applying q parallel applications of A and then m parallel applications of
properties of tensor products are listed below [11]:
1.
A\Omega (B\Omega
(A\Omega B)\Omega C;
2.
(A\Omega B)(C\Omega
AC\Omega BD; assume that the ordinary multiplications AC and BD are defined.
3.
4.
In contrast to the tensor products which can be used to describe various computations, the tensor
product of vector bases, called a tensor basis, can be used to describe the data access and storage patterns
of a multi-dimensional array. A vector basis e m
m, is a column vector of length m with a one at
position i and zeros elsewhere. We use e m
i to denote the ith index of a one-dimensional array of size m: Since
in+j , we can use e m
to represent the index [i; j] of a two-dimensional array if we assume a row-major
storage order of multi-dimensional arrays in the memory. In general, the tensor basis e m t
\Delta\Omega e m1
corresponds to index [i array. The indexing function needed
to access elements of a multi-dimensional array can be obtained by linearizing the tensor basis. For example,
linearizing the tensor basis e m t
\Delta\Omega e m1
results in the vector basis e m t \Delta\Delta\Deltam 1
. The index in
the linearized tensor basis is exactly the indexing function needed for accessing a t-dimensional array in the
row-major order. Equivalently, a vector basis e M
i can be factorized into a tensor product of vector bases
\Delta\Omega e m1
3). Factorization of a vector
basis corresponds to viewing a one-dimensional array as a multi-dimensional array. Using tensor bases, the
semantics of the stride permutation L mn
n can be formally expressed as: L mn
corresponds to matrix transposition of an m \Theta n array stored in the row-major order.
By appropriately factorizing the vector basis for an input vector, we can use the resulting tensor basis to
describe the data access pattern of a tensor product. For example, for a tensor product (I
m\Omega A
p;n\Omega I q ), the
input vector basis e mnq
s can be factorized to obtain the input tensor basis e m
j is known as
operator basis. The output vector basis e mnq
s can be factorized to obtain the output tensor basis e m
which can also be determined by the following identity,
m\Omega A
p;n\Omega I q )(e m
i\Omega A p;n e n
j\Omega I q e q
i\Omega A p;n e n
and by replacing A p;n e n
j by e p
Using these input and output bases, we can determine the input and output
data elements of an application of A p;n and derive a program with a two-dimensional iteration space:
We ignore the dimensions of matrices whenever they are clear from the context.
Code for (i; k)-th application of A p;n
The indices of the input data elements to the (i; k)-th application can be obtained from the linearized input
tensor basis e mnq
inq+qj+k as finq ng. Similarly, the output indices can be determined from
the linearized output tensor basis as fipq pg. Note that, there are no loops corresponding
to indices j and j 0 in the above program.
2.1 Tensor Product Formulation of Block Recursive Algorithms
A tensor product formulation of a block recursive algorithm has the following generic form:
Y
j\Omega A v
j\Omega I c j ); where A v j is a v j \Theta v j linear transformation, (6)
where
. Under this definition, a stride permutation can be understood as a special
case when are 1 \Theta 1 matrices, and A v j is a permutation matrix corresponding to a
stride permutation. The identity terms I r j allow decomposition of the computation into a set of smaller
size computations, which may be computed in the main memory. Similarly, the identity terms I c j allow
a decomposition of the computation into a set of subcomputations, where each subcomputation accesses
the data storage in a stride fashion. Although, these parallel and stride computational structures help in
decomposing the computation into smaller in-core computations, the task of combining these decompositions
with the goal of minimizing I/O for the entire computation is a challenging problem. We next present an
example to illustrate how to use tensor product formulas to represent the Cooley-Tukey FFT algorithm.
Fast Fourier Transform The Fourier transform can be denoted by the following matrix vector multiplications

where FN is a N \Theta N discrete Fourier matrix. FN (i;
and !N is the N-th
primitive root of unity
rs, then the discrete Fourier matrix FN can be factorized
as follows [11, 15],
r\Omega I s )T rs
rs
s is called twiddle factors and is defined as T rs
. By using the above factorization
recursively, we obtain the following tensor product representation of the Cooley-Tukey FFT algorithm,
Y
Y

Figure

1: The data organization for 4, and Each column is a disk. Each
box is a physical block. Each row consists of a physical track. The numbers in each box denote the record
indices.
2 i\Gamma1 is a diagonal matrix of constants and R 2 n permutes the input sequence to a bit-reversed order.
If we ignore the initial bit-reversal operation R 2 n and notice that I 2
2 i\Gamma1 is a diagonal matrix, then we
can see that the computational structure of F 2 n is captured by the first factor I 2
2\Omega I 2 i\Gamma1 . It can be
easily verified that this major computational portion is an example of Formula (6).
3 Parallel I/O Model with Block-Cyclic Data Distributions
We use a two-level model which is similar to Vitter and Shriver's two-level memory model [19]. However,
in our model the data on disks (called out-of-core data) can be distributed in different (logical) block sizes.
The model consists of a processor with an internal random access memory and a set of disks. The storage
capacity of each disk is assumed to be infinite. On each disk, the data is organized as physical block with
fixed size. Four parameters: N (the size of the input), M (the size of the internal memory), B d (the size of
each physical block), and D (the number of disks), are used in this model. We assume that M ! N ,
Bd .
In this model, disk I/O occurs in physical tracks (defined below) of size B d D. The physical blocks which
have the same relative positions on each disk constitute a physical track. The physical tracks are numbered
contiguously with the outermost track having the lowest address and the innermost track having the highest
address. The ith physical track is denoted by T i . Fig. 1 shows an example data layout with B
and Each parallel I/O operation can simultaneously access D physical blocks, one block from each
disk. Therefore parallelism in data access is at two levels: elements in one physical block are transferred
concurrently and D physical blocks can be transferred in one I/O operation. In this paper, we use the striped
disk access model in which physical blocks in one I/O operation come from the same track, as opposed to
the independent I/O model in which block can come from different tracks. We use the parallel primitives,
parallel read(i) and parallel write(i), to denote the read and write to the physical track T i , respectively. We
define the measure of I/O performance as the number of parallel I/Os required.
3.1 Block-Cyclic Data Distributions
Block-cyclic distributions have been used for distributing arrays among processors on a multiprocessor sys-
tem. A block-cyclic distribution partitions an array into equal sized block of consecutive elements and then
maps them onto the processors in a cyclic manner. If we regard the disks in the above model as processors,
D D D D
43
44
48

Figure

2: The data organization for 8. Each column is a disk.
The first left shadowed box denotes an example logical block. There are two logical tracks LT 0 and LT 1 each
of them consists of two physical tracks.
then the data organization described above (e.g. in Fig. 1) is exactly a block-cyclic distribution (denoted as
with the block size B d .
Moreover, we can assume that data can be distributed with an arbitrary block size 2 . Fig. 2 shows the
data organization for the same parameters as in Fig. 1, but with a cyclic(8) distribution. Notice that the
size of the physical track and the size of the physical block are not changed. However, they contain different
records. We will call B records in a block formed by a cyclic(B) distribution as a logical block. Similarly,
the logical blocks which have the same relative positions on each disk consist of a logical track. The ith
logical track is denoted as LT i . Note that each parallel I/O operation still accesses a physical track not
a logical track. Hence, several parallel I/O operations are needed to access a logical track. For example,
to load the logical track LT 1 in Fig. 2, two parallel read operations parallel read(2) and parallel read(3),
which respectively load the physical tracks T 2 and T 3 , are needed. We next use a simple example to show the
advantages of using logical distributions on developing I/O-efficient programs for block recursive algorithms.
Why Logical Data Distributions? Assume that we want to implement F
I 8 on our target model
under the parameters given in Fig. 1. Further, we assume that the size of the main memory is the half of the
size of the inputs. Because we are mainly interested in data access patterns, we ignore the real computations
conducted by F 8 . The only thing we need to remember is that F 8 needs eight elements with a stride of eight
because of the existence of the identity matrix I 8 .
We first consider implementing F
I 8 on the physical block distribution. From the above discussion,
we know that the first F 8 needs to be applied to eight elements: 0; 8; 16; 24; 28; 32; 40; 48, and 56. From
Fig. 1, we can see that these elements required by the F 8 computation are stored on four physical tracks.
However, our main memory can hold only two physical tracks, so that we can not simply load all of the four
physical tracks into the main memory and accomplish the computation in one pass of I/O. To get around
this memory limitation, we can use two different approaches.
First, we load the first physical track and keep the first half of the records in each physical block in that
loaded physical track and throw other half of the records. We do this for every other physical track. Then
we do the computation for half of the records in the main memory. After finishing computation for half of
the records, we write the results out. Then we repeat the above procedure. However, we now keep other half
2 Cormen has called this data organization on disks as a banded data layout [3] and studied the performance for a class of
permutations and several other basic primitives of NESL language[1].
of the records in the main memory for each loaded track. By doing computation in this way, it is obviously
that we need two passes to load out-of-core data.
Another method is that we use a logical block distribution. Suppose that the size of a logical block is
eight as shown in Fig. 2. Now, the eight records required by one F 8 are stored on two physical tracks, physical
track one and three, or physical track two and four. Therefore, if we can load physical tracks one and three
first and do the computation, then load physical track two and four and do the computation, then we can
finish the computation in one pass. This example clearly shows the advantages of using logical distributions
comparing with using only physical track distributions. However, there are several problems which we have
not addressed here, such as how to determine the block size of the logical distribution and how to determine
the data access patterns. We will discuss these issues in the rest of the paper.
Note that the striped I/O model with the data distribution cyclic(B d ) can be transformed to a model
with
1. However, this transformation may not be true for an arbitrary cyclic(B)
distribution. Therefore we normally can not reduce the problem to the simpler case, where
1. For simplicity, we make the following assumptions. The input and the output data are stored in
separate set of disks. All parameters are power of two 3 . The block size B of the distribution is a multiple
of B d .
3.2 Semantics of Data Distributions and Access Patterns
As discussed in [9], a block-cyclic distribution can be algebraically represented by a tensor basis. That
approach can be adopted to the disk model by substituting disks for processors. However, because of
the existence of the physical blocks and the physical tracks, the tensor basis used to define a block-cyclic
distribution for multiprocessors needs to be further factorized. We call this factorized tensor basis as a
(out-of-core) data distribution basis, which is defined as follows:
Definition 3.1 . If a vector N , with distributed according to the cyclic(B)
distribution on D disks, then its data distribution basis is defined as:
g\Omega e D
b\Omega e Bd
b d
We use D(s) to refer to the sth factor (from the left), e.g.,
d .
For example, the data distribution basis for Figure 2 is e 2
b d
. The data distribution basis
for

Figure

1 can be written as e 4
b d
selected portion of the distribution basis in
Formula (10) can be used to obtain the indexing function needed to denote a particular data unit such as a
logical track or a physical track. Let,
3 The results can be easily generalized to all parameters to be power of any integer.

Figure

3: An example data organization and access pattern. Each logical track consists of two physical
tracks. Each physical block is further decomposed into two sub-blocks. All of those shadowed sub-blocks
together can be used to form a data access pattern. One example is to access the shadowed sub-blocks in
the row-major order.
Then the indexing function for accessing the physical tracks can be obtained by linearizing the tensor
basis obtained from physical-track-basis(D). By taking the difference 4 of the data distribution basis
with each of them, we can have tensor bases which denote the records inside a logical track and a physical
track, respectively. These tensor bases are called a logical track-element basis (e D
b\Omega e Bd
b d
) and a physical
track-element basis (e D
d\Omega e Bd
b d
respectively. This concept can be further extended to denote other subsets of
data which we will discuss later.
On the other hand, different orders of instantiating the indices in the data distribution basis defined in
Formula (10) result in different access patterns of the out-of-core data. For example, if we instantiate the
indices in that formula in the order from left to right, i.e. g is the slowest and b d is the fastest changing
index, then we actually access data first in the first logical block in the first disk and then access the first
logical block in the second disk. After finishing the access to the first logical track sequentially, the second
logical track is accessed, and so on. It is obvious that if the order of the indices in a data distribution basis
is instantiated differently, then the out-of-core data will be accessed in a different pattern. For example,
if we instantiate the index b b in e B b
before the index d in e D
d in Formula (10), then it results in an access
pattern where in first the data along a physical track is accessed and then to the successive physical tracks
are accessed. This change in the instantiation order of the indices can be regarded as a permutation 5 of
the data distribution basis. We will call a permutation of a data distribution basis as a loop basis. In the
synthesized programs, every index in a loop basis may be used to generate a loop nest. The order of the
loop nests is determined by the order of the vector bases in the loop basis.
More generally, the out-of-core data can be viewed as organized as a multi-dimensional structure. For
example, the data layouts in Fig. 2 can be viewed as a three-dimensional array, where each logical block is
viewed as a B b \Theta B d matrix. Further we can convert B d to a two-dimensional structure. We can combine
records in some of the disks together as a submatrix. These other views can be denoted by factorizing and
regrouping data distribution bases and can be used to form different data access patterns. Fig. 3 shows an
example data organization and access pattern. We assume that each logical
4 Let S and G be two tensor bases. Their difference is denoted as S-G and is a tensor basis which is constructed from deleting
all the vector bases in G from S.
5 Let S be a tensor basis and S
=\Omega q
. Let ff be a permutation on a permutation of S is a tensor basis
defined as follows, ff(S)
=\Omega q
track consists of two physical tracks. We further decompose each physical block B d as two sub-blocks, each
of them has size B d1 . In order to reflect this data organization, we can factorize the data distribution basis
as follows,
g\Omega e D
b d
After we permute this factorized data distribution basis, we can have the loop basis which accesses subsets
of data in different patterns. Assume that we want to access data in Fig. 3 in the following order: first
accessing the darker shadowed sub-blocks in the row-major order and then accessing the lighter shadowed
sub-blocks in the row-major order. Then we can move D(3) and D(4) before D(2) and D(1), respectively.
This results in the following loop basis,
b d
b\Omega e D
We can verify the correctness of this loop basis as follows. (1) The index b d2 chooses the same shadowed
sub-blocks. (2) The index g chooses the logical tracks. (3) The index b b chooses the logical blocks. (4) The
index d chooses disks. (5) The index b d1 chooses the records inside a sub-block. Since the increasing of the
indices is in the reversed order of the above five steps. We first access records inside the first darker shadowed
sub-block in the first disk and then access the records in the same shadowed sub-block in the second disk,
and so on. After accessing all the records in the same darker shadowed sub-blocks in the row-major order,
we repeat the procedure for the lighter shadowed sub-blocks.
4 Overview of Program Synthesis
Efficient implementations of block recursive algorithms are obtained by using the properties of tensor products
to transform the tensor product representations for the block recursive algorithms. The transformations
will use the results of the performance of each individual tensor product. This performance of each tensor
product is obtained by using the method presented in Section 6. Fig. 4 shows the procedure of synthesizing
efficient out-of-core programs for block recursive algorithms. In that figure, the augmented tensor basis contains
the following information: data distribution bases, loop bases, sub-computations and memory-loads,
which are necessary for code generation and will be explained further in this section.
If the input tensor product formula consists of only a stride permutation, it will quickly pass through the
program transformation step. Then it will use an algorithm presented in Section 4 to generate an augmented
tensor basis. For other tensor product formulas, the successive steps will be applied to each tensor product
in the transformed tensor product formulas. To obtain an augmented tensor basis for each tensor product,
we use the method presented in Section 6. The code generation step for both tensor products and stride
permutations will use the same procedure presented in this section.
Our presentation of deriving efficient implementations for the block recursive algorithms is in a reversed
order of Fig. 4. We first present a procedure for code generation by using the information contained in the
Target machine model
Tensor Product Formula
Tensor Product Formula + Data Distribution
Augmented Tensor Basis
Code Generation
Parallel I/O Program
Program Transformation
Computation Partitioning
Access Pattern Analysis

Figure

4: The procedure of synthesizing efficient out-of-core programs for block recursive algorithms.
augmented tensor basis. Then, we determine efficient implementations for a stride permutation and a simple
tensor product with a given data distribution on a given model by determining the corresponding augmented
tensor bases. We develop an algorithm to determine the data distribution which can result in an efficient
implementation for a simple tensor product. Using the information obtained so far, we use a dynamic (or a
multi-step dynamic) programming algorithm to determine an efficient implementation for the block recursive
algorithms.
In the rest of this section, we summarize how to synthesize efficient programs for simple tensor products
and stride permutations on our target machine model. The further details and the performance of synthesized
programs are then discussed in the next two sections.
To minimize the number of I/O operations for a synthesized program, we need to exploit locality by
reusing the loaded data. This requires decomposing the computation and reorganizing data and data access
patterns to maximize data reuse. In the synthesized program, the same sub-computation is performed several
times over different data sets. Hence, the loop structure of the synthesized programs is constructed as follows.
An outer loop nest enclosing three inner loop nests: read loop nest for reading in the data, computation loop
nest for performing sub-computation on the loaded data, and write loop nest for writing the output data back
to the disk. The inner read loop nest should load the out-of-core data without overflowing the main memory.
We refer to each of this data sets as a memory-load. The inner computation loop nest should perform
subcomputation on a memory-load. Further, the data sets should be accessed using parallel primitives,
parallel read and parallel write, to load or store a physical track each time.
One of the main results in this paper is that an efficient tensor product decomposition of a computation
Generate loops for indices in -n
Generate loops for indices in -m
Parallel read using the input distribution basis
Construct a memory-load
End the loops corresponding to -m
Perform operations to a memory-load
Generate loops for indices in ' m
Parallel write using the output distribution basis
End the loops corresponding to ' m
End loops corresponding to -n

Figure

5: A procedure of code generation for a tensor product.
and an efficient data access pattern can be obtained by using the algebraic properties of data distribution
bases and loop bases. In other words, once we have determined data distribution bases and loop bases,
we can determine the memory-loads and the subcomputations (or operations) for each memory-load. For
a tensor product computation, the input and output data may be organized and accessed differently, we
therefore use input data distribution basis fi, output data distribution basis ffi, input loop basis -, and output
loop basis ' to denote them respectively. Data distribution bases can be obtained from input and output
bases by rewriting them into the form of Formula (10). However, it is a non-trivial task to determine the
loop bases with the goal of minimizing the number of I/O operations. We next present a generic synthesized
program and then summarize some general ideas for determining loop bases. We will discuss the details on
how to determine loop bases and therefore memory-loads and operations to each memory-load in the next
two sections.
Consider the task of generating target code assuming that the data distribution bases, loop bases,
memory-loads and operations to each memory-load have been already determined. As we discussed in
the previous section, under our striped I/O model, each I/O operation will read or store all the records in a
physical track each time. Hence, only part of loop basis will explicitly appear in the synthesized programs.
Moreover, the input (or output) loop basis can be separated into two parts such that the first part specifies
memory-loads and the second part specifies the records inside a memory-load. The second part can be
further separated into two parts: one part, denoted as -m , is used to construct a memory-load and another
part, denoted as - will not generate loop nests in the synthesized programs. In other words, we can write
the input and the output bases as follows:
where, we call -n (or ' n ) as a memory basis, since each instantiation of the indices in -n corresponds to a
memory-load. Using these loop bases, a generic program can then be obtained as described in Fig. 5. Notice
when parallel read a track, the track number is obtained from the indexing function of physical-track-basis(fi)
(part of the input data distribution basis) as defined in Formula (12). So does the parallel write. Fig. 6
shows an example synthesized program for I
I 4 . We assume that
F 2 a 2 \Theta 2 matrix, and data are distributed in cyclic(2) manner. It uses e 8
b d
as both the input and
the output distribution bases. The input and the output loop bases are also the same as e 2
b d
g1 is a factorization of e 8
.
Let us further examine Formulas (15) and (16). Obviously, if - and ' - consist of the physical track-
element bases for input and output data respectively, then the out-of-core data needs to be accessed only
once. In terms of memory-loads, each memory-load has the following properties. The input for each memory-
load occupies all of the locations in a set of physical tracks specified by the input data distribution basis.
And after computing these records in the main memory, they are organized to occupy all of the locations in
a set of physical tracks specified by the output data distribution basis. We call this type of memory-load as
a perfect memory-load. If we can construct memory-loads in this manner, then we can synthesize a program
which accesses out-of-core data only once (called a one-pass program).
However, it may not be possible to construct the perfect memory-loads for some computations. In that
case, we may need to keep only part of the records from a loaded physical track in the main memory and
discard the other records. Therefore, a multi-pass program needs to be synthesized in which the same physical
track is loaded several times. In terms of tensor bases, this corresponds to moving some of the vector bases
in the physical track basis out of the memory basis. Consider the example presented in Section 3.2, where we
moved fi(4) as the first factor. Since the unit of data access is still a physical track, this moving corresponds
to loading the same track Bd
times. However, for each loaded track, only half of the records can be kept in
the main memory if we assume that the size of the main memory is the half of the input size. The records
which should be kept for each loaded track can be determined by the vector bases which are moved before
the memory-load basis. In this case, it is e
. Each instantiation of the index b d2 determines the sub-blocks
in a physical track which should be kept for the current memory-load. The detailed program for loading
out-of-core data and constructing memory-loads is shown below,
ENDDO ENDDO ENDDO ENDDO ENDDO
where A is a temporary array for holding a physical track and X holds a memory-load.
In summary, in order to determine efficient loop bases, we construct initial loop bases - and ' such that
- and ' - consists of the physical track-element bases from the input and the output data distribution bases,
respectively. We then determine which of the vector bases, if any, need to be moved from - into -n . These
moved vector bases are used to determine which portions of a physical block should be kept for the current
memory-load. The size of these moved vector bases is equal to the number of times of the same physical
tracks to be loaded. Further, we may need to determine the order of the rest vector bases in - \Gamma - and
to reflect the order of accessing physical tracks.
Parallel read from a track
ENDDO
// Perform operations for a memory load
// Write the result back
Parallel write to a track
parallel
ENDDO ENDDO

Figure

Code for the tensor product I
2\Omega I 4 , where X is an array of size M and A =
I
I 4 .
Synthesizing Programs for Stride Permutations
In this section, we present a framework for synthesizing efficient out-of-core programs for stride permutations
using a cyclic(B) distribution. The performance of synthesized programs will be represented as a function
of the size of a sub-tensor basis, whose value can be obtained when the distribution size is given. We also
present an algorithm to determine the distribution which will optimize the performance.
5.1 Stride Permutations in Cyclic(B) Distribution
As we have mentioned, our goal is to decompose computations into a sequence of sub-computations operated
on perfect memory-loads. However, this may not be always possible because of the limited memory size. In
that case, we minimize the number of times the data is loaded for each memory-load as well as we ensure
that each physical track of the output is written only once in parallel. We will develop an approach to
determine the input and output loop bases for the given distribution cyclic(B). Based on these loop bases
and data distribution bases, we determine memory-loads and operations to the memory-loads. Following
which a program can be synthesized by using the procedure presented in Section 4. The cost of the program
can also be determined from the loop bases. We summarize our results as the following theorem and then
present a constructive proof.
Theorem 5.1 Let
Y be input and output vectors with length N ,
respectively. Let X and Y be distributed according to cyclic(B) and the data distribution bases be denoted
as fi and ffi, respectively. Further denote that
fi(2)\Omega fi(4) and '
ffi(2)\Omega ffi(4). Then a program can
be synthesized with N
operations for the stride permutation
Proof: We present an algorithm as shown in Fig. 7 for determining the input and the output loop bases. The
algorithm is further explained in Step 1 as shown below. In Step 2 and Step 3, we show how to construct
6 The notation j S j denotes the size of the tensor basis S, which is equal to the multiplication of the dimensions of each
vector basis in S.
Initialization
fi(1)\Omega fi(3)\Omega fi(2)\Omega fi(4)
// One-pass or multi-pass implementation
if (j ('
else
consists of the last factors of
the factorized tensor basis and
BdD .
// The final input and output loop bases

Figure

7: An algorithm for determining input and output loop bases.
memory-loads and operations for a memory-load. In Step 4, we show that I/O costs can be obtained from
those information.
1. Determine input and output loop bases. We begin with the following construction for the input
and the output loop bases,
where we use the convention that - appearing on the right hand side refers to the original representation,
which is equal to
fi(1)\Omega fi(3)\Omega fi(2)\Omega fi(4), and - appearing on the left hand side refers to an update. So
does '. Further, we assume that
ffi(2)\Omega ffi(4). It is easy to verify that (' - \Gamma-
is a permutation of (-
Therefore, they denote the same records. Thus, if the number of
records denoted by j
than the size of the main memory, then we can simply take
. However, the number of records denoted by j
may
exceed the size of the main memory. In that case, we want to construct memory-loads which can be
obtained by reading the input data several times however writing the output data only once. In terms
of tensor bases, as we discussed in Section 4, this reloading can be achieved by looping over part of
the indices in - . In other words, we need to factorize - as -2 and -1 such that -2 denotes which
sub-blocks should be kept for a loaded physical track and -1 denotes records inside each sub-block.
Further, j -2 j is equal to the number of times we will reload each physical track. This reloading is
achieved by taking moving -2 before - ' m . In summary, the input and output loop
bases in Formulas (17) and (18) are modified as follows:
consists of the last factors of the factorized tensor basis and the
size of ' m is equal to M
BdD .
ffl For input loop basis. Let -2
Therefore, now the input and output loop bases can be written as,
We further verify the following facts. First,
and '
contain the same vector bases, however, in the different order. (The proof is presented in


Appendix

B.) Therefore, they denote the same records, however, in different order. Second, from the
previous results, we have that j
. Therefore the records denoted by them
can fit into a memory-load. Third, j -m
DBd ), which means that we have loaded more
records which can fit into the main memory and we need to discard some of the records. The details
for determining which records to be discarded will be discussed in the next step. (4) -n and ' m contain
the same vector bases. We therefore can set that ' m=-n , which will only change the order of writing
results onto physical tracks.
2. Determine a memory-load. When j ('
Therefore, the records denoted by them can be used to form a perfect memory-load. However, when
this condition is not satisfied, we need to use Formula (19) and (20) as the input and output loop
bases, respectively. Because j
, the size of each memory-load can be
set to be equal to the size of the main memory. However, as we mentioned before, we need to discard
some records from each loaded track to form the memory-load. This can be done by linearize -2 .
Each instantiation of -2 will give a set of sub-blocks in a physical track which should be kept.
3. Determine operations for a memory-load. As we mentioned above, for each memory load, the
tensor vectors in the input and output loop bases which denote the records inside a memory-load are
the same, but in a different order. In other words, one is a permutation of another. Because the input
and output loop bases are permutations of the input and output data distribution bases, we actually
permute a memory-load of data each time. Therefore, each in-memory operation is nothing more than
a permutation for a subset of data distribution bases denoted by
-m\Omega
4. I/O cost of synthesized programs. It is readily to see that if j
program can be synthesized, i.e., the number of parallel I/Os is 2N
BdD . When the above condition does
not hold, we keep j records for each loaded physical track and load the same physical track j -2 j
times. Moreover, since
DBd , it can be easily determined that j -2
. Because we
write out each record only once, the number of parallel I/O operations is (1
BdD . Combining
these two cases together, we yield the performance results presented in the theorem. Further, a
program with this performance can be synthesized by using the procedure listed in Fig. 5.
5.2 Determining Efficient Data Distributions
In the previous subsection, we presented an approach for synthesizing efficient I/O programs for the given
data distribution. We now present an algorithm to determine the data distribution which optimizes the
performance of the synthesized program. The idea of the algorithm is as follows. We begin with the
physical track distribution cyclic(B d ), i.e., initially . If there is a one-pass algorithm under this
distribution, then B d is the desired block size for the data distribution. Otherwise, we double the value
of B. If the performance of the synthesized program under this distribution increases, we continue this
procedure. Otherwise, the algorithm stops and the current block size is the desired size of data distributions.
We formalize this idea in Fig. 8.
number of I/Os when using cyclic(B)
while (Cost 6= 2N
DBd and B - N
D ) do
number of I/Os when using cyclic(B)
If C new - Cost then Cost = C new else break
output distribution size = B=2, number of

Figure

8: Algorithm for computing the desired size of data distributions.
6 Synthesizing Programs for Tensor Products
For the tensor product, I
R\Omega A
V\Omega I C , the main computation matrix A V needs V records with the stride C.
We call each of these V records as a desired record. We first present a possible general form of the input
and output loop bases for the given distribution cyclic(B). The parameters in this form can be determined
by analyzing the relative values of the parameters. Based on these loop bases and also data distribution
bases, we can determine memory-loads and operations to each memory-load. Therefore a program can be
generated by using the procedure discussed in Section 4. The cost of the program can also be determined
from the loop bases.
Since the tensor product I
R\Omega A
V\Omega I C does not change the order of the inputs (or it can be computed
in-place), we will use the same input and output data distribution bases for the input and output data and
also the same input and output loop bases for programs synthesized in this section. Therefore, we will only
consider input, input distribution and input loop bases and assume that all of the output bases are the same
as their input bases respectively. We summarize our results as a theorem and then present a constructive
proof.
Theorem 6.2 Let the input data be distributed according to cyclic(B) and the input data distribution basis
be denoted as -. Let
fi(2)\Omega fi(4). Further assume that -1 denotes a subset of - and
is moved into the memory basis. Then for the tensor product I
R\Omega A
V\Omega I C , where RV
a program can be synthesized with 2N
BdD parallel I/O operations; otherwise a program can be
synthesized with j -2 j 3N
BdD parallel I/O operations.
Proof:
1. Determine input loop basis. If the desired records for an A V computation are stored in t physical
tracks and t - M
BdD , then we can simply load t tracks in parallel and therefore a one-pass program
can be generated. However, when t ? M
BdD , we can not keep all of the records in t tracks in the
main memory. We take a simple approach that we keep as many as possible records which follow
each desired records in the same track in the main memory. Then we reload these tracks to finish the
computations for the other records. In terms of the tensor basis, we need nothing more than factorizing
and permuting the input data distribution basis to reflect these data access patterns.
More specifically, we begin with
-(2)\Omega -(4), and
defined as the
same initial value as defined in Section 5. For a one-pass program, we factorize and permute
-n\Omega -m to
change the order of accessing physical tracks. However, for a multi-pass program, we need to factorize
and permute all of the -, since we need to keep part of the records loaded in the main memory and
discard the other records. The part of records to be kept or discarded can be denoted by a subset of
the vector bases in the physical track basis. One of the example factorizing and permuting and its
semantics has been discussed in Chapter 4. However, in general, in order to factorize and permute a
tensor basis to a desired form, we need to examine the relative values of the parameters in the targeted
I/O model, the tensor product and the size B of the data distribution. Because of the space limitation,
we present the major idea of the analysis in Appendix C. For the following analysis, it is enough to
say that we have found subsets of - , denoted as -1 and -2 . -2 is moved into the memory basis
and will generate loop nests for data access.
2. Determine a memory-load. For a one-pass program, we can simply factorize
-n\Omega -m and
BdD . For a multi-pass program, we to be
-n\Omega -m such that j -m
and all of the vector bases in -2 appear in -n . Moreover, for a multi-pass program, as discussed in
Section 5, we use -2 to determine which records should be kept for the current memory-load.
3. Determine operations for a memory-load. The original tensor product can be regarded as R
parallel applications of A V to the inputs with stride C. When data are distributed among disks and
loaded in units of physical tracks, the net effect is to possibly reduce the stride with which each A V will
access in the main memory. The operations to a memory-load have a general form of I M
Z\Omega A
V\Omega I Z .
However, the value of Z will depend on the relative values of the parameters. Appendix C presents
the major ideas on how to determine the value of Z.
4. I/O cost of synthesized programs. For a one-pass program which does not move any vector bases
in - , the number of parallel I/Os is simply equal to 2N
BdD . In other words, the synthesized program is
optimal in terms of the number of I/Os. For a multi-pass program, we need to read the inputs j -2 j
times. Therefore the number of parallel I/O operations is j -2 j 3N
BdD . The constant 3 can be explained
as follows. When we store a physical track, we need to read that physical track into the main memory
again, since part of the records in that physical track have been discarded. By reloading this physical
track, we can reassemble the physical track with the part of updated records and then write it out in
parallel. Otherwise, part of the records to be written out in that physical track may not be correct.
Further, "reassembling" the physical track needs to use the tensor basis ' -2 (notice that ' -2 is equal
to -2 ) to put the updated records into the correct locations of the physical track. This is similar with
how to use -2 to take sub-blocks out from a loaded physical track for the current memory-load.
Now, a program with the performance discussed above can be synthesized by using the procedure listed
in Fig. 5. However, to be accurate, when synthesizing a multi-pass program, we need to incorporate
the idea of "reassembling" a physical track into the write-out part of the procedure listed in Fig. 5,
which, as we discussed above, is nothing more than using the linearization of ' -2 to put sub-blocks in
the current memory-load into the correct locations of the reloaded physical track.7 Synthesizing Programs for Tensor Product Formulas
In this section, we discuss techniques of program synthesis for tensor product formulas. There are several
strategies for developing disk-efficient programs, such as exploiting locality and exploiting parallelism in
accessing the data. Similar ideas have been discussed in [13], where they use factor grouping to exploit
locality and data rearrangement to reduce the cost of I/O operations. We have also presented a greedy
method which uses factor grouping to improve the performance for the striped Vitter and Shriver's two-level
memory model with a fixed block size of data distribution [10].
Factor grouping combines contiguous tensor products in a tensor product formula together and therefore
reduces the number of passes to access the secondary storage. Consider the core Cooley-Tukey FFT
computation computation can be represented by Formula (9) by ignoring the
initial bit-reversal and the twiddle factor operations). For i=2 and 3, we have the following tensor products
I 2
2\Omega I 2 , and I 2
respectively. Assuming that each of these tensor products can be
implemented optimally, the number of parallel I/O operations required to implement these two steps individually
is 4N
DB . However, they are successive tensor products in Formula (6). Hence, by using the properties
of tensor products, they can be combined into one tensor product, I 2
2\Omega I 2 , which may also be
implementable optimally by using only 2N
DBd parallel I/O operations.
Data rearrangement uses the properties of tensor products to change the data access pattern. For
example, the tensor product I
R\Omega A
I C can be transformed into the equivalent form (I
C ). In the best case, the number of parallel I/Os required is 6N
DBd after using this transformation,
since at least three passes are needed for the transformed form. Because of the extra passes introduced by
this transformation, it is not profitable to use it for our targeted machine model. Further, the first and
the last terms in the transformed formula may not be implementable optimally. Therefore, we have not
incorporated this transformation into our current optimization procedures.
Minimizing I/O Cost by using Dynamic Programming Since factor grouping (as shown above)
and the size of the data distribution (as will be shown in the next section) have a large influence on the
performance of synthesized programs, we take the following approach for determining an optimal manner
in which a tensor product formula can be implemented. We use the algorithm for determining the optimal
data distribution presented in Fig. 8 as a main routine. However, for each cyclic(B) data distribution, we
use a dynamic programming algorithm to determine the optimal factor groupings. Hence, we also call this
method as a multi-step dynamic programming method.
Let C[i; j] be the optimal cost (the minimum number of I/O passes required to access the out-of-core
data) for computing (j \Gamma i) tensor factors from the ith factor to the jth factor in a tensor product formula.
Then C[i; j] can be computed as follows:
ae
In the above formula, C 0 denotes the cost for computing a tensor product. The method of determining
the cost of a tensor product has been discussed in Section 6. The values of C 0 for different cases can be
found in Table 3 and Table 4 presented in Section 8.2. A special case of needs to be further explained.
we assume that C[j and we use C[i; k] to represent the cost of grouping all the
tensor product factors from i to j together. Because the grouped tensor product is a simple tensor product,
the value of C[i; k] in this case can also be determined by using Table 3 and Table 4 presented in Section 8.2.
However, in this case, if or the size of grouped operations is larger than the size of the main
memory, we don't want to group all of k \Gamma i factors together. We assign a large value such as 1 to C[k; j]
to avoid it from being selected.
Performance Results
8.1 Performance of Synthesized Programs for Matrix Transposition
Given the flexibility of choosing different data distributions, we can synthesize programs with better performance
than those obtained using a fixed size of data distributions for stride permutations. We present
a set of experimental results for the number of I/O operations required by the cyclic(B d ) distribution and
distribution, where the size B of the distribution varies. These results are summarized in Table 1,
and

Table

2. From the tables, we can see that the number of passes is not a monotonically increasing or
decreasing function. However, it normally decreases and then increases as B is increased. Therefore the
algorithm in Fig. 8 has a good chance to find out an efficient size of data distributions. We also notice that
for any stride permutation, we can always find out a distribution which can implement the computation in
only one-pass.

Table

1: The number of I/O passes required for performing stride permutation L PQ
using various cyclic(B)
distributions (D

Table

2: The number of I/O passes required for performing stride permutation L PQ
using various cyclic(B)
distributions (D
8.2 Performance of Synthesized Programs for Tensor Products
The number of I/O passes required by the synthesized programs for a tensor product are summarized in

Table

3 and Table 4 by going through various cases using the approach presented in Appendix C. We
can verify that the results presented here are more comprehensive than the results presented in [10]. In
most cases, using the approach presented in Section 5.1, we can actually synthesize programs with better
performance. For example, when , from [10], a program with V BdD
passes will be synthesized. However, for those conditions, we can have that C ? B d , and V C ? M . If
we further assume that D, then from the results in Table 3, we can synthesize a program with V C
passes, which is less than V BdD
M .

Table

3: The number of I/O passes required for the tensor product I
R\Omega A
V\Omega I C using a cyclic(B) distribution,
BdD ) is the maximum number of the physical tracks in a memory-load.

Table

4: The number of I/O passes required for the tensor product I
R\Omega A
V\Omega I C with given size B of data
distributions, where N t (= M
BdD ) is the maximum number of the physical tracks in a memory-load.
We now show that by using an appropriate cyclic(B) data distribution, a better performance program
can be synthesized for most of the cases. Several typical examples are shown in Table 5. We notice that
when we increase B, we can reduce the number of passes of data access for most of the cases and the decrease
in the number of passes can be as large as eight times. The values in the table also suggest that we can use
the algorithm presented in Fig 8 to find out an efficient size of data distributions for a given tensor product.
We also notice that for some cases, such as C - B d , we can not improve the performance. The reason is that
the stride required by A V is less than the physical block size, we can not reduce it further by redistribution.

Table

5: The number of I/O passes required for various sizes of data distributions and the tensor product
I
R\Omega A
V\Omega I C . Let . Let the size N(= RV C) of the input vector be very large.
8.3 Performance of Synthesized Programs for Tensor Product Formulas
We show the effectiveness of the multi-step dynamic programming method by comparing the programs
synthesized by it with the programs synthesized by the greedy method and the dynamic programming
method (applied to the data distribution of fixed size), respectively. The example we use is the core Cooley-Tukey
FFT computation. The results for several typical sizes of inputs are shown in Table 6. We find out
that using dynamic programming for the fixed size of the cyclic(B d ) distribution normally can not improve
the performance over the greedy method. However, by using the multi-step dynamic programming, we can
reduce the number of passes for the synthesized programs by at least 1 if N is significantly large. Because
the input size is large, the performance gain by reducing even one pass to access the out-of-core data is
significant.
Greedy
D.P.
M.D.P. 4

Table

The number of I/O passes for the synthesized programs using Greedy, Dynamic programming
(D.P), and Multi-step dynamic programming(M.D.P) methods (D
9 Related Research
Tensor product algebra has been successfully used for synthesizing programs for block recursive algorithms for
various architectures such as vector, shared memory, and distributed memory machines [11, 9, 5]. Recently,
the tensor product framework is used for synthesizing programs for memory hierarchies. For example, a
method of program synthesis for a single disk system is discussed in [13]. However they have not addressed
the issues of data distributions on multiple disk systems. In [14], Kumar, Huang, Sadayappan and Johnson
discussed the method of program synthesis for cache memory, where they addressed the issue of data layouts
on a set-associated cache. In [10], we presented a framework of using tensor products for synthesizing efficient
programs for a deeper level of memory hierarchy modeled by Vitter and Shriver's two-level memory model.
However, we only considered the data distributions to be fixed physical track distributions. The programs
synthesized are also not as efficient as the programs synthesized from the approach presented in this paper.
There are also many other recent research efforts in other areas of I/O intensive applications, which include
out-of-core algorithms, languages, compilers, parallel file systems, and performance models [2]. For example,
Vitter and Shriver proved lower and upper bounds for matrix transposition and FFT graph computations for
the two-level memory model [19]. Cormen presented algorithms for BMMC permutations [3], which includes
stride permutations as a sub-class, on the two-level memory model.
Conclusions
We have presented a novel framework for synthesizing out-of-core programs for block recursive algorithms
using the algebraic properties of tensor products. We use the striped Vitter and Shriver's two level memory
model as our target machine model. However, instead of using the simpler physical track distribution normally
used by this model, we use various block-cyclic distributions supported by the High Performance Fortran to
organize data on disks. Moreover, we use tensor bases as a tool to capture the semantics of data distributions
and data access patterns. We show that by using the algebraic properties of tensor products, we can
decompose computations and arrange data access patterns to generate out-of-core programs automatically.
We demonstrate the importance of choosing the appropriate data distribution for the efficient out-of-core
implementations through a set of experiments. The experimental results also shows that our simple algorithm
for choosing the efficient data distribution is very effective. From the observations about the importance of
data distributions and factor grouping for tensor products, we propose a dynamic programming approach
to determine the efficient data distribution and the factor grouping. For an example FFT computation, this
dynamic programming approach can reduce the number of I/O passes by at least one comparing with using
a simpler greedy algorithm.

Acknowledgements

We thank Peter Mills for his comments on this paper.



--R

Vector models for data-parallel computing

Virtual Memory for Data-Parallel Computing
Integrating theory and practice in parallel file systems.
EXTENT: A portable programming environment for designing and implementing high performance block recursive algorithms.
A fast computer method for matrix transposing.
Parallel I/O systems and interfaces for parallel computers.
Kronecker Products and Matrix Calculus: With Applications.
Synthesizing Communication-efficient Distributed-Memory Parallel Programs for Block Recursive Algorithms
Generating efficient programs for two-level memories from tensor products
A methodology for designing
Efficient transposition algorithms for large matrices.
A methodology for generating efficient disk-based algorithms from tensor product formulas
An algebraic approach to cache memory characterization for block recursive algorithms.
Computational frameworks for the fast Fourier transform.
Parallel Algorithm Derivation and Program Transformation.
Parallel processing with the perfect shuffle.
Compilation of out-of-core data parallel programs for distributed memory machines
Algorithms for parallel memory I: Two-level memories
--TR
tensor product;block recursive algorithm;data distribution;program synthesis;fast Fourier transform;Parallel I/O