Technical Report 2012-001

Co-design of Extremely Scalable
Algorithms/Architecture for 3-Dimensional
Linear Transforms

Stanislav G. Sedukhin
July 2, 2012

Graduate School of Computer Science and Enginering
The University of Aizu
Tsuruga, Ikki-Machi, Aizu-Wakamatsu City
Fukushima, 965-8580 Japan

Technical Report 2012-001
Title:
Co-design of the Extremely Scalable Algorithms/Architecture for 3-Dimensional Linear
Transforms

Authors:
Stanislav G. Sedukhin
Key Words and Phrases:
3-dimensional linear transforms, multilinear tensor-by-matrix multiplication, computational
index space, cubical data decomposition, modular scheduling, data reusing, torus interconnect,
extremely scalable algorithms, algorithm/architecture co-design
Abstract:
The 3-dimensional (3D) forward/inverse separable discrete transforms, such as Fourier transform, cosine/sinetransform, Hartley transform, and many others, are frequently the principal limiters that prevent many practicalapplications from scaling to the large number of processors. Existing approaches, which are based on 1D or 2D
data decomposition, do not allow the 3D transforms to be scaled to the maximum possible number of computer
nodes. Based on the newly proposed 3D decomposition of an N ×N ×N initial data into P ×P ×P blocks, where
P = N/b and b ∈ [1, 2, ..., N ] is the blocking factor, we systematically design unified, highly scalable algorithms
for parallel implementation of any forward/inverse 3D transform on the one-to-one correspondent P×P×P torus
network of computer nodes. All designed algorithms require 3P “compute-and-roll” time-steps, where each stepis equal to the time of execution in each node b4 fused multiply-add (fma) operations and concurrent movementof O(b3 ) scalar data between nearest-neighbor nodes. The proposed 3D orbital algorithms gracefully avoid arequired 3D data transposition and can be extremely scaled up to the maximum number of N 3 simple nodes(fma-units) which is equal to the size of initial data.

Report Date:
7/2/2012

Written Language:
EnglishAny Other Identifying Information of this Report:
Revised and extended version of the Technical Report 2011-001, the University of Aizu, 2011,
24 p.
Distribution Statement:
First Issue: 10 copies

Supplementary Notes:

Distributed Parallel Processing Laboratory
The University of Aizu
Aizu-Wakamatsu
Fukushima 965-8580
JapanCo-design of Extremely Scalable Algorithms/Architecture
for 3-Dimensional Linear Transforms
Stanislav G. SedukhinDistributed Parallel Processing Lab., The University of Aizu,
Aizuwakamatsu City, Fukushima 965-8580, Japane-mail: sedukhin@u-aizu.ac.jp
July 2, 2012

Abstract
The 3-dimensional (3D) forward/inverse separable discrete transforms, such as Fourier transform,
cosine/sine transform, Hartley transform, and many others, are frequently the principal limiters thatprevent many practical applications from scaling to the large number of processors. Existing approaches,
which are based on 1D or 2D data decomposition, do not allow the 3D transforms to be scaled to the
maximum possible number of computer nodes. Based on the newly proposed 3D decomposition of an
N × N × N initial data into P × P × P blocks, where P = N/b and b ∈ [1, 2, ..., N ] is the blocking
factor, we systematically design unified, highly scalable algorithms for parallel implementation of any
forward/inverse 3D transform on the one-to-one correspondent P × P × P torus network of computer
nodes. All designed algorithms require 3P “compute-and-roll” time-steps, where each step is equal to
the time of execution in each node b4 fused multiply-add (fma) operations and concurrent movement ofO(b3 ) scalar data between nearest-neighbor nodes. The proposed 3D orbital algorithms gracefully avoid
a required 3D data transposition and can be extremely scaled up to the maximum number of N 3 simple
nodes (fma-units) which is equal to the size of initial data.

1

Introduction

Three-dimensional (3D) discrete transforms (DT) such as Fourier transform, cosine/sine transform, Hartleytransform, Walsh-Hadamard transform, etc., are known to play a fundamental role in many application
areas such as spectral analysis, digital filtering, signal and image processing, data compression, medical
diagnostics, etc. Increasing demands for high speed in many real-world applications have stimulated thedevelopment of a number of Fast Transform (FT) algorithms, such as Fast Fourier Transform (FFT), withdramatic reduction of arithmetic complexity [9]. These recursion-based FT-algorithms are deeply serialized
by restricting data reuse almost entirely to take advantage from the sequential processing.
A recent saturation of the performance of single-core processors, due to physical and technological limitations such as memory and power walls, demonstrates that a further sufficient increase of the speed of
FT-algorithms is only possible by porting these algorithms into massively-parallel systems with many-coreprocessors. However, by reason of complex and non-local data dependency between butterfly-connected
operations, the existing deeply serialized FT-algorithms are not well adapted for the massively-parallel implementation. For example, it was recently reported in [34] that, by using two quad-core Intel NehalemCPUs, the direct convolution approach outperforms the FFT-based approach on a 512 × 512 × 512 data
cube by at least five times, even when the associated arithmetic operation count is approximately two times
1higher. This result demonstrates that, because of the higher regularity and locality in the computation anddata access (movement), the convolution achieves significantly better performance than FFT even with ahigher arithmetic complexity. Because cost of arithmetic operations becomes more and more negligible
with respect to the cost of data access or data movement, the conventional matrix-based DT-algorithms witha simple, regular and local data movement are expected to be more suitable for scalable massively-parallel
implementation. As three-dimensional, N ×N ×N discrete transforms are computationally intensive problems with respect to N they are often calculated on large, massively parallel networks of computer nodes,
each of which includes at least one processor that operates on a local memory. The computation of a transform on a network having a distributed memory requires appropriate distribution of the data and work among
the multiple nodes in the network as well as orchestrating of data movement during processing.
For parallel implementation of the 3D discrete transforms (at least for the very popular Fourier transform), there are two distinct approaches which are differ in the way of data decomposition over the physicalnetwork of computer nodes. One approach is the 1D or “slab” decomposition of a 3D N × N × N initial
data which allows scaling (or distribution of work) among np = P = N/b nodes, where a 2D slab of size
N × N × b, is assigned to each node (see Fig. 1 (a)). Here and below, b ∈ {1, 2, ..., N } is the blocking
factor. Different implementations of the 3D FFT with a “slab” decomposition can be found, for example, in[14, 16, 17]. Although this approach has relatively small inter-node communication cost, the scalability ofthe slab-based method or the maximum number of nodes is limited by the number of data elements along a= N when b = 1.
single dimension of the 3D transform, i.e. nmax
p
(a)
N

(b)

(c)

b

b

N
b
b

N

N

b

b
2 ≤ np ≤ N

22 ≤ np ≤ N2

23 ≤ n p ≤ N3

Figure 1: 3D data distribution over np computer nodes: (a) 1D or “slab’ decomposition, (b) 2D or “pencil”
decomposition, and (c) 3D or “cube” decomposition.
Another approach, shown in Fig. 1 (b), is the 2D or “pencil” decomposition of a 3D N ×N ×N initial
data among a 2D array of np = P × P nodes (P = N/b) where a 1D “pencil” or “rod” of size N × b × b,b ∈ {1, 2, ..., N }, is assigned to each node. Parallel 3D FFT implementations with a 2D data decomposition
are discussed in [13, 40, 32]. This approach overcomes the scaling limitation inherent into previous method
since it increases the maximum number of nodes in the system from N to N 2 . As a consequence, each noderequires less memory than in a “slab”-based approach. However, this increasing the number of nodes leads
to rising of the communication cost.
It is important to note that in both of these so-called “transpose” approaches, the computational partand inter-node communication part are separated [6]. Moreover, a computational part on an assigned data
is implemented inside each node of a network by using either 2D or 1D fast (recursive) algorithms for“slab”- or “pencil”-based decomposition, respectively, without any inter-node communication. However,on completion of each computational part, in order to support contiguity of memory access, a transpositionof the 3D data array is required to put data in an appropriate dimension(s) into each node. At least one or2two transpositions would be needed for the 1D or 2D data decomposition approaches, respectively. Eachof this 3D data transposition is implemented by “all-to-all” inter-node, message-passing communication.
This global data exchange imposes an overhead which is proportional to the number of nodes and can bea dominant factor in the total time of processing for even small number of nodes [12]. As it was recently
reported [41], by using a 2D data decomposition approach, the communication part of a 5123 FFT on the
BlueGeene/P computer with up to 512 × 512 = 262, 144 cores (processors) may require more than 95% of
the runtime.
The proposed in this paper a new, transpose-free approach for parallel implementation of the 3-dimensional
discrete transforms improves even further a (strong) scalability of the 3D transforms by increasing the maximum number of computer nodes from N 2 to the extreme number of N 3 . It becomes possible due to
1. the 3D or “cubic” decomposition of an N ×N ×N initial data among np = P × P × P computer
nodes where a 3D data “cube” of size b × b × b is assigned to each node (see Fig. 1 (c));
2. computing of the basic one-dimensional N -size transform not on a single, but on the P = N/bcyclically interconnected (for data reuse) nodes of a 3D torus network by employing not recursive,
but blocked matrix multiplication based algorithms with a well-structured data access/movement at
the expense of increasing the number of arithmetic operations;3. integration of a local, intra-node computation with a nearest-neighbor inter-node communication ateach step of 3-dimensional processing.
A 3D transform is represented as three chained sets of the cubical tensor-by-matrix or matrix-by-tensor
multiplications which are executed in a 3D torus network of computer nodes by the fastest and extremely
scalable orbital algorithms. We call an algorithm as extremely scalable if it can be scaled to the maximum
number of computer nodes which is limited only by the size of initial data, i.e. in our case, N 3 . In otherwords, the number of simultaneously executed scalar (basic for a given algorithm) operations may be equal
to the amount of initial data. Of course, exascale does not mean exaflop.
The paper is organized as follows. In Section 2, the scalar and block notations of the 3D forward andinverse, separable transforms are described. It is shown that these notations are based on the multilinear
matrix multiplication. Section 3 shortly discusses a systematic way for selection the fastest and extremely
scalable matrix-matrix multiplication algorithms and its chaining. Section 4 introduces a 3D block data
partition of the three-way tensor and proposes orbital, highly-scalable algorithms for parallel implementationof the 3D forward and inverse transforms on a 3D network of toroidally interconnected computer nodes. The
paper concludes in Section 5.

2

3D Separable Transforms

Let X = [x(n1 , n2 , n3 )], 0 ≤ n1 , n2 , n3 < N , be an N ×N × N cubical grid of input data or three-way
data tensor [24]. ...
A separable forward 3D transform of X is another cubical grid of an N ×N ×N data or...
three-way tensor X = [ x(k1 , k2 , k3 )] where for all 0 ≤ k1 , k2 , k3 < N :...
x(k1 , k2 , k3 ) =

N
−1 N
−1 N
−1
X
X
X

x(n1 , n2 , n3 ) · c(n1 , k1 ) · c(n2 , k2 ) · c(n3 , k3 )n3 =0 n2 =0 n1 =0

(1)

...
...
In turn, a separable inverse or backward 3D transform of three-way tensor X = [ x(k1 , k2 , k3 )] is expressed
as:
N
−1 N
−1 N
−1
X
X
X
...
x(k1 , k2 , k3 )c(n1 , k1 ) · c(n2 , k2 ) · c(n3 , k3 )x(n1 , n2 , n3 ) =(2)
k3 =0 k2 =0 k1 =03

where 0 ≤ n1 , n2 , n3 < N and X = [x(n1 , n2 , n3 )] is an output N ×N ×N cubical tensor.
There is a direct correspondence between the equations (1), (2) and the so-called
symmetric multilinear...
matrix multiplication, which is used to represents three-way data tensor X or X in different bases and where
an N × N matrix C = [c(ns , ks )] = [c(n, k)], s = 1, 2, 3 is a (non-singular) change-of-basis matrix(see [24, 33, 22] for more details). It is also interesting to mention that equations (1) and (2) can be viewed
as the so called Tuker’s 3D tensor decomposition which is represented in the form of three-way tensor-bymatrix multiplication [24, 30]. Moreover, the equation (1) or (2) represents the so-called three-way tensor
contraction which is widely used in ab initio quantum chemistry models [5, 35, 47].
The separable transforms differ only by the transform coefficient (change-of-basis) matrix C = [c(n, k)]which can be• symmetric, i.e. C = C T , and unitary, i.e. C −1 = C ∗T , C ∗ is a complex conjugate of C, like in the
2πnk
2πnk
Discrete Fourier Transform (DFT), where c(n, k) = exp [− 2πi
N (n · k)] = cos( N ) − i sin( N )√2πnk
and i = −1, or in the Discrete Hartley Transform (DHT), where c(n, k) = cos( 2πnk
N ) − sin( N );• unitary and real, i.e. orthogonal, like in the Discrete Cosine Transform (DCT), where coefficientπc(n, k) = cos [ 2N
(2n + 1) · k] and C 6= C T ;• consists only ±1 and be symmetric and orthogonal, like in the Discrete Walsh-Hadamard Transform(DWHT).
We will abbreviate the generic form of a discrete transform as DXT without taking into account the specificfeatures of a coefficient matrix, i.e. we will use direct algorithms for the 3D transforms (1) and (2) with anarithmetic complexity of O(N 4 ) instead of using the so-called “fast” DXT algorithms, like 3D FFT, withO(N 3 log N ) complexity.
n3

n3
N

N

N

N
n2
(a)

n1

P

P

N

n1

N3
P

P
b

b

N2 N1
(b)

b

N2

(c)

Figure 2: Partitioning of an N × N × N initial data (a) into P slabs along N2 -axis with the size of an
N × b × N each (b) and P × P × P cubes with the size of b × b × b each (c).

The first step to design a scalable 3D DXT is to represent it in a block matrix notation where theconventional scalar form (1) or (2) is the only one extreme case of a block notation. To formulate ourproblem in this block notation we firstly divide an N × N × N input data volume or three-way tensor
X = [x(n1 , n2 , n3 )] into P × P × P data cubes, where each cube X(N1 , N2 , N3 ), 0 ≤ N1 , N2 , N3 < P,
has the size of b × b × b, i.e. b = N/P and 1 ≤ b ≤ N/2 is the blocking factor. Then the forward 3D DXT

4

(3D FDXT) can be expressed as a block version of the multilinear matrix multiplication:
P
−1 P
−1 P
−1
X
X
X
...
X(N1 , N2 , N3 ) × C(N1 , K1 )×C(N2 , K2 )×C(N3 , K3 ),X(K1 , K2 , K3 ) =(3)

N3 =0 N2 =0 N1 =0

where 0 ≤ K1 , K2 , K3 < N and C(Ns , Ks ) = C(N, K), s = 1, 2, 3, is an (Ns , Ks )-th block of thetransform matrix C.
It is clear that a 3D block inverse transform (3D IDXT) can be written as:
X(N1 , N2 , N3 ) =

P
−1 P
−1 P
−1
X
X
X

...
X(K1 , K2 , K3 ) × C(N1 , K1 )×C(N2 , K2 )×C(N3 , K3 ),(4)

K3 =0 K2 =0 K1 =0

where 0 ≤ N1 , N2 , N3 < P . As it can be seen from (3) and (4), both 3D transforms are implemented ina 6D computational index space I6D = {(N1 , N2 , N3 , K1 , K2 , K3 ) : 0 ≤ N1 , N2 , N3 , K1 , K2 , K3 < P }with the number of index points kI6D k = P 6 . Each index point in this 6D space is associated with three
cubical tensor-by-matrix multiplications.
Due to separability of the linear transforms, the dimension of a computational index space can be reduced from six to four by splitting a 3D transform into three data dependent sets of 1D transforms as it isshown below for the 3D FDXT (3).
At the first stage, the P × P 1D FDXT of X(N1 , N2 , :) are performed for all (N1 , N2 ) pairs, 0 ≤N1 , N2 < P, as block cubical tensor-by-matrix multiplication:Ẋ(N1 , N2 , K3 ) =P−1
X

X(N1 , N2 , N3 ) ×C(N3 , K3 ),

(5)

N3 =0

where 0 ≤ N1 , N2 , K3 < N . It is clear from (5) that this stage is implemented in a 4D computational index
space
I
I4D
= {(N1 , N2 , N3 , K3 ) : 0 ≤ N1 , N2 , N3 , K3 < P }with an embedded 3D input data tensor {X(N1 , N2 , N3 ), 0 ≤ N1 , N2 , N3 < P } as one out of eight 3D faces
I . This cubical tensor is shown in red color in Fig. 3(a). For simplicity, a graphicalof 4D index space I4D
example here is provided for P = 2, i.e. blocking factor b = N/2.
At the second stage, the P × P 1D FDXT of Ẋ(:, N2 , K3 ) are implemented for all (N2 , K3 ) pairs,
0 ≤ N2 , K3 < P, as second block tensor-by-matrix multiplication:
Ẍ(K1 , N2 , K3 ) =P−1X

Ẋ(N1 , N2 , K3 ) ×C(N1 , K1 ),

(6)

N1 =0

where 0 ≤ K1 , N2 , K3 < N . This stage is also implemented in a 4D index spaceII
I4D= {(N1 , N2 , K3 , K1 ) : 0 ≤ N1 , N2 , K3 , K1 < P }T II
I/II = I I
with a common to the previous stage 3D face I3D
I4D = {(N1 , N2 , K3 ) : 0 ≤ N1 , N2 , K3 < P }4Dwhich is shown in yellow in Fig. 3(b).
At the third stage, the P × P 1D FDXT of Ẍ(K1 , :, K3 ) are implemented for all (K1 , K3 ) pairs, 0 ≤K1 , K3 < P, as third block tensor-by-matrix multiplication:
P
−1
X
...
Ẍ(K1 , N2 , K3 ) ×C(N2 , K2 ),X(K1 , K2 , K3 ) =
N2 =0

5

(7)

(a)

(b)
K1

N3

N3

N2

N2
N1
N1
K3

K3
K1

(c)

K1

(d)

N2

N2

K2

K2

N1

K3

K3

Figure 3: United 4D (mesh-based) computational index space of the 3D DXT: (a) index space for the firststage with an initial cubical data grid shown in red; concatenation of the 4D index spaces for the (b) firstand second and (c) second and third stages with a common intermediate cubical data shown in yellow andgreen, respectively; (d) the final index space with a computed data cube shown in blue.

6

where 0 ≤ K1 , K2 , K3 < N . A 4D computational index space for this final stage isIII
I4D= {(K1 , N2 , K3 , K2 ) : 0 ≤ K1 , N2 , K3 , K2 < P }.
T III
II/III = I IIThis index space includes a common to the previous stage 3D face I3D
I4D
4D... = {(K1 , N2 , K3 ) :0 ≤ K1 , N2 , K3 < P } (shown in green in Fig. 3(c)) and a 3D output data tensor {X(K1 , K2 , K3 ), 0 ≤K1 , K2 , K3 < P } which is depicted in blue color in Fig. 3(d).
It is clear that, in general, these three stages can be implemented in any of six possible orders of summation and the total number of block cubical tensor-by-matrix multiplications is reduced from P 6 , as for the
direct computing by (3) or (4), to 3P 4 which is equal to the number of index points in a combined 4D index
space
[
[
I
II
III
I4D = I4D
I4D
I4D
.
(8)

By slicing cubical data, i.e. representing three-way tensors as the set of matrices, it is possible toformulate a 3D transform (3) or (4) as conventional block matrix-by-matrix multiplication. In this case, an
initial P×P×P data grid or three-way tensor {X(N1 , N2 , N3 ), 0 ≤ N1 , N2 , N3 < P }, is divided into P 1D
“slices” or “slabs” or “matrices-of-cubes” along one of axises, for example, along N2 -axis, such that eachb × b × b data cube X(N1 , N2 , N3 ) can be referred as a block element X(N1 , N3 )N2 of the N2 -th P × P
matrix, where N2 ∈ [0, P ) (see Fig. 2). Then a 3D FDXT (3) can also be computed in three data-dependentstages as chaining sets of the block matrix-by-matrix products with a 3D computational index space for eachproduct.
(b)

(a)

N3

(c)
K1

K1
N2

N2

N2

K2
N1
N1
K3

K3

K3

Figure 4: Partition of the mesh-based 4D computational index space into the set of 3D disjoint index spaces(opposite parallel cubical sides) for the first (a), second (b), and third (c) stages of computing.

At the first stage, since there is no dependency in 1D transforming data between all N2 -slabs, thecubical tensor-by-matrix multiplication (5) can be presented as P, independent on N2 , block matrix-bymatrix multiplications:Ẋ(N1 , K3 )N2 =P−1
X

X(N1 , N3 )N2 × C(N3 , K3 )

N3 =0

7

(9)

where the same change-of-basis matrix C = [C(N3 , K3 )], 0 ≤ N3 , K3 < P, is used for all N2 -slices,N2 ∈ [0, P ). As it follows from (9), each block matrix-by-matrix product in a N2 -th slab is implemented in
its own 3D index space with P 3 index poins
I
I3D(N2 ) = {(N1 , N3 , K3 ) : 0 ≤ N1 , N3 , K3 < P }.

It is clear that all these cubical index spaces are disjoint, i.e.
\I
I3D
(N2 ) = Ø,0≤N2 <P

and, therefore, an associated with each 3D index space set of computing can be implemented independently
or in parallel.
I (0)For the smallest 2×2 matrix case, i.e. when b = N/2 and, therefore, P = 2, the 3D index spaces I3DI
and I3D (1) can be seen in Fig. 4(a) as opposing parallel sides (cubes in bold lines) along the N2 -axis. In
I (N ) is associated with a conventional
these 3D grids of index points, each point p = (N1 , N3 , K3 )T ∈ I3D
2
block matrix-by-matrix multiply-add
Ẋ(N1 , K3 )N2 ← X(N1 , N3 )N2 × C(N3 , K3 ) + Ẋ(N1 , K3 )N2 ,(10)

where an addition or accumulation is performed along N3 -direction and a block matrix Ẋ(N1 , K3 )N2 is a
cubical b × b × b data tensor which should be zeroed initially. Now, for the original first stage (5), the 4Dcomputational index space can be defined as
[
I
I
I4D
=
I3D
(N2 ).
0≤N2 <P

A 3D index space for the computed on this stage intermediate cubical tensor Ẋ(N1 , N2 , K3 ) can be seen in
Fig. 3(b) in yellow color.
At the second stage, the original tensor-by-matrix product (6) is computed by the set of P, also independent on N2 , block matrix-by-matrix multiplications in the following form which keeps the required
row-by-column index agreement:
Ẍ(K1 , K3 )N2 =P−1X

C(N1 , K1 )T × Ẋ(N1 , K3 )N2 .

(11)

N1 =0

It can be seen from (11) that the same coefficient matrix C = [C(N1 , K1 )], 0 ≤ N1 , K1 < P, is used
for all N2 -slices, N2 ∈ [0, P ). As in the previous stage, computing for each N2 -slab is implemented in a 3D
index spaceII
I3D(N2 ) = {(N1 , K1 , K3 ) : 0 ≤ N1 , K1 , K3 < P }.
The same as above, all these 3D index spaces are disjoint, i.e.
\II
I3D
(N2 ) = Ø,0≤N2 <P

and, therefore, the related to each 3D index space set of computing can also be implemented independently.
II (0) and I II (1) are shown for P = 2 in Fig. 4(b) as two opposing parallel sides
The 3D index spaces I3D
3D
II (N ) is
(cubes) along the same as in the prior stage N2 -axis. Each index point p = (N1 , K1 , K3 )T ∈ I3D2
associated here with a block multiply-add computingẌ(K1 , K3 )N2 ← C(N1 , K1 )T × Ẋ(N1 , K3 )N2 + Ẍ(K1 , K3 )N2 .
8(12)

In this stage, an addition is performed along N1 -direction and a block matrix Ẍ(K1 , K3 )N2 is a b × b × b
data tensor which should be zeroed initially. For the original stage (6), the 4D computational index space isdefined as
[
II
II
I4D
=
I3D
(N2 ).
0≤N2 <P

A 3D index space for the computed on this stage intermediate cubical tensor Ẍ(K1 , N2 , K3 ) can be seen in
Fig. 3(c) in a green color.
At the third stage, the final tensor-by-matrix multiplication (7) should be implemented inside the 4D
index space by summation data along N2 -direction. This can be done by 1D slicing or partition of thecomputed at the previous stage (see Fig. 3(c)), cubical tensor Ẍ(K1 , N2 , K3 ), 0 ≤ K1 , N2 , K3 < P, along
K3 -axis, such that each K3 -th slab will keep the all nedeed data at the index points along N2 -direction, i.e.
[
Ẍ(K1 , N2 , K3 ) =Ẍ(K1 , N2 )K3 .
0≤K3 <P

With this data partition, the tensor-by-matrix product (7) is implemented as the set of P, independent on K3 ,block matrix-by-matrix multiplications:
P
−1
X
...
(K
,
K
)
=
Ẍ(K1 , N2 )K3 × C(N2 , K2 ),X 1 2 K3

(13)

N2 =0

where the same P × P matrix C = [C(N2 , K2 )], 0 ≤ N2 , K2 < P, is used for all K3 -slices, K3 ∈ [0, P ).
A computing in each K3 -slab is implemented in its own 3D index spaceIII
I3D(K3 ) = {(K1 , N2 , K2 ) : 0 ≤ K1 , N2 , K2 < P }.

Again, all these 3D index spaces are disjoint, i.e.
\

III
I3D
(K3 ) = Ø,0≤K3 <P

and, therefore, the associated with each 3D index space set of computations can also be implemented independently.
III (0) and I III (1) are shown in Fig. 4(c) as opposing
For the case of P = 2, the 3D index spaces I3D
3D
III (K ) is associated withparallel sides along the K3 -axis. Here, each index point p = (K1 , N2 , K2 )T ∈ I3D3a block matrix-by-matrix multiply-add
...
...
(14)
X(K1 , K2 )K3 ← Ẍ(K1 , N2 )K3 × C(N2 , K2 ) + X(K1 , K2 )K3 .
...
In this final stage, an addition is performed along N2 -direction and a block matrix X(K1 , K2 )K3 is a b×b×b
data tensor which should be zeroed initially. It is clear that for the original stage (7), the 4D index space isdefined as
[
III
III
I4D
=
I3D
(K3 ).
0≤K3 <P...
A 3D index space for the resulting cubical tensor X(K1 , K2 , K3 ) is shown in Fig. 3(d) in a blue color .
Note that in (9), (11) and (13) each of P products is three-way b × b × b tensor by b × b matrix multiplication [4]. Moreover, for each N2 -slice of a cubical data X, the first two stages perform independently
the canonical 2D FDXT in the form of triple matrix multiplication with one matrix transposed, i.e., in the9matrix form, Ẍ = C T ×X×C. It is clear that in the first stage (9), N3 -axis from the original basis is converted into K3 -axis of the new basis, i.e. transforming [N3 → K3 ] is implemented. In the second stage(11), changing [N1 → K1 ] is performed (see Fig. 5 (a)). Before the third stage, where the final changing [N2 → K2 ] should be executed, an 1D slicing or working partition of the intermediate cubical tensor
Ẍ along N2 -axis is algorithmically changed to the required by (13) N2 -summation by using an 1D partition along K3 -axis (see Fig. 5 (b)). This changing of partition and summation direction is possible becauseẌ(K1 , K3 )N2 ≡ Ẍ(K1 , N2 )K3 , i.e. the same b×b×b data cube simultaneously belongs to N2 - and K3 -slices
or slabs. Actually, a slab’s sub-index in the above equations is used only to represent a 3D DXT computing
in the form of canonical (planar) matrix-by-matrix multiplication where the row-by-column summation index agreement is required. It is clear that this slab’s sub-index can be included in the main 3-dimensional
index when direct multilinear multiplication of a three-way tensor by a 2D matrix is utilized [33, 39].
N3

N1

K3

K1

K3

N2

N2

K1

(a)

K2

(b)

Figure 5: Partition of a 3D data tensor into slabs (slices) of cubes (a) along N2 -axis for the first and second
stages and (b) along K3 -axis for the third stage.
It is important to mention that an initial slicing or 1D partition of a given cubical tensor affects the orderin which the data dependent sets of 1D transforms are implemented. Excluding an existing in our cubicaltensor (super)symmetry [24], the selected above initial data partition along N2 -axis or the so-called lateralslicing [24], leads to the following 3-stage transform order:[N3 → K3 ] ⇒ [N1 → K1 ] ⇒ [N2 → K2 ].
A partition of the same cubical tensor along N1 -axis or frontal data slicing, will result the order:[N3 → K3 ] ⇒ [N2 → K2 ] ⇒ [N1 → K1 ].
Finally, a cubical data partition along N3 -axis or horizontal slicing, will demand the order:[N2 → K2 ] ⇒ [N1 → K1 ] ⇒ [N3 → K3 ].
It is clear that an inverse 3D transform (IDXT) (4) is implemented in the reverse order,...i.e. as rollingback of a forward 3D DXT. By keeping the slicing of an initial cubical P × P × P tensor X(K1 , K2 , K3 )along K3 -axis, a 3D IDXT would require implementation of the following three stages:
Stage I: for all pairs (K1 , N2 ) at slabs K3 ∈ [0, P ) doẌ(K1 , N2 )K3 =

P
−1
X

...
T
X(K1 , K2 )K3 × C(N2 , K2 ) , 0 ≤ K1 , N2 < P.

K2 =0

10

This first stage for a 3D IDXT is computed in the same as for the third stage of a 3D FDXT (7) 4D index
space
III
I4D
= {(K1 , K2 , K3 , N2 ) : 0 ≤ K1 , K2 , K3 , N2 < P }.
After completion, a resulting in this stage cubical tensor Ẍ(K1 , N2 , K3 ) is used as sliced into 1D slabs along
N2 -axis.
Stage II: for all pairs (N1 , K3 ) at slabs N2 ∈ [0, P ) do
Ẋ(N1 , K3 )N2 =P−1X

C(N1 , K1 ) × Ẍ(K1 , K3 )N2 , 0 ≤ N1 , K3 < P.

K1 =0

The same as for the equation (6), this stage is implemented in a 4D index spaceII
I4D= {(K1 , N2 , K3 , N1 ) : 0 ≤ K1 , N2 , K3 , N1 < P }T II
III/II = I III
with a common 3D face I3D
I4D = {(K1 , N2 , K3 ) : 0 ≤ K1 , N2 , K3 < P }.
4D

Stage III: for all pairs (N1 , N3 ) at slabs N2 ∈ [0, P ) do
X(N1 , N3 )N2 =P−1X

Ẋ(N1 , K3 )N2 × C(N3 , K3 )T , 0 ≤ N1 , N3 < P.

K3 =0

The final third stage of the 3D IDXT is implemented in the same as for the first stage of a 3D FDXT (5) 4D
index spaceI
I4D= {(N1 , N2 , K3 , N3 ) : 0 ≤ N1 , N2 , K3 , N3 < P }T I
II/I = I II
with a common 3D face I3D
I4D = {(N1 , N2 , K3 ) : 0 ≤ N1 , N2 , K3 < P }.
4D
It is easy to verify that the total number of P × P block matrix-by-matrix multiplications in a forward
or inverse 3D DXT is 3P 4 . Each block matrix-by-matrix multiplication is a b × b × b tensor-by-matrix
product which require an execution of b4 scalar fused (indivisible) multiply-add (fma) operations, whereb = N/P is a blocking factor. The total number of such scalar fma-operations for a 3D DXT is, therefore,
3P 4 · b4 = 3N 4 , i.e. blocking does not change an arithmetic complexity of the transformation. Obviously,
this is because our 3D DXT implementation is totally based on a matrix-by-matrix multiplication. At eachstage of computing, the maximal degree of reusing of the N 3 data elements or an arithmetic density can be
estimated as N 4 /N 3 = N scalar multiply-add operations per data element. Note that the existing so-called“fast” algorithms for the 3D DXT, like 3D FFT, require an execution of O(N 3 log N ) scalar arithmetic
operations, i.e. an arithmetic density is only O(log N ) operations per data element. This low degree of datareusing is one of the main reasons of low efficiency and poor scalability of parallel implementations of thefast DXT algorithms.

3

Exa-Scalable Matrix-by-Matrix Multiplication and 3D Transforms

3.1

Fast Matrix-by-Matrix Multiplication

It is clear from the above discussion that a 3D FDXT of a cubical data array or three-way tensor X can be
expressed as multilinear matrix multiplication
3D FDXT: N2 →K2z}|z }| {X×C{z{

1D FDXT: N3 →K3...
T
X= C ×|

2D FDXT: N1 →K1

11


}

×C,

(15)

and a 3D IDXT as

3D IDXT: K3 →N3z}|z... }| {
T
X×C
{z

{

1D IDXT: K2 →N2


X= C×
|

2D IDXT: K1 →N1



×C T .

(16)}

Recall that this fixed order of parentetization directly corresponds to 1D initial partition (slicing) of the input
three-way tensor along N2 -axis which has been selected above as one of three possible 1D partitions (seeFig. 2 (b)).
From the equations (15) and (16) it is evident that to implement a 3D transform under unlimited or
extreme parallelism, i.e. when the maximum number of simultaneous operations is limited only by the sizeof data, the answers to the following questions should be found:• What is (are) the fastest and technologically justified parallel matrix-by-matrix multiply-add (MMA)
algorithm(s) with the maximal data reuse and extreme scalability?
• How to design a collection of parallel MMA algorithms with different transposed/non-transposed
scenarios, like in a famous GEneral Matrix-matrix Multiplication (GEMM) from the Level-3 BLAS
[11]?
• How the different sets of MMA operations can be chained in time-space to effectively implement fastmultilinear matrix multiplication (15) and (16) without any data transposition?
• How it can S
be possible
S IIIto map three, concatenated by a common 3D data, 4D computational index
IIIspaces I4D I4D I4D of a 3D DXT (see Figures 3 and 4) into a 3-dimensional physical (spatial)
space of execution or processor space?
Some answers to these questions can be found in our previous paper [46] where extremely scalable scalar
MMA algorithms for the two-dimensional DXT on a 2D torus array processor have been systematicallydesigned and evaluated. The main results related to the systematic design of the scalable MMA algorithms
are described below. In [46] , it was formally shown that for the following three forms of the matrix-bymatrix multiply-add:
C ← A × B + C, where c(i, j) ←n−1
Xa(i, k) · b(k, j) + c(i, j);

k=0
n−1
X

A ← C × B T + A, where a(i, k) ←(17a)c(i, j) · b(k, j) + a(i, k);(17b)a(i, k) · c(i, j) + b(k, j);

(17c)

j=0

B ← AT × C + B, where b(k, j) ←

n−1
X
i=0

where A = [a(i, k)], B = [b(k, j), and C = [c(i, j)] are matrices of the same n × n size, there exist threeclasses of the fastest MMA algorithms which are efficiently implemented on a planar n×n mesh/torus arrayprocessor in n time-steps. These three classes are differ in linear arrangements (time-step scheduling) of
the partially-ordered set of n3 scalar indivisible (fused) multiply-add operations which are algorithmically“located” in a cubical n × n × n computational index space I3D = {(i, j, k)T : 0 ≤ i, j, k < n}. Actually,
every scheduling class defines for each step(p) ∈ {0, 1, 2, ..., n − 1} the set of n2 “active” index points
p = (i, j, k)T ∈ I3D with an associated data-triplet {a(i, k), b(k, j), c(i, j)}, where partial computing of
12

(17) is implemented by concurrent execution of n2 scalar multiply-add operations from the correspondingset:c(i, j) ← a(i, k) · b(k, j) + c(i, j);(18a)a(i, k) ← c(i, j) · b(k, j) + a(i, k);(18b)b(k, j) ← a(i, k) · c(i, j) + b(k, j).

(18c)

The scheduled at each time-step set of n2 “active” index points directly defines an associated data distribution in 3D index space including an initial data placement, which is equal to the data allocation at thefirst time-step. Additionally, this set of index points also defines, but indirectly, a specific for each class ofscheduling way of data reuse. As a result of indirect definition of data reuse, the different scenarios of its
physical implementation may exist as it will be shown below.
The time-step scheduling function is used in the linear or modular form:
step(p) = (αT · p) or step(p) = (αT · p) mod n,
respectively, where α = (αi , αj , αk )T is a scheduling vector and (·) is a scalar product of two vectors.
Based on the different space-time data arrangements in a 3D index space I3D , the following classificationof MMA algorithms can be established (see [46] for details).
• Broadcast-Broadcast-Compute (BBC) class:– the rank-one [48] or outer-product [19] update is used here where computational index space
I3D is represented as a 3D mesh of points;– there is only one linear scheduling function for this class: step(p) = k, i.e. α = (0, 0, 1)T ;– linear projection of the 3D index space along summation axis into a 2D physical space gives awell-known broadcast- or replication-based MMA implementation on a planar array processor[2, 48];• Broadcast-Compute-Roll (BCR) class:– a computational index space I3D is a 3D cylindrically interconnected points;– there are two (orthogonal) modular scheduling functions: step(p) = (k − i) mod n andstep(p) = (k − j) mod n, i.e. α = (−1, 0, 1)T or α = (0, −1, 1)T ;– a 3D → 2D projection of the cylindrical index space also along summation axis into a 2D
physical space gives the well-known Fox’s algorithm on a non-planar, semi-toroidal, array processor [15];• Compute-Roll-All (CRA) or orbital class:– an index space I3D is a 3D toroidally interconnected points;– excluding a symmetry, there are four distinct modular scheduling functions:step(p) = (±i ±j ±k) mod n, i.e. α = (±1, ±1, ±1)T ;– different admissible linear 3D → 2D projections of the torus index space into a 2D physicalspace produce a variety of Cannon-like algorithms [7, 46] on a 2D torus array processor.

13

Recall that all three MMA forms (17) are equally scheduled, i.e. if transposition is required, a corresponding
operation (17b) or (17c) is implemented with the same initial data distribution in a 3D index space and withthe same data movement pattern as for (17a), i.e. actual matrix transposition is avoided.
Different linear or modular arrangements (schedules) of n3 data dependent scalar operations impose a
different style in reusing of the matrices A and B for concurrent updating n2 elements of a matrix C at eachtime-step 1 . That is, at each of n iterations or consecutive time-steps of computing:• an algorithm from the BBC class reuses (by data replication or copying, since αi = 0 and αj = 0)n data elements (one vector-column) of a matrix A and n data elements (one vector-row) of a matrix
B, i.e. totally 2n elements are reused for updating n2 elements of a matrix C which are then shifted(since αk = 1) within a cubical index space I3D to the next iteration step;• algorithms from the BCR class reuse n data elements (vector-diagonal) of a matrix A or B (by replication, since αj = 0 or αi = 0) and all n2 elements of a matrix B or A, respectively, i.e. n2 + nelements are reused for updating n2 elements of a matrix C and, after updating, the all 2n2 elementsof matrices C and B or A, respectively, are circularly shifted or rolled (since αk = 1 and αi = −1or αj = −1) inside cylindrical index space I3D to the next step along or opposite corresponding
direction (orbit) according to the sign;• algorithms from the CRA class reuse all elements of matrices A and B, i.e. totally 2n2 data elements
are reused here for updating n2 elements of a matrix C and, after updating, the all 3n2 elements of thematrices A, B and C are circularly shifted (since αi = ±1, αj = ±1, αk = ±1) inside torus index
space I3D to the next step also along or opposite corresponding direction (orbit) according to the sign.
Note that the well-known planar systolic MMA algorithms [26, 28, 27], which are belong to the so-called
All-Shift-Compute (ASC) class, require 3n − 2 time-steps on a mesh array processor and, therefore, they
are not as fast as other discussed here n-steps algorithms. Moreover, in the systolic array processors, theinitial data streams should be appropriately skewed and, in some cases, scattered before main processing isstarted. An additional, non-trivial parallel memory, surrounding a systolic array processor, is usually needed
to keep and form these initial/final data streams. These challenging requirements sufficiently complicate apractical implementation of systolic array processors. Note that ASC class is a localized (pipelined) versionof the BBC class where the time-step function is defined as step(p) = i + j + k, i.e. α = (1, 1, 1)T .
The provided classification makes it clear that MMA algorithms from BBC class are the most economical in terms of the quantity of reusing data but, at the same time, this class require a global, i.e. depending on
the size of the problem, data replication among n2 different operations. This global data replication can bephysically implemented by either data broadcast (“one-to-all”) or by its antipode, concurrent access (“allto-one”) to this reusable data in a shared memory. If MMA computing is implemented on a synchronous
system where the signal propagation time between different components is no longer than the clock period
then signal/data replication by broadcast or by shared access can be implemented, in principle, within asingle cycle. In other words, MMA algorithms from the BBC class would be effective in any synchronous
system which size is limited and adjusted for the single clock period. This synchronous system acts as asingle node. Because technically, concurrent access to the shared memory is much faster than data replication by broadcast within the same limited size of system, a hierarchical shared memory approach is widelyused in today high-performance CPU and GPU processors. However, it is clear from the existing physicalconstrains, such as finiteness of the speed of light, that a global data broadcast or its alternative, global data
sharing, will prevent a computer system to be extremely scaled to any relatively big size (far beyond the1

An explanation here is related to the canonical form (17a), but form (17b) or (17c) can be equally considered.

14

synchronization area or volume) when the number of concurrent operations is equal to the size of problem,
i.e. the size of matrices.
From the above classification it is also follows that the only Compute-Roll-All schedule defines at eachtime-step a global circulation by local movement of all 3n2 matrix elements without any redundant data
replication. Therefore, only MMA algorithms from the CRA class, which are supported by appropriate
asynchronous local data exchange, can be extremely scaled on the system without a global clock. Note that
due to cyclical (modular) nature of the orbital CRA-based MMA algorithms, the all rolled at each timestep data elements are finally (after n “compute-and-roll” time-steps) returned to the originating (initial)
positions in a 3D torus index space. This unique for the orbital processing property is effectively used forchaining of the different matrix products (see [46] for details). It is clear that BCR class of parallel MMA
algorithms is a combination of n-data replication and n2 -data circulation, i.e. mixing of the BBC and CRA
classes.
Note that to effectively realize the set of highly-repetitive in MMA algorithms scalar operations (18), thefollowing forms of fused multiply-add (fma) instructions should be considered: form 123 does c ← a · b + c,
form 321 does a ← c · b + a and form 132 does b ← a · c + b. The ordering number is just the order of theoperands on the right side of the expression.

3.2

Exa-Scalable Implementation of the 3D Transforms in a 4D Computational SpaceDue to extreme scalability, a Compute-and-Roll (CRA) approach will be used for the orbital implementationof a block multilinear matrix multiplication (15) and (16). As it was previously shown (see Fig. 4), a 4D
computational index space I4D of the 3D DXT from (8) can be represented in each of three data dependentstages as disjoint sets of the 3D index spaces. Each 3D index space represents an independent matrix-bymatrix multiplication. All these disjoint index spaces should be properly, i.e. conflict-free, planned forexecution by applying one of a few possible for CRA class modular scheduling functions which will define
an orbital way of processing. A scheduled for execution 4D torus computational space, which is composedof many 3D tori of index points, should be finally mapped into a 3D space of physical implementation byusing a corresponding 4D → 3D linear projection for each consecutive stage of processing.
At the first stage of a 3D DXT computing (9), the index points p = (N1 , N2 , K3 )T ∈ I3D from eachN2 -th 3D computational index spaceI
I3D(N2 ) = {(N1 , N3 , K3 ) : 0 ≤ N1 , N3 , K3 < P }, N2 ∈ [0, P ),are scheduled for execution and data transferring by using one of the admissible modular functions. This
I (N ) into 3D torus or
modular (circular) scheduling converts each cubical, mesh-based index space I3D2
orbital index space. The same scheduling function is used as for all P disjoint 3D index spaces as well as
for all stages of processing. As an example, the modular scheduling function step(p) = (αT · p) mod P
with the scheduling vector α = (1, 1, 1)T is used here for the smallest case when P = N/b = 2. For
the first stage (see equation (9)) of orbital computing, the result of this space-time scheduling is shown
schematically in Fig. 6 for each of two disjoint 3D orbital index spaces by indicating “active” index points,
which are filled in black, and directions of data movement in two different time-steps (a) and (b). At every
time-step, each “active” index point is associated with a matrix-by-matrix multiply-add which is definedby (10). After computing part, blocks of the matrix data X(N1 , N3 )N2 , C(N3 , K3 ) and Ẋ(N1 , K3 )N2 arerolled along K3 -, N1 -, and N3 -orbits, respectively. These positive directions of data rolling are defined by
the given scheduling vector α = (1, 1, 1)T . Recall that in the first stage, an accumulation of the matrix-bymatrix products is implemented along N3 -axis or orbit.
I (N ),
It is clear that at each time-step only P 2 index points are active in each 3D torus index space I3D23
and, therefore, totally, P index points are active in all index spaces. As it can be seen later, this condition15

holds for all three stages of orbital processing. Therefore, for every stage of the 3D DXT processing, nomore than P 3 computer nodes would be needed at each of P time-steps to implement all P 4 operations
I/II/III . For the first stage, the P unconnected (disjoint) orbital index
associated with all index points in each I4D
I (N ) is scheduled such that “compute-and-roll” processing in each N -th 3D torus index space
spaces I3D2
2
begins at step(p) = N2 , i.e. starting time is skewed for each index space. As it can be seen from Fig. 6
I (0) is scheduled to start from step(p) = 0 while an index spaceon time-step (a), a cubical index space I3DI
I3D (1) is planned to start from step(p) = 1. This initial time skewing in the orbital scheduling guaranties
that all index points with the same coordinates in different 3D orbital spaces will have different time-steps
during all processing. After P “compute-and-roll” time-steps, i.e. after finishing of the first stage, blocks ofintermediate matrix data Ẋ(N1 , K3 )N2 will be located in a 4D index space as it is required by the second
stage. Actually, due to the cyclical nature of data movement, these updated blocks are finally, i.e. after P
time-steps, returned to the original positions (indexes).
The second stage of orbital computing (11) is performed by following the same as for the first stagescheduling policy for each of P disjoint by N2 -axis 3D torus index spaces
II
I3D(N2 ) = {(N1 , K1 , K3 ) : 0 ≤ N1 , K1 , K3 < P }, N2 ∈ [0, P ).

A time-space scheduling of “active” index points is shown in Fig. 6 for time-steps (c) and (d). It can be seenthat there is no difference between two stages as in positioning of the independent 3D torus index spaces
in a 4D space as well as in scheduling of the index points with the same coordinates, despite the fact thatN3 -axis was replaced after the first stage by K1 -axis. Here, at each consecutive time-step, every “active”
index point is associated with a given by (12) matrix-by-matrix multiplication. For this stage, according
to the scheduling vector α = (1, 1, 1)T , a corresponding initial distribution and cyclical movements of anintermediate matrix Ẍ(K1 , K3 )N2 along N1 -axis and a coefficient matrix C(N1 , K1 )T along K3 -axis are
required. A previously computed matrix Ẋ(N1 , K3 )N2 is rolled at this stage along the “same” K1 -axis. An
updating of the intermediate matrix Ẍ(K1 , K3 )N2 is implemented by accumulation of the matrix-by-matrix
products along N1 -axis. It is clear that P additional “compute-and-roll” time-steps are needed to finish
second stage of orbital processing in a 4D index space. On completion of this stage, N1 -axis is replaced by
K2 -axis.
For the third stage of a 3D DXT computing, we keep the same scheduling function for all index points
III as for two previous stages. However, a required by (13) accumulation of products
in a 4D index space I4D
along N2 -axis has to be implemented in the differently oriented 3D torus index spaces
III
I3D(K3 ) = {(K1 , N2 , K2 ) : 0 ≤ K1 , N2 , K2 < P },which are disjoint now along K3 -axis, K3 ∈ [0, P ). The “active” index points and directions of data
movement are shown schematically in Fig. 6 for the final time-steps (e) and (f ) . Each “active” by the given
scheduling index point from a 3D torus index space is associated here with a matrix-by-matrix multiplication(14). Besides, a selected at the beginning of processing modular scheduling ...function defines not only an
initial data distribution, but also directions of moving of a resulting matrix X(K1 , K2 )K3 along N2 -axis
and a coefficient matrix C(N2 , K2 ) along K1 -axis. An updated at previous stage matrix Ẍ(K1 , N2 )K3 isrolled along the “same” K2 -axis. It is clear that P additional “compute-and-roll” time-steps is required
for this final stage and, therefore, totally, 3P such steps are needed for computing of a 3D DXT in a 4D
P × P × P × P orbital index space.

4

A Scalable 3D DXT in 3D Torus Network of Nodes

The next step in the Algorithm/Architecture co-design is mapping a described above 4D logical implementation into the 3D physical space. This mapping is achieved by proper projection of the 4D scheduled orbital
16

Stage I

(a)

N3

N2

X(N1,K3)N2

(b)

N3

N2

C(N3 , K3)
X(N1,N3)N2

N1
K3

(d)

K1

N2

X(N1,K3)N2

K3

Stage II

(c)

K1

N1

N2

X(K1,K3)N2
C(N1 , K1) T

T

N1

N1
K3

T

C(N2 , K2)

K3
Stage III

(e)

K1

(f)

K1

N2

N2

X(K1,K2)K3
X(K1,N2)K3

T

K2

T

K3

K2
K3

Figure 6: A 4D computational index space for the 3D DXT as disjoint sets of the 3D orbital index spaces;
distribution of active index points (filled in black) as well as directions (orbits) of circular data movement
at each “compute-and-roll” time-step (a), (b), (c), (d), (e), (f ) for the scheduling vector α = (1, 1, 1)T and
P = 2; orbits for products accumulation are shown in red; for the given smallest 2 × 2 × 2 torus grid,
wrap-around connections are just opposite to specified direction of data circulation.

17

index space to the 3D computer space where no more than one “active” index point is projected into a single computer node. We used the well-known 3D → 2D linear projection method [28, 31, 45] which isextended here to support a 4D → 3D transform. For the first and second stages of processing, each N2 -th
3D torus index space is projected onto 2D space along K3 -axis (see Fig. 6 (a), (b), (c), (d)) This parallel
mapping results in P different 2D torus networks of P × P computer nodes each such that N2 -slices of a
combined 3D torus of P 3 simple nodes independently and concurrently implements two transforms in 2P“compute-and-roll” time-steps. The third stage is implemented on all P parallel K3 -slices (which are orthogonal to N2 -slices) of this 3D torus of computer nodes by mapping each K3 -th 3D torus index space (seeFig. 6 (e), (f )) into corresponding K3 -slice along K1 -axis. The location of computer nodes in a 3D physicalspace are indicated by triplet (Q, R, S), where 0 ≤ Q, R, S < P. Below, an orbital implementation of the
Forward and Inverse N × N × N linear transforms on a 3D torus of P × P × P computer nodes, where
P = N/b and b ∈ [1, N ] is a blocking factor, is formally described.

S

Q

R

Figure 7: A 4 × 4 × 4 array of computer nodes with toroidal interconnection in a 3D (Q, R, S) space;different colors are used to indicate nodes with distinct initial values of the given orbital function T =(Q + R + S) mod 4 ∈ {0 (red), 1 (blue), 2 (green), 3 (yellow)}.

4.1

3D Block Forward Transform

At the beginning, each computer node CN(Q, R, S) in a P ×P ×P torus array (see Fig. 7 for P = 4) holds
in a local memory the four b × b × b data cubes:• X = X(Q, R, S) = X(N1 , N2 , N3 ),• Ẋ = Ẋ(Q, R, T ) = Ẋ(N1 , N2 , T ) = O,
• Ẍ = Ẍ(S, R, T ) = Ẍ(N3 , N2 , T ) = O,... ...
...
• X = X(S, Q, T ) = X(N3 , N1 , T ) = O,as well as the three b × b change-of-basis matrices of transform coefficients:18

• CI = C(S, T ), CII = C(Q, S), and CIII = C(R, Q),
where O is a cubical b × b × b tensor with all its entries being zero. The initially selected and used modular
function for space-time scheduling of a 4D computational index space defines one of possible functions
T = (αT · p) mod P = (±Q ± R ± S) mod P,
where α = (αR , αQ , αS )T = (±1, ±1, ±1)T is the scheduling vector and p = (Q, R, S)T is a node’s
position in a 3D physical index space (see Fig. 7). Each computer node CN(Q, R, S) in a 3D torus network
has six bi-directional links labeled as ±Q, ±R, ±S. During processing some blocks of tensor and matrix
data are rolled, i.e. cyclically shifted, along (+) or opposite (−) axis (orbit) according to the schedulingvector α (see [46] for more details). As it can be seen from the assignment above, the P × P matrices CI
and CII are replicated among P parallel along R-axis slabs of computer nodes while P × P matrix CIII
is duplicated among P parallel along S-axis slabs. These change-of-basis matrices CI , CII , and CIII are
preloaded into torus array only once for practically possible many 3D transforms.
A three-stage orbital implementation of the 3D forward DXT in a 3-dimensional network of toroidallyinterconnected nodes {CN(Q, R, S) : 0 ≤ Q, R, S < P } under the scheduling function T = (Q + R +S) mod P , i.e. α = (αQ , αR , αS )T = (1, 1, 1)T , is described below (see also Fig. 8 for the smallest case:P = 2, where red colored links show an accumulation or summation direction).
Stage I. Ẋ(Q, R, T ) =P

0≤S<P X(Q, R, S) × C(S, T ) :• for all P 3 CN(Q, R, S) do P times:
1. compute: Ẋ ← X × CI + Ẋ
2. data roll:+S:αS =1

⇐=

−S

Ẋ ⇐=

k

+Q:αQ =1

⇐=

−Q

CI ⇐=

Stage I is an orbital implementation in a physical space of the equation (5) or (9). This implementation is
I (N ), 0 ≤ N < P, which are shownachieved by linear projection of the disjoint 3D torus index spaces I3D2
2
for different time-steps in Fig 6 (a) and (b), onto planar processor or execution space along K3 -axis. Each
I (N ) defines its own N -th 2D torus processor.
mapped 3D torus index space I3D2
2
Stage II. Ẍ(S, R, T ) =P

0≤Q<P C(Q, S)

T × Ẋ(Q, R, T ) :• for all P 3 CN(Q, R, S) do P times:
1. compute: Ẍ ← CIIT × Ẋ + Ẍ
2. data roll:+Q:αQ =1

⇐=

−Q

Ẍ ⇐=

k

+S:αS =1

⇐=

−S

Ẋ ⇐=

Stage II is an orbital implementation of the equation (6) or (11) by projection of the 3D torus index spaces
II (N ), 0 ≤ N < P , shown in Fig 6 (c), (d), onto the same as above 2D processor space along K -axis.
I3D
2
2
3
...
P
Stage III. X(S, Q, T ) = 0≤R<P Ẍ(S, R, T ) × C(R, Q) :• for all P 3 CN(Q, R, S) do P times:...
...
1. compute: X ← Ẍ × CIII + X+Q:αQ =1
+R:αR =1 ... −R−Q
2. data roll: ⇐=⇐= Ẍ ⇐=X ⇐= k19

S

Stage I

S

R

R

Q

(a)
S

Q
Stage II

(b)
S

R

R

Q

(c)
S

Q
Stage III

(d)
S

R

Q

(e)

R

Q

(f)

Figure 8: An orbital or circular data communication between nodes for the 3D FDXT.

Stage III represents an orbital implementation in a physical space of the equation (7) or (13) by using linear
III (K ), 0 ≤ K < P , shown in Fig 6 (e), (f, ) onto a new 2D
projection of the 3D torus index spaces I3D33
execution space along K1 -axis.
Recall that due to the orbital (cyclical or rotational) nature of processing, after completion of eachstage, i.e. after P “compute-and-roll” steps, all rotated data are returned to the same originated nodes and,
therefore, a computed intermediate cubical tensor can be immediately used for the next stage of a cyclical
3D processing. Note that initial tensor X = X(N1 , N2 , N3 ) = X(Q, R, S) is assigned
to the nodes in the...
canonical (in-order) layout whereas all intermediate and final tensors, Ẋ, Ẍ and X, will be distributed in theskewed (out-of-order) layouts. It is also important to note that because at any stage each cubical tensor-bymatrix multiplication is separated into the set of independent matrix-by-matrix products, the “compute” and“data roll” parts can also be splitted and effectively overlapped by software pipelining.

20

S

Stage I

S

R

R

Q

(a)
S

Q
Stage II

(b)
S

R

R

Q

(c)
S

Q
Stage III

(d)
S

R

R

Q

(e)

Q

(f)

Figure 9: An orbital or circular data communication between nodes for the 3D IDXT.

The first two stages implement the set of P space-independent 2D FDXTs on P parallel along R-axis(orbit) slabs, 0 ≤ R < P , with the P × P toroidally interconnected computer nodes {CN(Q, ∗, S)R :0 ≤ Q, S < P } each (see Fig. 8 for the stages I and II). Totally, 2P “compute-and-roll” time-steps areneeded for each R-th slab-of-nodes to independently implement a 2D FDXT (stages I and II) of the R-thslab-of-cubes X(N1 , N3 )R=N2 .
As it was discussed, to make data distribution of the intermediate cubical tensor Ẍ consistent with therequired distribution of this tensor on Stage III, each R-th 2D FDXT is implemented in time-skewed mannersuch that after 2P time-steps, the b×b×b blocks of a cubical tensor Ẍ will be correctly distributed among all
P 3 computer nodes in a 3D torus network. This required initial, intermediate and final data distributions arecontrolled by the modular (orbital) scheduling function T . By using this scheduling, no data redistribution
is needed during processing and the final P independent sets of 1D FDXTs in Stage III can be immediately
started by P, now parallel along S-axis (0 ≤ S < P ), slabs of the P ×P toroidally interconnected computer
nodes {CN(Q, R, ∗)S : 0 ≤ Q, R < P } each (see Stage III in Fig. 8 for P = 2).
21

4.2

3D Block Inverse Transform

An orbital computing of the 3D inverse DXT is implemented as rolling back of...the described above 3D
forward DXT where skewed (out-of-order) distribution of the final cubical tensor X among computer nodes
will correspond to the initial three-way tensor distribution for the 3D IDXT (see also Fig. 9 for P = 2).
Obviously, the change-of-basis matrices CI , CII , and CIII are the same as for the 3D FDXT.
Stage I. Ẍ(S, R, T ) =...

P

0≤Q<P X(S, Q, T ) × C(R, Q)

T :• for all P 3 CN(Q, R, S) do P times:...
T + Ẍ
1. compute: Ẍ ← X × CIII
+R ... −R+Q−Q
2. data roll: ⇐= X ⇐= k ⇐= Ẍ ⇐=Stage II. Ẋ(Q, R, T ) =P

0≤S<P C(Q, S) × Ẍ(S, R, T ) :• for all P 3 CN(Q, R, S) do P times:
1. compute: Ẋ ← CII × Ẍ + Ẋ−Q+Q

2. data roll: ⇐= Ẍ ⇐=Stage III. X(Q, R, S) =+S−Sk ⇐= Ẋ ⇐=P

0≤T <P Ẋ(Q, R, T ) × C(S, T )T :• for all P 3 CN(Q, R, S) do P times:
1. compute: X ← Ẋ × CIT + X+S−S

2. data roll: ⇐= Ẋ ⇐=+Q−Qk ⇐= CI ⇐=Note that at the last stage of the 3D IDXT, an accumulation or T -summation is implemented inside eachcomputer node. After 3P “compute-and-roll” time-steps, i.e. after completion of a 3D IDXT, the resulting
cubical tensor X will be distributed in a 3D processor space in the canonical (in-order) layout. As a consequence, the frequently required by many real-world applications Forward/Backward looping, 3D FDXT
3D IDXT, can be implemented on a 3D torus network of computer nodes in a very natural and efficient way.

4.3

Complexity and Scalability Analysis

An orbital implementation of the forward or inverse block N × N × N 3D DXT on a P × P × P torus
network of computer nodes requires totally 3P block “compute-and-roll” time-steps, where P = N/b and1 ≤ b ≤ N is the blocking factor. Each “compute-and-roll” time-step involves the left or right multiplicationof a 3D b × b × b tensor by a b × b coefficient matrix with a possible matrix transpose. In turn, each of thistensor-by-matrix multiplication requires execution by every computer node exactly b4 scalar fused multiplyadd (fma) operations and rolling of either b3 + b2 or 2b3 scalar data. Each node implements a data rollingby concurrently sending and receiving data to/from two predefined nearest-neighbor nodes. Actually, eachcomputer node doesn’t “feel” the size of a network, because at each time-step only local computing andtopologically nearby (contact) data movement operations are implemented which can be done independent
from the total size of a system. In this respect and by analogy with Strassen’s 2 × 2 matrix multiplication22

algorithm 2 , it is enough to organize a parallel data processing for the smallest, but fundamentally basic,
2 × 2 × 2 case which can be recursively scaled then to the bigger sizes with a proportional escalating ofalgorithm’s time and space complexity.
The required memory space in one computer node, which is proportionally to 4b3 + 3b2 , might be sufficiently reduced by using the well-known “compute-in-place” approach. Note that for simplicity, we assumed
here that there is no shared memory between nearest-neighbor nodes, i.e., after computing, a local data
must be exchanged (rolled) by using a message-passing. However, if upgrading of data (a “compute” part)precedes a movement of this data (a “roll” part) and there exist a shared memory between nearest-neighbornodes then data upgrading can be directly done in this shared memory, which will replace a time-consuming,
message-passing, “roll” part of this data by quick inter-node synchronization. For the other unchangeable,
but moved data, “roll” part(s) can be implemented concurrently with a “compute” part (see description of
the 3D DXT algorithms above), making all nearest-neighbor data movements totally overlapped with a local
computing. This way of sharable data exchange should be very promising for massively-parallel solution ofrelatively small problems.
If blocking factor b = 1, i.e. P = N , our 3D transform algorithms achieve the highest or extreme
degree of parallelization with a total processing time equals to 3N scalar “compute-and-roll” steps on an
N × N × N torus array of very simple computer nodes. This fastest, fine-grained implementation will
require an fma-unit in each computer node to be with one-step or one-cycle latency which is possible for
fixed-point arithmetic, like in “Anton” supercomputer [49]. Moreover, only in this extreme case, the sizeof required in each computer node memory is independent from the size of problem. For floating-point
arithmetic, however, the pipelined fma-units usually require a few cycles of latency as in today advanced
microprocessors [42, 20, 10, 8, 38, 37], which force the size of blocking factor b to be more then one inorder to hide this pipeline latency by concurrent execution of a few independent fma-instructions on thesame fma-unit (usually, b ≥ 4).
If, in contrast, a blocking factor b = N , i.e. P = 1, the proposed implementation of a 3D transform isdegenerated to a serial algorithm which is executed on a single computer node in 3N 4 fma-steps. It is clearthat independently on the size of the blocking factor b ∈ [1, 2, ..., N ] and, therefore, degree of parallelization,
any implementation of a 3D transform in the form of three-way tensor-by-matrix product requires executionof 3N 4 fma-operations.

5

Conclusions

We have systematically designed massively-parallel and extremely-scalable orbital block algorithms to efficiently perform any 3D separable transform and its inverse in a 3D network of toroidally interconnectedcomputer nodes. The designed orbital algorithms preserve all advantages of the previously popular systolic array processing [28, 26] such as simplicity, regularity, nearest-neighbor data movement, scalable
parallelism, pipelining, etc. However, unlike systolic processing, our orbital computing keeps all initial,
intermediate, and final data inside the same 3D torus array processor allowing a simultaneous access to all
data, reusing data by its circulation, and an effective space-time chaining of different parts of processing.
Our basic 3D DXT algorithms are structurally represented in a 4-dimensional computational index space
which is partitioned into the 3D disjoint sets of index points. A modular (orbital) scheduling of these index
points for execution and data movement as well as appropriate 4D→3D projections allow to map a complex three-stage data-dependent processing into a single 3D physical network of toroidally interconnected
2

In Strassen’s algorithm, the time of a 2 × 2 matrix-by-matrix multiplication is decreased by reducing the number of multiplications in the price of increasing the number of more cheap addition/subtractions operations and irregular data access; in our case,however, we directly reduce this time by overlapping (indivisible) multiply-add operations combined with a local and regular data
movement.

23

nodes. A multidimensional torus interconnection has always been widely used in the past and present
distributed-memory supercomputers, for example, in Cray T3E (3D torus) [44], IBM Blue Gene/L (3Dtorus) [1] and Blue Gene/Q (5D torus), D.E. Shaw Research “Anton” (3D torus) [49], Fujitsu K-computer(6D torus) [18, 3]. It is expected that direct implementation of our coarse-grained (b  1) forward andinverse three-dimensional DXT algorithms on supercomputers with a multidimensional torus interconnect
would be beneficial with respect to other existing implementations which are based on a 1D or 2D data
decomposition. We reserve this porting and performance evaluation as a future work. Note that instead ofusing widely available and deeply optimized “fast”, but internally serial, DXT algorithms (as in 1D or 2D
decomposition approaches), our 3D decomposition technique allows to use within each node, also widely
available and practically even more deeply optimized and, therefore, more fast and efficient, GEMM-basedalgorithms [29, 36].
Moreover, our extremely scalable solution (b ∼= 1) for a fine-grained 3D DXT implementation on a
cubical “mesh-of-tori” network 3 [43] of very simple execution nodes (each with one multiply-add unit anda few, inter-node connected, registers) is perfectly suited for forthcoming 3-dimensional VLSI technology[25], Giga-size sensor arrays with parallel read-out [21] and many real-world applications (including embedded) with multidimensional data. Note that preliminary, but promising results of FPGA implementationof a 3D torus array processor for the fine- and coarse-grained 3D Discrete Cosine Transform (3D DCT) has
recently been reported [23].
One of the notable characteristics of our unified and extremely scalable 3D DXT algorithms is in thepossibility to implement concurrently a few distinct transforms on the same cubical data by using different
matrices of transform coefficients. Additionally, if it is required by application, the different transforms fordifferent cubical data dimensions can be easily implemented at the same time.

References[1] N. R. Adiga, M. A. Blumrich, D. Chen, P. Coteus, A. Gara, M. E. Giampapa, P. Heidelberger,
S. Singh, B. D. Steinmacher-Burow, T. Takken, M. Tsao, and P. Vranas, “Blue Gene/L torus
interconnection network,” IBM J. Res. Dev., vol. 49, pp. 265–276, March 2005. [Online]. Available:
http://dx.doi.org/10.1147/rd.492.0265[2] R. Agarwal, F. Gustavson, and M. Zubair, “A high performance matrix multiplication algorithm on a
distributed-memory parallel computer, using overlapped communication,” IBM J. of Res. and Develop.,
vol. 38, no. 6, pp. 673–681, 1994.
[3] Y. Ajima, S. Sumimoto, and T. Shimizu, “Tofu: A 6D Mesh/Torus Interconnect for Exascale Computers,” Computer, vol. 42, pp. 36–40, 2009.
[4] B. W. Bader and T. G. Kolda, “Algorithm 862: Matlab tensor classes for fast algorithm
prototyping,” ACM Trans. Math. Softw., vol. 32, pp. 635–653, December 2006. [Online]. Available:
http://doi.acm.org/10.1145/1186785.1186794[5] G. Baumgartner, A. Auer, D. Bernholdt, A. Bibireata, V. Choppella, D. Cociorva, X. Gao, R. Harrison,
S. Hirata, S. Krishnamoorthy, S. Krishnan, C. Lam, Q. Lu, M. Nooijen, R. Pitzer, J. Ramanujam,
P. Sadayappan, and A. Sibiryakov, “Synthesis of high-performance parallel programs for a class of abinitio quantum chemistry models,” Proceedings of the IEEE, vol. 93, no. 2, pp. 276 –292, feb. 2005.
3

i.e. a tile-based regular network without “long” wrap-around connections.

24[6] E. Brachos, “Parallel FFT Libraries,” Master’s thesis, The University of Edinburgh, 2011. [Online].
Available: http://www.brachos.gr/files/E-Brachos FFT.pdf[7] L. Cannon, “A cellular computer to implement the Kalman filter algorithm,” Ph.D. dissertation, Montana State University, 1969.
[8] S. Chatterjee, L. R. Bachega, P. Bergner, K. A. Dockser, J. A. Gunnels, M. Gupta, F. G. Gustavson,
C. A. Lapkowski, G. K. Liu, M. Mendell, R. Nair, C. D. Wait, T. J. C. Ward, and P. Wu, “Designand exploitation of a high-performance SIMD floating-point unit for Blue Gene/L,” IBM Journal of
Research and Development, vol. 49, no. 2.3, pp. 377 –391, march 2005.
[9] J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series,”
Mathematics of Computation, vol. 19, pp. 297–301, 1965.
[10] M. Cornea, J. Harrison, and P. T. P. Tang, “Intel Itanium floating-point architecture,” in Proceedingsof the 2003 workshop on Computer architecture education: Held in conjunction with the 30th
International Symposium on Computer Architecture, ser. WCAE ’03. New York, NY, USA: ACM,
2003. [Online]. Available: http://doi.acm.org/10.1145/1275521.1275526[11] J. J. Dongarra, J. Du Croz, S. Hammarling, and I. S. Duff, “A set of level 3 basic linear algebra
subprograms,” ACM Trans. Math. Softw., vol. 16, pp. 1–17, March 1990. [Online]. Available:
http://doi.acm.org/10.1145/77626.79170[12] M. Eleftheriou, B. G. Fitch, A. Rayshubskiy, T. J. C. Ward, and R. S. Germain, “Performance measurements of the 3D FFT on the Blue Gene/L supercomputer,” in Euro-Par, ser. Lecture Notes in Computer
Science, J. C. Cunha and P. D. Medeiros, Eds., vol. 3648. Springer, 2005, pp. 795–803.
[13] M. Eleftheriou, J. E. Moreira, B. G. Fitch, and R. S. Germain, “A Volumetric FFT for BlueGene/L,”
in HiPC, ser. Lecture Notes in Computer Science, T. M. Pinkston and V. K. Prasanna, Eds., vol. 2913.
Springer, 2003, pp. 194–203.
[14] S. Filippone, “The IBM parallel engineering and scientific subroutine library,” in PARA, ser. Lecture
Notes in Computer Science, J. Dongarra, K. Madsen, and J. Wasniewski, Eds., vol. 1041. Springer,
1995, pp. 199–206.
[15] G. Fox, S. Otto, and A. Hey, “Matrix algorithms on a hypercube I: Matrix multiplication,” Parallel
Computing, vol. 4, pp. 17–31, 1987.
[16] M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE,
vol. 93, no. 2, pp. 216–231, 2005, special issue on “Program Generation, Optimization, and Platform
Adaptation”.
[17] ——. (2012, June) FFTW for version 3.3.2. [Online]. Available: http://www.fftw.org[18] FUJITSU, “K computer,” Website, June, 2012. [Online]. Available: http://www.fujitsu.com/global/
about/tech/k/[19] G.H.Golub and C. Loan, Matrix Computations (3rd ed.).

Johns Hopkins University Press, 1996.

[20] F. G. Gustavson, J. E. Moreira, and R. F. Enenkel, “The fused multiply-add instruction leads
to algorithms for extended-precision floating point: applications to java and high-performancecomputing,” in Proceedings of the 1999 conference of the Centre for Advanced Studies on25

Collaborative research, ser. CASCON ’99. IBM Press, 1999, pp. 4–. [Online]. Available:
http://portal.acm.org/citation.cfm?id=781995.781999[21] E. H. Heijne, “Gigasensors for an Attoscope: Catching Quanta in CMOS,” IEEE Solid-State Circuits
Newsletter, vol. 13, pp. 28–34, 2008.
[22] C. Hillar and L.-H. Lim, “Most tensor problems are NP hard,” 2009. [Online]. Available:
http://www.citebase.org/abstract?id=oai:arXiv.org:0911.1393[23] Y. Ikegaki, T. Miyazaki, and S. G. Sedukhin, “3D-DCT Processor and its FPGA Implementation,”
IEICE Transactions on Information and Systems, vol. E94.D, no. 7, pp. 1409–1418, 2011. [Online].
Available: http://dx.doi.org/10.1587/transinf.E94.D.1409[24] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3,
pp. 455–500, September 2009.
[25] M. Koyanagi, T. Fukushima, and T. Tanaka, “Three-dimensional integration technology and integrated
systems,” in Design Automation Conference, 2009. ASP-DAC 2009. Asia and South Pacific, Jan. 2009,
pp. 409–415.
[26] H. T. Kung, “Why systolic architectures?” Computer, vol. 15, no. 1, pp. 37–46, Jan. 1982. [Online].
Available: http://dx.doi.org/10.1109/MC.1982.1653825[27] H. Kung and C. Leiserson, “Systolic array apparatuses for matrix computations,” Patent US 4,493,048,
Jan. 8, 1985. [Online]. Available: http://www.patentlens.net/patentlens/patent/US 4493048/[28] S. Kung, VLSI Array Processors.

Prentice Hall, 1988.

[29] J. Kurzak, S. Tomov, and J. Dongarra, “Autotuning GEMMs for Fermi,” University ofTennessee, Tech. Rep. UT-CS11-671, April 2011. [Online]. Available: http://icl.cs.utk.edu/news pub/submissions/lawn245.pdf[30] L. D. Lathauwer, B. D. Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM
J. Matrix Anal. Appl, vol. 21, pp. 1253–1278, 2000.
[31] D. Lavenier, P. Quinton, and S. Rajopadhye, “Advanced systolic design,” in Digital Signal Processing
for Multimedia systems, ser. Signal Processing Series. Marcel Dekker, 1999, ch. 23, pp. 657–692.
[32] N. Li and S. Laizet, “2DECOMP&FFT - A Highly Scalable 2D Decomposition: Library andFFT Interface,” in Cray User Group 2010 conference, 2010, pp. 1–13. [Online]. Available:
http://www.2decomp.org/pdf/17B-CUG2010-paper-Ning Li.pdf[33] L.-H. Lim, “Numerical Multilinear Algebra: a new beginning?” in Matrix Computations andScientific Computing Seminar, 2006, pp. 1–51. [Online]. Available: http://galton.uchicago.edu/∼lekheng/work/mcsc1.pdf[34] O. Lindtjorn, R. Clapp, O. Pell, H. Fu, M. Flynn, and O. Mencer, “Beyond traditional microprocessors
for geoscience high-performance computing applications,” IEEE Micro, vol. 31, pp. 41–49, March
2011. [Online]. Available: http://dx.doi.org/10.1109/MM.2011.17

26

[35] Q. Lu, X. Gao, S. Krishnamoorthy, G. Baumgartner, J. Ramanujam, and P. Sadayappan, “Empiricalperformance model-driven data layout optimization and library call selection for tensor contraction
expressions,” J. Parallel Distrib. Comput., vol. 72, no. 3, pp. 338–352, Mar. 2012. [Online]. Available:
http://dx.doi.org/10.1016/j.jpdc.2011.09.006[36] K. Matsumoto, N. Nakasato, T. Sakai, H. Yahagi, and S. G. Sedukhin, “Multi-level optimization of
matrix multiplication for GPU-equipped systems,” Procedia CS, vol. 4, pp. 342–351, 2011.
[37] R. K. Montoye, E. Hokenek, and S. L. Runyon, “Design of the IBM RISC System/6000
floating-point execution unit,” IBM J. Res. Dev., vol. 34, pp. 59–70, January 1990. [Online]. Available:
http://dx.doi.org/10.1147/rd.341.0059[38] H.-J. Oh, S. Mueller, C. Jacobi, K. Tran, S. Cottier, B. Michael, H. Nishikawa, Y. Totsuka, T. Namatame, N. Yano, T. Machida, and S. Dhong, “A fully pipelined single-precision floating-point unit in
the synergistic processor element of a CELL processor,” Solid-State Circuits, IEEE Journal of, vol. 41,
no. 4, pp. 759 – 771, april 2006.
[39] I. V. Oseledets, D. V. Savostyanov, and E. E. Tyrtyshnikov, “Linear algebra for tensor problems,”
Computing, vol. 85, pp. 169–188, July 2009. [Online]. Available: http://dx.doi.org/10.1007/
s00607-009-0047-6[40] D. Pekurovsky, “P3DFFT,” Website, August, 2011. [Online]. Available: http://www.sdsc.edu/us/
resources/p3dfft/[41] M. Pippig and D. Potts, “Scaling Parallel Fast Fourier Transform on BlueGene/P,” Jülich Supercomputing Centre, Tech. Rep. FZJ-JSC-IB-2010-03, May 2010.
[42] E. Quinnell, E. Swartzlander, and C. Lemonds, “Floating-point fused multiply-add architectures,” in
Signals, Systems and Computers, 2007. ACSSC 2007. Conference Record of the Forty-First Asilomar
Conference on, nov. 2007, pp. 331 –337. [Online]. Available: http://dx.doi.org/10.1109/ACSSC.2007.
4487224[43] A. A. Ravankar and S. G. Sedukhin, “Mesh-of-tori: A novel interconnection network for frontal plane
cellular processors,” in Proceedings of the 2010 First International Conference on Networking andComputing, ser. ICNC ’10. Washington, DC, USA: IEEE Computer Society, 2010, pp. 281–284.
[Online]. Available: http://dx.doi.org/10.1109/IC-NC.2010.30[44] S. Scott and G. Thorson, “The Cray T3E network: Adaptive routing in a high performance 3D torus,”
in Proceedings of HOT Interconnects IV, 1996, pp. 147–156.
[45] S. G. Sedukhin and I. S. Sedukhin, “Systematic approach and software tool for systolic design,” in
CONPAR, ser. Lecture Notes in Computer Science, B. Buchberger and J. Volkert, Eds., vol. 854.
Springer, 1994, pp. 172–183.
[46] S. G. Sedukhin, A. S. Zekri, and T. Myiazaki, “Orbital algorithms and unified array processor forcomputing 2D separable transforms,” Parallel Processing Workshops, International Conference on,vol. 0, pp. 127–134, 2010. [Online]. Available: http://dx.doi.org/10.1109/ICPPW.2010.29[47] E. Solomonik, J. Hammond, and J. Demmel, “A preliminary analysis of cyclops tensor framework,”
EECS Department, University of California, Berkeley, Tech. Rep. UCB/EECS-2012-29, Mar 2012.
[Online]. Available: http://www.eecs.berkeley.edu/Pubs/TechRpts/2012/EECS-2012-29.html
27

[48] R. van de Geijn and J. Watts, “SUMMA: scalable universal matrix multiplication algorithm,” The
University of Texas, Tech. Rep. TR-95-13, Apr. 1995.
[49] C. Young, J. A. Bank, R. O. Dror, J. P. Grossman, J. K. Salmon, and D. E. Shaw, “A 32x32x32,
spatially distributed 3D FFT in four microseconds on Anton,” in Proceedings of the Conference on
High Performance Computing Networking, Storage and Analysis, ser. SC ’09. New York, NY, USA:ACM, 2009, pp. 23:1–23:11. [Online]. Available: http://doi.acm.org/10.1145/1654059.1654083

28