Optimized M2L Kernels for the Chebyshev Interpolation based Fast
Multipole Method

arXiv:1210.7292v2 [cs.NA] 20 Nov 2012

Matthias Messner, Berenger Bramas, Olivier Coulaud, Eric Darve
November 21, 2012
Abstract
A fast multipole method (FMM) for asymptotically smooth kernel functions (1/r, 1/r4 , Gauss and
Stokes kernels, radial basis functions, etc.) based on a Chebyshev interpolation scheme has been introduced in [5]. The method has been extended to oscillatory kernels (e.g., Helmholtz kernel) in [12].
Beside its generality this FMM turns out to be favorable due to its easy implementation and its high
performance based on intensive use of highly optimized BLAS libraries. However, one of its bottlenecks
is the precomputation of the multiple-to-local (M2L) operator, and its higher number of floating point
operations (flops) compared to other FMM formulations. Here, we present several optimizations for that
operator, which is known to be the costliest FMM operator. The most efficient ones do not only reduce
the precomputation time by a factor up to 340 but they also speed up the matrix-vector product. We
conclude with comparisons and numerical validations of all presented optimizations.

1

Introduction

The fast multipole method (FMM) is a method first designed in [7] to reduce the cost of matrix-vector
products from O(N 2 ) to O(N ) or O(N log N ) depending on the underlying kernel function. Most FMM
variants have been developed and optimized for specific kernel functions [13, 8, 3, 2]. However, some have
a formulation that is independent of the kernel function [14, 10, 6, 5, 12]. The paper at hand addresses the
optimization of one of these formulations, the so called black-box FMM (bbFMM), presented in [5]. It is
based on the approximation of the kernel function via a Chebyshev interpolation and is a black-box scheme
for kernel functions that are asymptotically smooth, e.g., 1/(r2 + c2 )n/2 with r = |x − y|, c a given parameter
and n ∈ N. The bbFMM has been extended to the directional FMM (dFMM) for oscillatory kernels in [12].
It is suitable for any kernel function of the type g(r)eıkr where g(r) is an asymptotically smooth function
(ı2 = −1 is the imaginary unit and k the wavenumber).
The main idea of the FMM is to properly separate near-field (|x − y| → 0) and far-field (|x − y| → ∞).
The near-field is evaluated directly and the far-field can be approximated and thus computed efficiently. In
this paper, we will denote M2L the operator that transforms a multipole expansion into a local expansion.
The M2L operator is the costliest step in the method, since it needs to be applied to all clusters in the
interaction list, that is 189 times for each cluster for a typical configuration in bbFMM. In this paper we
focus on various optimizations of this operator for both the bbFMM and the dFMM.
First, we address the optimization proposed in [5, 12]. In that paper, the M2L operator between a cluster
and its interaction list is further compressed using a singular value decomposition (SVD). The key idea is that
the SVD provides the smallest possible rank given an error ε (therefore leading to the smallest computational
cost). In [5], a single SVD is used to compress all the M2L operators at a given level using a single
singular vector basis. However, the singular values of individual M2L operators decay at different rates; it
is often the case for example that cluster pairs that are separated by a small distance have higher rank than
pairs that are separated by a greater distance. If we denote w the diameter (or width) of a cluster and d
the distance
√ smallest distance corresponds to dmin = 2w while the largest is
√ between two clusters, then the
dmax = 3 3w. The ratio dmax /dmin = 3 3/2 is in fact relatively large. This can be taken advantage of to
1

further compress the M2L operator on a cluster-pair basis leading to an even smaller computational cost
than the original algorithm [5]. This point is investigated in details in this paper.
Another bottleneck in the FMM of [5] is the precomputation time of the M2L operators. We therefore
introduce a new set of optimizations that exploit symmetries in the arrangement of the M2L operators. For
example, for bbFMM, this allows us to express all M2L operators (316 = 73 − 33 ) as permutations of a subset
with only 16 unique M2L operators. These permutations correspond to reflections across various planes (for
example x = 0, y = 0, z = 0, etc). Table 5 for example reports reductions in precomputation time by a
factor of about 50–200x.
Besides drastically reducing precomputation time and memory requirement this approach paves also the
road for various algorithmic improvements. Modern processors use memory cache to enhance performance
and, as a result, schemes that block data (thereby reducing the amount of data traffic between main memory
and the floating point units) can result in much higher performance. In our case, the use of symmetries
allows calling highly optimized matrix-matrix product implementations (instead of matrix-vector) that have
better cache reuse. Overall, this may result in acceleration by a factor 4–6x for smooth kernels (Tables 6–8).
We also present results using oscillatory kernels of the type mentioned above (g(r)eıkr ). In that case, the
acceleration is more modest and ranges from about 1.2 to 2.7x (Table 9–11).
In this paper we therefore focus both on performance and precomputation time. Both are important
factors for the FMM: the precomputation is crucial if we are interested in only one matrix-vector product,
while in other cases, such as the solution of a linear system, the fast application of matrix-vector products
is central.
The paper is organized as follows. In Sec. 2, we briefly recall the bbFMM and the dFMM and introduce
notations that are needed for explanations later in this paper. In Sec. 3, we address the separation of nearand far-field, introduce the notion of transfer vector to uniquely identify M2L operators, and describe how
the kernel (smooth or oscillatory) affects the interaction list. We start Sec. 4 with a brief recall of the known
optimization of the M2L operator (see [5, 12]) and suggest further improvements. Then, we present a new
set of optimizations and explain them in details. These include using a low-rank approximation of M2L for
individual interactions (individual cluster pairs), along with the definition of symmetries and permutations
and how they are used to reduce the computational cost. Finally, in Sec. 5, we present numerical results.
We compare all variants for bbFMM and dFMM with three datasets corresponding to a sphere, an oblate
spheroid, and a prolate spheroid. The measured running time of the FMM, along with the precomputation
time, and its accuracy as a function of various parameters are reported. The efficiency of the algorithms
as a function of the target accuracy is considered. We also performed tests using two different kinds of
linear algebra libraries for the matrix-matrix products and other operations (an implementation of BLAS
and LAPACK vs the Intel MKL library).

2

Pairwise particle interactions

The problem statement reads as follows. Assume, the cells X and Y contain source {yj }N
j=1 and target
M
particles {xi }i=1 , respectively. Compute the interactions
fi =

N
X

K(xi , yj )wj

for i, . . . , M.

(1)

j=1

The kernel function K(x, y) describes the influence of the source particles onto the target particles. The cost
of directly evaluating the summation in Eqn. (1) grows like O(M N ) which becomes prohibitive as M, N → ∞
and it is why we need a fast summation scheme.

2.1

Fast summation schemes based on Chebyshev interpolation

For a detailed derivation and error analysis of the FMM based on Chebyshev interpolation we refer the reader
to [5, 12]. We adapt most of their notations and repeat only concepts which are necessary to understand
2

explanations hereafter.
Let the function f : R3 → C be approximated by a Chebyshev interpolation scheme as
X
Sℓ (x, x̄α )f (x̄α )
f (x) ∼

(2)

|α|≤3ℓ

with the 3-dimensional multi-index α = (α1 , α2 , α3 ) and |α| = max(α1 , α2 , α3 ) with αi ∈ (1, . . . , ℓ). The
interpolation points x̄α = (x̄α1 , x̄α2 , x̄α3 ) are given by the tensor-product of the Chebyshev roots x̄αi of the
Chebyshev polynomial of first kind Tℓ (x) = cos(arccos x) with x ∈ [−1, 1]. The interpolation operator reads
as
(3)
Sℓ (x, x̄α ) = Sℓ (x1 , x̄α1 ) Sℓ (x2 , x̄α2 ) Sℓ (x3 , x̄α3 ).
For the interpolation on arbitrary intervals, we need the affine mapping Φ : [−1, 1] → [a, b]. We omit it
hereafter for the sake of readability.
2.1.1

Black–box FMM (bbFMM)

If two cells X and Y are well separated, we know from [5] that asymptotically smooth kernel functions can
be interpolated as
X
X
K(x̄α , ȳβ ) Sℓ (y, ȳβ ).
(4)
Sℓ (x, x̄α )
K(x, y) ∼
|β|≤ℓ

|α|≤ℓ

We insert the above approximation into Eqn. (1) and obtain
fi ∼

X

Sℓ (xi , x̄α )

X

K(x̄α , ȳβ )

Sℓ (yj , ȳβ ) wj

(5)

j=1

|β|≤ℓ

|α|≤ℓ

N
X

which we split up in a three-stage fast summation scheme.
1. Particle to moment (P2M) or moment to moment (M2M) operator: equivalent source values are
anterpolated at the interpolation points ȳβ ∈ Y by
Wβ =

N
X

S(yj , ȳβ ) wj

j=1

for

|β| ≤ ℓ.

(6)

2. Moment to local operator (M2L): target values are evaluated at the interpolation points x̄α ∈ X by
X
K(x̄α , ȳβ ) Wβ for |α| ≤ ℓ.
(7)
Fα =
|β|≤ℓ

3. Local to local (L2L) or local to particle (L2P) operator: target values are interpolated at final points
xi ∈ X by
X
S(xi , x̄α ) Fα for i = 1, . . . , M.
(8)
fi ∼
|α|≤ℓ

Recall, the cells X and Y are well separated and thus all contributions of fi can be computed via the above
presented fast summation scheme (no direct summation is necessary).

3

2.1.2

Directional FMM (dFMM)

Whenever we deal with oscillatory kernel functions, such as the Helmholtz kernel, the wavenumber k comes
into play. Depending on the diameter of the cells X and Y and the wavenumber they are either in the
low-frequency or in the high-frequency regime. In the low-frequency regime the fast summation schemes of
the dFMM and the bbFMM are the same. In the high-frequency regime the fast summation scheme becomes
directional. From [12] we know that any oscillatory kernel function K(x, y) = G(x, y)eık|x−y| , where G(x, y)
is an asymptotically smooth function, can be rewritten as
K(x, y) = K u (x, y)eıku·(x−y)

K u (x, y) = G(x, y)eık(|x−y|−u·(x−y)) .

with

(9)

We assume that the cells X and Y of width w are centered at cx and cy and cy lies in a cone of direction u
being centered at cx (think of the domain around X being virtually subdivided in cones given by directional
unit vectors {uc }C
c=1 , where C is determined by their aperture). If the cell pair (X, Y ) satisfies the separation
criterion O(kw2 ) and the cone-aperture criterion O(1/kw), the error of the Chebyshev interpolation of the
kernel function
X
X
K u (x̄α , ȳβ ) Sℓ (y, ȳβ )
(10)
Sℓ (x, x̄α )
K u (x, y) ∼
|β|≤ℓ

|α|≤ℓ

decays exponentially in the interpolation order ℓ (independent of the wavenumber k; see [11, 12]). We insert
the above interpolated kernel function in Eqn. (1) and obtain
fi ∼ eıku·xi

X

Sℓ (xi , x̄α )e−ıku·x̄α

X

K(x̄α , ȳβ )eıku·ȳβ

Sℓ (yj , ȳβ )e−ıku·yj wj

(11)

j=1

|β|≤ℓ

|α|≤ℓ

N
X

for all i = 1, . . . , M . Similarly as with the bbFMM, a three-stage fast summation scheme for oscillatory
kernels in the high-frequency regime can be constructed.
1. Particle to moment (P2M) or moment to moment (M2M) operator: equivalent source values are
anterpolated at the interpolation points ȳβ ∈ Y by
Wβu = eıku·ȳβ

N
X

S(yj , ȳβ ) e−ıku·yj wj

j=1

for |β| ≤ ℓ.

(12)

2. Moment to local operator (M2L): target values are evaluated at the interpolation points x̄α ∈ X by
X
K(x̄α , ȳβ ) Wβu for |α| ≤ ℓ.
(13)
Fαu =
|β|≤ℓ

3. Local to local (L2L) or local to particle (L2P) operator: target values are interpolated at final points
xi ∈ X by
X
(14)
S(xi , x̄α ) e−ıku·x̄α Fαu for i = 1, . . . , M.
fi ∼ eıku·xi
|α|≤ℓ

Even though the bbFMM and the dFMM are here presented as single-level schemes, they are usually
implemented as multilevel schemes. Strictly speaking, the steps one and three of both schemes are the P2M
and L2P operators. Let us recall briefly on the directional M2M and L2L operators of dFMM: based on
the criterion O(1/kw), the aperture of the cones at the upper level is about half the aperture at the lower
level. Due to a nested cone construction along octree levels, we preserve the accuracy of the Chebyshev
interpolation within the multilevel scheme. For a detailed description of all operators we refer to [5, 12]).
Note, the similarity of the M2L operators (step two of both schemes) of bbFMM and dFMM. In fact, the
only difference in their implementation is that in the bbFMM case we have one loop over all cell pairs,
whereas in the dFMM case we have two loops: the outer loop over all existing cones of direction {uc }C
c=1
and the inner loop over all cell pairs lying in the current cone. In this paper, we focus on the M2L operators
and their efficient numerical treatment.
4

3

M2L operators

The first step of any FMM consists in a proper separation of near- and far-field. After that, the near-field
is evaluated directly and the far-field efficiently using a fast summation scheme. In this section, we focus
on the first step. The union of near- and far-field of a target cell X is spatially restricted to the near-field
of its parent cell. Algorithm 1 explains how these interactions are computed for dFMM [12]. The recursive
partitioning starts with the two root cells X and Y of the octrees for source and target particles. If a pair
of cells satisfies the separation criterion in the high- or low-frequency regime, Y is a far-field interaction of
X. Else, if they are at the leaf level of the octree, Y is a near-field interaction of X. If none is true, the cell
is subdivided and the tests are repeated recursively. In line 3 in Alg. 1 we use the term directional far-field,
Algorithm 1 Separate near- and far-field in the low- and high-frequency regime
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:

function SeparateNearAndFarField(X, Y )
if (X, Y ) are admissible in the high-frequency regime then
add Y to the directional far-field of X return
else if (X, Y ) are admissible in the low-frequency regime then
add Y to the far-field of X return
else if (X, Y ) are leaves then
add Y to the near-field of X return
else
for all Xchild ∈ X and all Ychild ∈ Y do
SeparateNearAndFarField(Xchild , Ychild )
end for
end if
end function

a concept explained in detail in [12]: In the high-frequency regime the far-field is subdivided into cones of
direction u needed by the directional kernel function K u (x, y). Each source cell Y is assigned to a cone and
there are as many directional far-fields as there are cones.

3.1

Transfer vectors

In order to address interactions uniquely we introduce transfer vectors t = (t1 , t2 , t3 ) with t ∈ Z3 . They
describe the relative positioning of two cells X and Y and are computed as t = (cx −cy )/w where cx and cy
denote the centers of the cells X and Y and w their width. In the following we use transfer vectors to
uniquely identify the M2L operator that computes the interaction between a target cell X and a source cell
Y . In matrix notation an M2L operator reads as Kt of size ℓ3 × ℓ3 and the entries are computed as
(Kt )m(α)n(β) = K(x̄α , ȳβ )

(15)

with the interpolation points x̄α ∈ X and ȳβ ∈ Y . Bijective mappings
m, n : {1, . . . , ℓ} × {1, . . . , ℓ} × {1, . . . , ℓ} → {1, . . . , ℓ3 }

(16)

with m−1 (m(α)) = α and n−1 (n(β)) = β provide unique indices to map back and forth between multi-indices
of x̄α and ȳβ and rows and columns of Kt . We choose them to be m(α) = (α1 − 1) + (α2 − 1)ℓ + (α3 − 1)ℓ2 + 1
and same for n(β).

3.2

Interaction lists

Very commonly fast multipole methods are used for translation invariant kernel functions K(x, y) = K(x +
v, y + v) for any v ∈ R3 . Because of that and because of the regular arrangement of interpolation points x̄
and ȳ in uniform octrees it is sufficient to identify unique transfer vectors at each level of the octree and to
5

compute the respective M2L operators. In the following we refer to such sets of unique transfer vectors as
interaction lists T ⊂ Z3 .
If we
√ consider asymptotically smooth kernel functions the near-field is limited to transfer vectors satisfying
|t| ≤ 3; it leads to 33 = 27 near-field interactions (see [5]). In a multi-level scheme, these 27 near-field
interactions contain all 63 = 216 near- and far-field
interactions of its eight child-cells. Far-field interactions
√
are given by transfer-vectors that satisfy |t| > 3. This leads to a maximum of 63 − 32 = 189 interactions per
cell and we end up with the usual overall complexity of O(N ) of fast multipole methods for asymptotically
smooth kernel functions. The union of all possible far-field interactions of the eight child cells gives 73 − 33 =
316 interactions. That is also the largest possible number of different M2L operators have to be computed per
octree level. Most asymptotically smooth kernel functions happen also to be homogeneous K(αr) = αn K(r)
of degree n. In other words, if we scale the distance r = |x − y| between source and target by a factor of α the
resulting potential is scaled by αn , where n is a constant that depends on the kernel function. The advantage
of homogeneous kernel functions is that the M2L operators need to be precomputed only at one level and can
simply be scaled and used on other levels. This affects the precomputation time and the required memory.
If we consider oscillatory kernel functions we need to distinguish between the low- and high-frequency
regime [12]. The admissibility criteria in the low-frequency regime are the same as those for asymptotically
smooth kernel functions. However, in the high-frequency regime the threshold distance between near-field
and far-field is O(kw2 ); nonetheless, as shown in [12] we end up with the usual complexity of O(N log N )
of fast multipole methods for oscillatory kernel functions. It depends on the wavenumber k (a property
of the kernel function). Thus, the size of near- and far-field is not a priori known as it is in the case of
asymptotically smooth kernel functions. Table 1 summarizes the number of near and far-field interactions
to be computed depending on different kernel functions.
Table 1: Number of near- and far-field interactions depends on the kernel function
type of kernel function
smooth
smooth and homogeneous
oscillatory

4

cells in near-field

cells in far-field

≤ 27 per leaf
≤ 27 per leaf
depends on k

≤ 316 per level
≤ 316 for all levels
depends on k

Optimizing the M2L operators

In all fast multipole methods the M2L operator adds the largest contribution to the overall computational
cost: for bbFMM it grows like O(N ) and for dFMM like O(N log N ). In this section, we first briefly recall
the optimization that was used up to now and suggest an improvement. Then, we present a new set of
optimizations that exploit the symmetries in the arrangement of the M2L operators.

4.1

Single low-rank approximation (SA)

In the following we explain the basics of the optimization used in [5, 12] and we refer to it as the SA variant
hereafter. The idea is based on the fact that all M2L operators Kt with t ∈ T can be assembled as a single
big matrix in two ways: either as a row of matrices K(row) = [K1 , . . . , Kt , . . . , K|T | ] or as a column of matrices
K(col) = [K1 ; . . . ; Kt ; . . . ; K|T | ] of M2L operators. The cardinality |T | gives the number of transfer vectors in
the interaction list T . Next, both big matrices are compressed using truncated singular value decompositions
(SVD) of accuracy ε as
K(row) ∼ UΣV∗ and K(col) ∼ AΓB∗
(17)

with the unitary matrices U, B of size ℓ × r and V, A of size |T |ℓ3 × r and the r singular values in Σ, Γ. With
a few algebraic transformations each M2L operator can be expressed as
Kt ∼ UCt B∗

where Ct = U∗ Kt B of size r × r is computed as
6

Ct = ΣV∗t B or

Ct = U∗ At Γ.

(18)

The advantage of this representation is that the cost of applying the M2L operator gets reduced from
applying a matrix of size ℓ3 × ℓ3 to a matrix of only r × r. Moreover, less memory is required. However, the
precomputation time grows cubically with the accuracy of the method due to the complexity of the SVD.
In [12] the SVD has been substituted by the adaptive cross approximation (ACA) followed by a truncated
SVD [1]. The precomputation time has been cut down drastically due to the linear complexity of the ACA.
4.1.1

SA with recompression (SArcmp)

If we use the SA variant the achieved low-rank r is the same for all M2L operators given by Ct . To a large
extent, r is determined by the greatest individual low-rank of the M2L operators. This means that most of
the matrices Ct of size r × r have effectively a smaller low-rank rt ≤ r. We exploit this fact by individually
approximating them as
Ct ∼ Ūt V̄t∗

with Ūt , V̄t of size r × rt and the constraint rt < r/2.

(19)

Without the constraint the low-rank representation is less efficient than the original representation. The
effects of the recompression are studied in Sec. 5.2.

4.2

Individual low-rank approximation (IA)

As opposed to the SA approach, an individual low-rank approximation of the M2L operators as
Kt ∼ Ut Vt∗

with Ut Vt of size ℓ3 × rt

(20)

directly leads to the optimal low-rank representation of each of them. As in the previous section, the
approximation can be performed by either a truncated SVD or the ACA followed by a truncated SVD (note,
the rank rt in the Eqns. (19) and (20) might be similar but has not to be the same). All these variants (SA,
SArcmp and IA) still require the approximation and storage of all M2L operators. In terms of time and
memory however, it would be desirable to come up with a method that requires the approximation and the
storage of a subset of operators only. Let us present a set of optimizations that fulfill these two requests.
4.2.1

Symmetries and permutations

Here, we illustrate how the full set of M2L operators can be expressed by a subset only. The idea is based on
symmetries in the arrangement of M2L operators and exploits the uniform distribution of the interpolation
points. We start by presenting the idea using a model example. After generalizing this idea, we demonstrate
that M2L operators can be expressed as permutations of others.
Model example The target cell X in Fig. 1 has three interactions Yt with the transfer vectors t ∈
{(2, 1), (1, 2), (−2, 1)}. We choose the reference domain to be given by t1 ≥ t2 ≥ 0. The goal is to express the
M2L operators of all interactions via M2L operators of interactions that lie in the reference domain only. In
our example this is the interaction with the transfer vector (2, 1). The two transfer vectors (1, 2) and (−2, 1)
can be shown to be reflections along the lines given by t1 = 0 and t1 = t2 , respectively. We easily find that
any t ∈ Z2 can be expressed as a reflection (or a combination of reflections) of transfer vectors that satisfy
t1 ≥ t2 ≥ 0.
We claim that any reflection of a transfer vector corresponds to a permutation of the respective M2L
operator. Recall that the evaluation of K(x̄α , ȳβ ) gives the entry from row m(α) and column n(β) of the
M2L operator. K(2,1) is the only M2L operator whose transfer vector satisfies t1 ≥ t2 ≥ 0. The multi-indices
are constructed as presented in Fig. 1. As can be checked, the entry (K(2,1) )mn is not the same as (K(1,2) )mn
or (K(−2,1) )mn . However, if we use the permuted multi-indices from Fig. 2a for K(−2,1) or those from Fig. 2b
for K(1,2) they are the same. The logic behind this can be summarized as follows.

7

β = (β1 , β2 )

axial sym.

α = (α1 , α2 )

(1,3) (2,3) (3,3)

(1,2) (2,2) (3,2)

g.

.
m

(1,1) (2,1) (3,1)

Y(1,2)

(1,3) (2,3) (3,3)

(1,2) (2,2) (3,2)

a
di

sy

β2

Y(−2,1)

(1,1) (2,1) (3,1)

(1,3) (2,3) (3,3)

(1,2) (2,2) (3,2)

Y(2,1)

(1,1) (2,1) (3,1)

α2

β1

(1,3) (2,3) (3,3)

e2

e1

(1,2) (2,2) (3,2)

X
(1,1) (2,1) (3,1)

α1
Figure 1: Axial and diagonal symmetries of interactions. The interpolation points x̄α and ȳβ are indexed by
the multi-indices α and β, respectively (interpolation order ℓ = 3). The only transfer vector that satisfies
t1 ≥ t2 ≥ 0 is t = (2, 1). In that case, we claim that K(2,1) is the only M2L operator we need to compute.
The direction of the arrows indicates the growth of the respective multi-index component.
• If an axial symmetry is given by t1 = 0 as shown in Fig. 2a, we invert the corresponding component of
the multi-index as
α ← (ℓ − (α1 − 1), α2 ) and β ← (ℓ − (β1 − 1), β2 ).
(21)
• If the diagonal symmetry is given by t1 = t2 as shown in Fig. 2b, we swap the corresponding components
as
α ← (α2 , α1 ) and β ← (β2 , β1 ).
(22)
Sometimes it is necessary to combine axial and diagonal permutations. Take as example the transfer vector (−1, 2): we need to flip it along t1 = 0 and then along t1 = t2 to get (2, 1). Note that reflections are
non commutative, i.e., the order of their application matters. This is also true for permutations of the M2L
operators.
Generalization Let us extend the above concept to the three dimensional case. We start by introducing
three axial and three diagonal symmetries in Z3 .
• Axial symmetry planes are given by t1 = 0, t2 = 0 and t3 = 0 (see Fig. 3). Each of the three planes
divides Z3 in two parts, i.e., the negative part ti < 0 and the positive part ti ≥ 0. By combining all
three planes Z3 is divided into octants. In the following we use Z3+ , i.e., the octant with t1 , t2 , t3 ≥ 0
as reference octant.
• Diagonal symmetry planes are given by t1 = t2 , t1 = t3 and t2 = t3 (see Fig. 4). In Z3 there are six
diagonal symmetries; however, we restrict ourselves to the symmetries affecting the reference octant
Z3+ .
By combining the three diagonal symmetries and the three axial symmetries we obtain the cone shown in
Fig. 5. We refer to it as Z3sym ; it is given by

Z3sym = Z3sym ⊂ Z3 : t1 ≥ t2 ≥ t3 ≥ 0 with t ∈ Z3 .
8

(23)

α = (ℓ − (α1 − 1), α2 )

(3,3) (2,3) (1,3)

Y(−2,1)

(3,2) (2,2) (1,2)

β2

β = (β2 , β1 )

axial sym.

β = (ℓ − (β1 − 1), β2 )

α = (α2 , α1 )

(3,1) (3,2) (3,3)

β1

(2,1) (2,2) (2,3)

(1,1) (1,2) (1,3)

β2

a
di

(3,1) (2,1) (1,1)

β1

Y(1,2)

(3,3) (2,3) (1,3)

α2

α1

g.

sy

.
m

(3,1) (3,2) (3,3)

(3,2) (2,2) (1,2)

(2,1) (2,2) (2,3)

(3,1) (2,1) (1,1)

(1,1) (1,2) (1,3)

X

X
α1

α2

(a) Invert the component α1 and β1 due to the axial symmetry
t1 = 0.

(b) Swap the components α1 ↔ α2 and β1 ↔ β2 due
to the diagonal symmetry t1 = t2 .

Figure 2: The direction of the arrows indicates the growth of the respective multi-index component such
that the M2L operators K(−2,1) and K(1,2) become the same as K(2,1) . In other words, the mapping from
the arrows in Fig. 1 to the arrows here is analog to the mapping of the multi-indices and corresponds to the
permutations of K(2,1) such that it coincides with K(−2,1) and K(1,2) .
By its means we can identify the subset of transfer vectors Tsym ⊂ T ⊂ Z3 as
Tsym = T ∩ Z3sym

(24)

such that all others T \Tsym can be expressed as reflections of this fundamental set. Next, we claim that
these symmetries are also useful for M2L operators.
Permutation matrices Any reflection of a transfer vector along a symmetry plane determines the permutation of its associated M2L operator as
Kt = Pt Kp(t) P⊤
t .

(25)

The permutation matrix Pt depends on the transfer vector t ∈ T . We also need the surjective mapping
p : T → Tsym ; it associates every transfer vector in T to exactly one in Tsym . The left application of Pt ,
essentially, corresponds to the permutation of α and its right application to the permutation of β, affecting
rows (respectively columns) of the original matrix Kp(t) . Note, the permutation matrices Pt depend only on
the transfer vector t. How do we construct them? For some t we introduce axial and diagonal permutations
πtA and πtD that read as follows.
• Axial symmetries: multi-index permutations are computed as
(
πtA (α1 , α2 , α3 ) = (ᾱ1 , ᾱ2 , ᾱ3 )

with

ᾱi =

αi
if ti ≥ 0,
ℓ − (αi − 1) else.

(26)

There exist 8 different possibilities that correspond to the octants presented in Fig. 3. Note, πtA (α) = α
is only true for transfer vectors with t1 , t2 , t3 ≥ 0.
9

e3

e2
e1

Figure 3: Three axial symmetry planes split Z3 in octants. The reference octant is given by t1 , t2 , t3 ≥ 0.

e3

e3

e2
e1

e3

e2
e1

(a) t1 = t2

e2
e1

(b) t1 = t3

(c) t2 = t3

Figure 4: Three diagonal symmetry planes in the reference octant.

e3

e2
e1

Figure 5: The final cone Z3sym (t1 ≥ t2 ≥ t3 ≥ 0) is obtained by combining axial and diagonal symmetries.

10

• Diagonal symmetries: multi-index permutations are computed as
πtD (α1 , α2 , α3 ) = (αi , αj , αk )

such that |ti | ≥ |tj | ≥ |tk |.

(27)

There exist 6 different possibilities that correspond to the 6 different cones if we consider Fig. 5. Note
again, πtD (α) = α is only true for transfer vectors with t1 ≥ t2 ≥ t3 ≥ 0.
Given these multi-index permutations and the mapping functions m(α) and n(β) we can define a permutation
matrix Pt of size ℓ3 × ℓ3 . Its entries are 0 except in column j the entry i = m(πt (m−1 (j))) is 1. Let us go
through the computation of this index: first, we compute the multi-index α = m−1 (j), then, we permute the
multi-index ᾱ = πt (α) and last, we compute the row-index i = m(ᾱ). Permutation matrices may be written
as

Pt = em(πt (m−1 (0))) , em(πt (m−1 (1))) , . . . , em(πt (m−1 ((ℓ−1)3 ))) ,
(28)

where ej denotes a column unit vector of length ℓ3 with 1 in the jth position and 0 elsewhere. Permutation
−1
matrices are orthogonal Pt P⊤
= P⊤
t = I, hence, the inverse exists and can be written as Pt
t . Note that the
A
combination of permutations is non commutative. Given these permutations πt and πtD we setup PA
t and
PD
and
construct
the
permutation
matrix
as
t
A
Pt = PD
t Pt .

(29)

The permutation for the multi-index β is the same.
4.2.2

IA with symmetries (IAsym)

By exploiting the above introduced symmetries we end up with an optimization we refer to as the IAsym
variant. We individually approximate and store only M2L operators with t ∈ Tsym and express all others via
permutations as shown Eqn. (25). The IAsym variant for an arbitrary transfer vector t ∈ T consists of the
following three steps.
1. Permute multipole expansions
w t = P⊤
t w

(30)

ft = Kp(t) wt

(31)

f = Pt ft

(32)

2. Compute permuted local expansions
3. Un-permute local expansions
Note that the permutation matrix is not applied to the actual M2L operator (remains unchanged as can be
seen in step 2). Its application is implemented as a reordering of vector entries (step 1 and 3). Depending
on whether the M2L operator exist in its full-rank or in its low-rank representation (see Eqn. (20)) the
application corresponds to one or two matrix-vector products. In the following we introduce a blocking
scheme that leads to a faster execution on a computer.
4.2.3

IAsym with blocking (IAblk)

We know from Sec. 4.2.1 that based on the consideration of symmetries and permutations, many interactions share the same M2L operators. This paves the road for blocking schemes. Essentially, the idea is to
substitute many matrix-vector products by a few matrix-matrix products. Blocking schemes do not change
the overall complexity of the algorithm, but they allow for the use of highly optimized matrix-matrix product implementations. Such achieve much higher peak performances than optimized matrix-vector product
implementations due to better cache reuse and less memory traffic [4, 9].
In our concrete case, we use the blocking scheme to block multipole and local expansions. Instead
of permuting them and applying the M2L operators individually (matrix-vector products), we assemble
11

Algorithm 2 Blocking scheme with |Tsym | matrix-matrix products
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:

function BlockedM2L(target cell X and all far-field interactions IY )
allocate Fp and Wp for p = 1, . . . , |T |sym
retrieve f from X
set all cp = 0
for all source cells Y in IY do
retrieve w from Y and compute t from cell-pair (X, Y )
column cp(t) of Wp(t) gets P⊤
⊲ Permute multipole expansions
t w
increment cp(t)
end for
for all {Kp } do
Fp ← Kp Wp
⊲ Compute permuted local expansions
end for
set all cp = 0
for all source cells Y in IY do
compute t from cell-pair (X, Y )
retrieve ft from column cp(t) of Fp(t)
increment cp(t)
f ← f + Pt ft
⊲ Permute permuted local expansions
end for
end function

those that share the same M2L operator as distinct matrices to whom we apply the M2L operators then
(matrix-matrix products). Algorithm 2 explains this in details. We need the matrices Wp and Fp of size
ℓ2 × np for p = 1, . . . , |Tsym |. Their columns store the permuted multipole and the resulting (also permuted)
local expansions. The values for np indicate how many interactions in T share the same M2L operator of
interactions in Tsym , in other words, n1 + · · · + np + · · · + n|Tsym| = |T | is true. In the case of bbFMM this is a
priori known, since the full interaction list is a priori given (see Sec. 3.2). That is not the case for dFMM and
the values for np have to be determined during a precomputation step. We also need counters cp to indicate
the position of the currently processed expansions in Wp and Fp . As opposed to IAsym, here we split up
the single loop over all interactions into three loops. In the first one, we assemble the set of matrices Wp .
At the end cp ≤ np is true for all p. In the second loop, we perform at most |Tsym | matrix-matrix products.
And in the last loop, we increment the local expansion with the expansions from all Fp .
Blocking along multiple target cells Algorithm 2 proposes to use the blocking scheme for all interactions of only one target cell. In the worst case no M2L operator is shared and the algorithm coincides with
IAsym. Moreover, the size of the matrices Wp and Fp might vary to a large extent. That is why we pursued
the blocking idea further to come up with a more efficient scheme. Instead of using individual np we choose
it to be a constant n for all p = 1, . . . , |Tsym |. Then we keep on blocking expansions using interactions lists
of multiple (as opposed to one) target cells. Once cp = n is true for some p, we apply the M2L operator as
Fp = Kp Wp where Wp , Fp are both of size ℓ3 × n. In our numerical studies we use this blocking scheme with
n = 128.

5

Numerical results

In the previous sections we introduced various optimizations of the M2L operators for bbFMM and dFMM.
As representative kernel functions we use
the Laplace kernel K(x, y) =

1
4π|x − y|

and the Helmholtz kernel K(x, y) =

12

eık|x−y|
.
4π|x − y|

(33)

In the numerical studies, hereafter, we use the same parameter setting as are used in the respective publications and for dFMM we use the wavenumber k = 1. All computations are performed on a single CPU of
a Intel Core i7–2760QM CPU @ 2.40GHz × 8 with 8GB shared memory. We used the compiler gcc 4.6.3
with the flags “-O2 -ffast-math”. The source files (C++) of the program we used can be downloaded via
http://github.com/burgerbua/dfmm; they are distributed under the “BSD 2-Clause” license.
M2L optimizations We show and analyze results for the following eight variants:
1. NA: the full interaction list T is represented by full-rank (not approximated) M2L operators
2. NAsym: the reduced interaction list Tsym is represented by full-rank (not approximated) M2L operators
3. NAblk: same as NAsym but with additional blocking of multipole and local expansions
4. SA: the variant presented in [5] and briefly sketched in Sec. 4.1
5. SArcmp: same as SA but with additional recompression of all Ct
6. IA: same as NA but with low-rank (individually approximated) M2L operators
7. IAsym: same as NAsym but with low-rank (individually approximated) M2L operators
8. IAblk: same as NAblk but with low-rank (individually approximated) M2L operators (see Alg. 2)
Moreover, we study two different low-rank approximation schemes for the SA and IA variants: on one
hand we use a truncated singular value decomposition (SVD) and on the other hand the adaptive cross
approximation (ACA) followed by a truncated SVD [1].
Example geometries We use the three benchmark examples shown in Fig. 6 and described in the listing
below. The depth of the octree is chosen such that near- and far-field are balanced, i.e., we want the fastest
possible matrix-vector product.
1. The sphere from Fig. 6a is contained in the bounding-box 64 × 64 × 64; 168 931 particles are randomly
scattered on its surface. The octree has 6 levels.
2. The oblate sphere from Fig. 6b is contained in the bounding-box 6.4 × 64 × 64; 125 931 particles are
randomly scattered on its surface. The octree has 6 levels.
3. The prolate sphere from Fig. 6c is contained in the bounding-box 6.4 × 6.4 × 64; 119 698 particles are
randomly scattered on its surface. The octree has 7 levels.
In approximative sense, the sphere is a three-dimensional, the oblate sphere a two-dimensional and the
prolate sphere a one-dimensional object in R3 . We choose these three geometries to study their influence on
the performance on the dFMM. Table 2 shows the size of the full and the reduced interaction lists (|T | and
|Tsym |) per level (lf stands for low-frequency, hf for high-frequency regime) for all three geometries. The size
of the interaction lists clearly grows with the dimensionality of the geometry. We report on the impact of
this behavior later.

3(hf)
|T |
|Tsym |

668
20

sphere
4(hf) 5(hf)
18710
518

2666
93

6(hf)

3(hf)

418
21

60
6

oblate sphere
4(hf) 5(hf)
2336
203

1502
89

6(hf)

4(hf)

400
21

214
35

prolate sphere
5(hf) 6(hf)
738
61

382
21

7(lf)
424
22

Table 2: Size of interactions lists per level for dFMM (hf and lf stands for high- and low-frequency regime,
respectively)

13

❨❳

(a) sphere

❨

❨❳
❩

❩

(b) oblate spheroid

❩

❳

(c) prolate spheroid

Figure 6: The example geometries are centered at (0, 0, 0)

5.1

Accuracy of the method

Both, the bbFMM and the dFMM, have two approximations: 1) the interpolation of the kernel functions
determined by interpolation order ℓ, and 2) the subsequent low-rank approximation of the M2L operators
determined by the target accuracy ε. The final relative error is a result of both approximations. We compute
it as
P
¯ 2 1/2
i∈M |fi − fi |
P
εL 2 =
(34)
i∈M |fi |

where M is the number of particles x in an arbitrary reference cluster at the leaf level; f and f¯ are the exact
and approximate results, respectively. In Fig. 7 we compare achieved accuracies for the bbFMM and the
dFMM with the IAblk variant (other variants produce identical results). Both plots show the behavior of the
relative error εL2 depending on the interpolation order ℓ and the target accuracy ε. Evident is the matching
result between the left and right figure. All curves show an initial plateau and then, after a sharp knee a
slope of roughly 1. The knee occurs approximately at (ℓ, ε) = (Acc, 10−Acc ). In the rest of the paper we use
this convention to describe the accuracy Acc of bbFMM and dFMM. The low-rank approximations for the
computations whose accuracies are shown in Fig. 7 were conducted with the ACA followed by a truncated
SVD. By just using a truncated SVD we obtain identical results.

5.2

Reducing the cost with the SArcmp variant

The cost of applying an approximate M2L operator mainly depends on its rank k. In Tab. 3 we compare
the average rank of M2L operators for the SA and the IA variants at all levels that have expansions. Let us
explain how we computed the average rank for SA. Recall, when we use that variant, all M2L operators from
one interaction list posses the same rank; the bbFMM and the dFMM in the low-frequency regime have one
and the dFMM in the high-frequency regime has potentially multiple (directional) interaction lists per level.
Thus, the average rank per level is the average of all ranks used in all interactions at that level.
The application of one M2L operator from SA and IA requires O(k 2 ), respectively O(2kℓ3 ) operations.
Note, the ranks k of the M2L operators are different; for a given accuracy Acc they are normally significantly
lower for IA than for SA. This can be seen in Tab. 3. The large ranks at level 7 (first level in the low-frequency
regime) of SA are noteworthy. There, the lower bound for the separation criterion is given by the usual lowfrequency criterion saying that non touching cells are admissible. Hence, the smallest possible transfer
14

Relative error εL2

10−1

10−3

10−5

10−1

ℓ=3
ℓ=4
ℓ=5
ℓ=6
ℓ=7
ℓ=8
ℓ=9

10−3

10−5

10−7

ℓ=3
ℓ=4
ℓ=5
ℓ=6
ℓ=7
ℓ=8
ℓ=9

10−7

10−9
10−12 10−10

−8

−6

−4

10
10
10
10
Target accuracy εACA

−2

10

0

10−9
10−12 10−10

(a) bbFMM for the smooth kernel function

10−8 10−6 10−4 10−2
Target accuracy εACA

(b) dFMM for the oscillatory kernel function

Figure 7: Accuracies for the prolate sphere from Fig. 6c. The target accuracy εACA refers to the accuracy
of the approximate M2L operators (see Eqn. (20)). Here, we used the ACA followed by a truncated SVD.

Acc

4(hf)

3
5
7
9

9.8
21.7
38.2
57.8

SA
5(hf) 6(hf)
12.3
30.8
58.6
96.5

12.8
39.2
80.0
138.7

7(lf)

4(hf)

19
71
163
296

5.4
11.3
18.9
28.6

IA
5(hf) 6(hf)
5.7
12.3
21.5
33.2

5.7
13.5
24.7
39.7

7(lf)
5.0
12.6
24.1
40.1

Table 3: Comparison of average ranks k for the SA and IA variants of the dFMM (prolate sphere)

15

100

vectors have a length of mint∈T |t| = 2. The slowly decaying singular values of associated M2L operators are
responsible for the large ranks. On the other hand, the upper bound for the separation criterion coincides
with the lower bound of level 6 (parent level), which is in the high-frequency regime. Hence, the largest
possible transfer vectors have a length of maxt∈T |t| ∼ 4k. M2L operators whose transfer vectors are in that
range have much faster decaying singular values. This fact explains the efficiency of SArcmp (individual
recompression of each M2L operator).
In Tab. 4 we analyze the cost savings of SArcmp compared to SA. The left values in each column
multiplied by 106 give the overall number of floating point operations per level for SA. The right values (in
brackets) indicate the respective cost ratio of SArcmp to SA. The recompression reduces the cost remarkably
(see also Fig. 8). At the low-frequency level 7, SArcmp reduces the cost by more than a factor of 2. This is
almost twice as much as in high frequency levels. For the impact on timing results we refer to Sec. 5.4.

Acc
3
5
7
9

4(hf)
0.8
4.9
19.3
55.0

cost(SA)/106
5(hf)

(0.98)
(0.97)
(0.97)
(0.97)

20.3
161.4
687.8
2138.2

(cost(SArcmp)/ cost(SA))
6(hf)

(0.93)
(0.89)
(0.88)
(0.87)

31.5
339.3
1590.2
5237.6

(0.93)
(0.86)
(0.83)
(0.83)

7(lf)

204.3
2889.3
15420.0
51505.7

(0.62)
(0.47)
(0.40)
(0.38)

Table 4: Comparison of cost in terms of floating point operations for SA and SArcmp (prolate sphere)
In Fig. 8 we compare the cost of SA, SArcmp and IA for bbFMM and dFMM for the prolate sphere and
accuracy Acc = 9. The bbFMM in Fig. 8a has expansions at the levels 2 − 7. The reason for the jump
between level 4 and 5 is because the levels 2 − 4 have a maximum of |T | = 16 and the levels 5 − 7 a maximum
of |T | = 316 (common maximum size for bbFMM) possible M2L operators. The jump in the case of dFMM
in Fig. 8b at level 5 can be explained in the same way: if we look at Tab. 2 we see that |T | is about twice as
large as it is at the other levels. The levels 4, 5, 6 of dFMM are in the high-frequency regime. There, the cost
of IA is approximately 12, 5, 3 times greater than the cost of SA. However, at level 7 (low-frequency regime)
of dFMM and at the levels 5 − 7 of bbFMM the cost of IA is only about 2/3 compared to SA. The reason is
the size of the interaction lists |T |. The larger they become the larger the span between the smallest and the
largest rank and that favors individual approximations. The SArcmp is computationally the least expensive
variant. Similar results are obtained for all other accuracies.

5.3

Speeding up the precomputation with IA variants

The bottleneck of SA and SArcmp is the relatively large precomputation time. This is a minor issue in the
case of homogeneous kernel functions. In that case the approximated M2L operators can be stored on the disk
and loaded if needed for further computations. However, if we deal with non-homogeneous kernel functions,
such as the Helmholtz kernel, IAsym is the way to go. In Tab. 5 we compare the precomputation time of
SA, IA and IAsym (we do not report on SArcmp; due to the additional recompression its precomputation
time is higher than the one for SA). For the low-rank approximation we use a truncated SVD or the ACA
followed by a truncated SVD. In both cases we end up with the same rank. We get remarkable speedups in
the precomputation. Let us look at an extreme example: for the sphere and an accuracy Acc = 7 we have
tSVD = 4740.9 s for SA versus tACA = 1.8 s for IAsym. This corresponds to a speedup greater than 2600.
For the oblate and prolate sphere we obtain similar results.

5.4

Comparison of all M2L variants

In the previous sections we have revealed the impact of SArcmp and IAsym on the computational cost and
the precomputation time of the M2L operators, respectively. In this section we compare all variants and
focus on their impact on memory requirement and running time (of applying M2L operators during the
16

floating point operations

1011

1011
SA
SArcmp
IA

SA
SArcmp
IA

1010

1010

109

109

108

108

107

2

3

4

5

6

octree level

7

107

4

5

6

7

octree level

(a) bbFMM

(b) dFMM

Figure 8: Comparison of the M2L cost (floating point operations) growth per level for the SA, SArcmp and
IA variants of bbFMM and dFMM (Acc = 9 and prolate sphere)
actual matrix-vector product). Tables 6–8 (respectively, 9–11) present running times of the M2L operators
for bbFMM (respectively, dFMM) and all three geometries. The upper set of rows reports on results obtained
with a BLAS and LAPACK implementation (libblas 1.2.20110419-2, which is sub-optimal in our case). The
lower set shows the times obtained with the Intel MKL [9] (Intel MKL 10.3.11.339 (intel64)), which proved
faster for our purpose.
Impact of using symmetries on the memory requirement Missing times in the tables indicate that
the required memory exceeded the available memory of 8 GB. Clearly, the NA variant (no approximation, no
use of symmetries) requires the most memory. On the other hand, IAsym and IAblk are the most memory
friendly variants. Both require only memory for storing the low-rank representations of M2L operators from
Tsym (see Tab. 2 for a comparison of T and Tsym ).
Impact of the blocking scheme on runtime Bold numbers indicate the fastest variants. Two results
attract our attention. 1) If we look at the times in the Tabs. 6–11 we notice that in four cases (bbFMM
with the sphere, the oblate sphere and the prolate sphere and dFMM with the prolate sphere) IAblk, and
in two cases (dFMM with the sphere and the oblate sphere) SArcmp is fastest. To be more general, IAblk
wins at levels having non-directional expansions (all levels of bbFMM and low-frequency levels of dFMM)
and SArcmp at levels having directional expansions (high-frequency levels of dFMM). Why is that? The
reason follows from the cost studies in Sec. 5.2. Let us take for example Acc = 5. Recall, we measure the
computational cost based on the size of the approximate M2L operators, ie., O(k 2 ) for SA and O(2kℓ2 ) for
IA. The respective ranks k are given in Tab. 3. The ratio of these costs for SA to IA is 0.46 at the highfrequency level 6 and 1.60 at the low-frequency level 7. As a matter of fact, at high-frequency levels wins SA
and at low-frequency levels IA. Even savings of 47% at the low-frequency level due to the recompression of
SArcmp (see Tab. 4) are not sufficient to outweigh the advantage of the blocking scheme of IAblk. If there
is no low-frequency level, such as for the sphere in Tab. 9 and the oblate sphere in Tab. 10, the SArcmp
outperforms all other variants. For example, if we would repeat the computations for the prolate sphere
with an octree of depth 6 (no low-frequency level) the resulting timing patterns would follow those from

17

Acc

SA
tSVD [s] tACA [s]

tSVD [s]

IA
tACA [s]

IAsym
tSVD [s] tACA [s]

sphere
3
4
5
6
7

7.0
75.0
317.2
1336.4
4740.9

4.6
20.2
69.1
197.0
435.9

6.2
37.0
188.8
790.5
-

1.9
4.9
12.4
28.6
-

0.2
1.1
5.6
22.9
84.0

0.1
0.2
0.4
0.9
1.8

0.4
1.1
2.7
6.4
13.3
29.0
-

0.1
0.5
2.8
11.1
40.9
129.3
369.2

0.0
0.1
0.2
0.5
0.9
2.0
3.9

0.3
0.7
2.1
5.0
11.6
25.0
50.2

0.1
0.4
2.5
15.7
66.0
322.3
807.4

0.0
0.1
0.1
0.4
0.8
1.7
3.7

oblate sphere
3
4
5
6
7
8
9

1.4
13.4
59.7
299.4
1021.2
-

0.9
3.8
13.5
39.4
97.6
217.7
444.6

1.2
7.1
37.0
150.1
551.0
1751.4
prolate sphere

3
4
5
6
7
8
9

0.6
7.5
64.7
358.2
1374.8
-

0.7
4.1
18.3
73,3
204.3
549.3
1273.2

0.8
5.3
31.0
199.4
808.9
-

Table 5: Precomputation times (SVD versus ACA) for the SA, IA and IAsym variants. Missing numbers
mean that the available memory of 8 GB has not been sufficient for the respective computation.
the sphere and the oblate sphere (the overall application time would increase too, since the tree-depth is
based on our choice of balancing near- and far-field, i.e., the shortest overall application time). 2) Evident
is the wide margin in the speedup of variants that use blocking and those which do not. If we use the MKL
(as opposed to libblas) for NA, NAsym, SA, SArcmp, IA and IAsym we end up with 1.5 − 2 times faster
application times. However, if we use the MKL for NAblk and IAblk we achieve 3 − 4 times faster times.
Even though these speedups are greatest with the MKL library, it highlights the benefits of the blocking
scheme presented in Alg. 2.
Varying growth of application times In Fig. 9 we visualize the different growth rates of M2L application
times for the bbFMM with increasing accuracies Acc. We are interested in the growth rates due to algorithmic
changes. That is why we only study those variants that do not use blocking. Since no approximation is
involved the times for NAsym grows the fastest. The times for SA grow slower but still faster than those for
IAsym. SArcmp features the slowest growth, it is the optimal variant in terms of computational cost (see
Tab. 4).

6

Conclusion

The fast multipole method based on Chebyshev interpolation, first presented in [5] for smooth kernel functions
(bbFMM) and extended in [12] to oscillatory ones (dFMM), is a convenient-to-use algorithm due to its easy

18

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

0.4
1.3
3.4
9.0
20.2
40.7
74.0

0.4
1.4
3.6
8.7
18.4
35.9
68.6

0.3
0.7
1.8
4.2
8.6
16.2
28.4

0.4
0.8
2.0
4.0
8.2
16.2
32.9

0.2
0.4
1.0
1.7
3.5
6.0
9.9

libblas 1.2.20110419-2
3
4
5
6
7
8
9

0.7
3.6
14.2
42.8
102.7
229.3
-

0.8
3.8
12.8
38.2
101.7
234.2
484.5

0.4
1.8
5.9
16.5
40.7
89.7
180.0

0.5
1.4
4.6
11.5
25.5
47.8
83.8

0.3
0.9
2.2
4.8
9.6
18.3
30.0

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7
8
9

0.3
1.8
9.4
28.2
72.6
127.6
-

0.3
1.2
6.3
14.1
57.7
117.0
260.5

0.2
0.6
1.9
4.8
12.0
25.6
50.8

0.2
0.6
2.0
6.9
19.3
34.3
60.4

0.4
0.7
1.4
2.6
5.8
12.0
20.6

0.4
0.7
2.0
5.0
12.7
24.6
44.6

Table 6: M2L timings for the bbFMM (sphere). In this table and below as well, bold numbers correspond to
the smallest entry in a row. In some cases, two columns use bold font when the running times are sufficiently
close that the difference is not significant.

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

0.2
0.6
1.6
4.0
9.0
18.2
32.4

0.2
0.6
1.6
3.9
9.0
16.5
30.4

0.1
0.3
0.8
1.9
3.8
7.2
12.8

0.2
0.4
0.9
1.8
4.4
7.2
14.2

0.1
0.2
0.4
0.8
1.6
2.7
4.4

libblas 1.2.20110419-2
3
4
5
6
7
8
9

0.3
1.6
6.4
19.1
46.8
102.8
-

0.4
1.6
5.8
16.7
44.8
101.9
212.2

0.2
0.8
2.6
7.4
18.3
40.1
81.1

0.2
0.8
1.9
5.4
11.2
21.3
37.0

0.1
0.4
0.9
2.1
4.5
8.4
13.5

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7
8
9

0.1
0.7
4.5
13.1
33.7
57.9
117.4

0.1
0.5
2.0
6.7
27.4
52.7
118.2

0.1
0.3
0.9
2.2
5.4
11.5
22.9

0.1
0.3
1.0
3.2
8.2
14.6
27.3

0.2
0.3
0.7
1.2
2.7
5.2
9.1

0.2
0.3
0.9
2.6
5.7
10.9
20.0

Table 7: M2L timings for bbFMM (oblate sphere)

19

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

0.1
0.4
1.2
2.6
7.2
11.6
21.1

0.1
0.4
1.0
2.4
5.2
10.8
19.9

0.1
0.2
0.6
1.4
2.9
4.9
8.7

0.1
0.2
0.6
1.1
2.4
5.5
11.2

0.1
0.1
0.3
0.5
1.0
1.7
2.8

libblas 1.2.20110419-2
3
4
5
6
7
8
9

0.2
1.0
4.0
12.0
29.3
68.4
131.7

0.2
1.1
4.0
11.8
31.4
71.4
137.3

0.1
0.5
1.9
4.7
11.7
25.5
50.7

0.1
0.4
1.2
3.4
6.9
13.1
22.3

0.1
0.3
0.6
1.3
2.8
5.2
8.6

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7
8
9

0.1
0.4
2.9
8.2
20.5
35.7
73.3

0.1
0.3
1.3
4.2
16.3
34.1
73.8

0.1
0.2
0.6
1.4
3.4
7.2
14.4

0.1
0.2
0.6
2.1
5.9
9.1
17.6

0.1
0.2
0.4
0.7
1.7
3.7
5.9

0.1
0.2
0.6
1.5
3.7
7.1
12.8

Table 8: M2L timings for bbFMM (prolate sphere)

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

3.2
8.0
19.8
43.0
-

2.9
8.1
19.4
42.8
86.5

3.5
8.3
19.7
40.4
78.1

2.4
5.9
14.0
29.7
60.8

2.2
4.8
10.0
17.6
30.8

libblas 1.2.20110419-2
3
4
5
6
7

6.3
25.7
113.0
-

5.9
25.1
89.5
275.9
-

4.8
20.1
71.4
202.2
-

2.0
4.3
7.4
14.5
-

2.0
3.9
7.1
11.6
-

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7

4.8
21.5
109.1
-

3.9
17.1
61.5
162.2
-

2.3
7.5
22.8
58.8
-

1.6
3.3
6.3
11.6
-

1.7
3.4
6.3
9.8
-

3.2
6.7
10.6
35.1
-

Table 9: M2L timings for dFMM (sphere); high-frequency leaf level

20

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

0.9
2.7
6.9
15.9
32.1
59.2
-

0.9
2.6
6.8
15.9
32.3
59.6
105.1

1.1
2.7
6.9
14.8
28.9
52.7
90.1

0.8
1.8
4.5
10.3
21.8
42.0
79.6

0.6
1.3
3.1
5.9
11.0
18.8
30.9

libblas 1.2.20110419-2
3
4
5
6
7
8
9

1.8
8.9
31.4
91.8
-

1.8
8.6
30.6
95.9
228.4
-

1.6
6.8
24.8
71.6
176.6
-

0.6
1.3
2.8
5.5
9.6
16.4
25.5

0.6
1.1
2.3
4.1
7.3
11.2
21.9

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7
8
9

1.4
7.2
25.8
62.5
-

1.0
5.0
20.7
56.7
141.3
-

0.7
2.4
7.9
20.4
48.8
-

0.4
0.9
2.2
4.3
7.8
13.6
21.0

0.4
1.0
1.9
3.4
5.5
9.5
14.7

0.7
2.0
5.5
12.6
24.7
45.3
-

Table 10: M2L timings for dFMM (oblate sphere); high-frequency leaf level

Acc

NA

NAsym

NAblk

SA

SArcmp

IA

IAsym

IAblk

0.3
0.9
2.4
6.0
12.1
24.3
43.5

0.3
0.9
2.4
5.6
11.8
23.3
43.5

0.3
1.0
2.5
5.7
11.1
21.2
37.2

0.2
0.6
1.5
3.4
7.5
15.8
32.1

0.2
0.4
1.0
2.0
3.8
7.0
11.6

libblas 1.2.20110419-2
3
4
5
6
7
8
9

0.6
3.2
11.9
32.5
80.3
-

0.6
2.8
10.0
30.8
77.0
-

0.5
2.3
8.7
25.1
61.1
-

0.3
1.0
2.7
6.5
13.2
24.8
47.5

0.2
0.5
1.2
2.6
5.1
9.1
18.6

Intel MKL 10.3.11.339 (intel64)
3
4
5
6
7
8
9

0.3
2.1
8.7
20.2
48.4
-

0.3
1.4
5.8
17.9
46.5
-

0.2
0.7
2.5
6.7
16.3
-

0.2
0.5
2.0
5.1
10.3
16.6
31.6

0.2
0.4
1.0
2.0
4.0
7.2
14.8

0.2
0.6
2.0
4.6
9.6
18.3
33.7

Table 11: M2L timings for dFMM (prolate sphere); low-frequency leaf level

21

NAsym (s)

M2L application time [s]

NAsym (ps)

102

SA (s)
SA (ps)
SArcmp (s)
SArcmp (ps)

101

IAsym (s)
IAsym (ps)

100

10−1
3

4

5

6

7

8

9

Accuracy Acc
Figure 9: Running times versus accuracy for NAsym, SA, SArcmp and IAsym for bbFMM taken from the
Tabs. 6, 7 and 8; (s) stands for sphere and (ps) for prolate sphere
implementation and its kernel function independence. In this paper we investigated algorithms to reduce the
running time of the M2L operator. We proposed several optimizations and studied their respective strengths
and weaknesses.
On one hand we proposed SArcmp, which uses an individual recompression of the suboptimally approximated M2L operators obtained via SA (the variant presented in [5]). We have shown that this new variant
reduces the computational cost noticeably. In some settings it even provides the fastest M2L application
times. On the other hand we also proposed a new set of optimizations based on an individual low-rank
approximation of the M2L operators; we refer to them as IA variants. As opposed to SA they directly lead
to the optimal low-rank representation for each operator. The overall number of flops is greater than for
SArcmp (which is strictly a lower bound on the number of flops). However, the advantage of the individual
treatment of the M2L operators is that we can exploit symmetries in their arrangement. This means that
all operators can be expressed as permutations of a subset. For example, in the case of the bbFMM (in
which the full interaction list has a constant size), we need to approximate and store only 16 instead of
316 operators. The remaining ones can be expressed as permutations thereof. This has a great impact on
the precomputation time and the memory requirement. Moreover it allows to express (again in the case
of the bbFMM) the at most 189 matrix-vector products (applications of the M2L operators) as at most 16
matrix-matrix products. We referred to this approach as the IAblk variant. It can then take advantage of
highly optimized implementations of matrix-matrix operations (e.g., the MKL [9]).
Let us conclude by comparing SArcmp and IAblk, the two variants that have the fastest running times.
IAblk wins if we compare precomputation time, required memory and runtime at levels having non-directional
expansions (bbFMM and low-frequency levels in dFMM). SArcmp wins only if we compare the runtime at
levels having directional expansions (high-frequency levels in dFMM). However, in order to identify the
optimal variant we have to distinguish two potential uses of the FMM as a numerical scheme to perform
fast matrix-vector products. 1) If we are interested in the result of a single matrix-vector product, a quick
precomputation is essential. However, 2) if we are looking for the iterative solution of an equation system
(many matrix-vector products), a fast running time of the M2L operator is crucial. Let us explain this
with an example. We take dFMM (with MKL) with accuracy Acc = 5 for the sphere. IAblk wins if we
are interested in the former use. The precomputation takes 0.4 s versus 69.1 s (for SArcmp) and the M2L
application takes 10.0 s versus 6.3 s, which sums up to 10.4 s versus 75.4 s. All other operators (P2P, P2M,
M2M, L2L and L2P) have nearly the same runtime in both cases, and their runtimes are negligible compared
to M2L. Looking at the latter use, SArcmp starts being faster if the iterative solution requires more than 19
22

matrix-vector products. For higher accuracies this threshold rises, e.g., for Acc = 6 it lies at 26 matrix-vector
products. In the case of bbFMM, IAblk is always optimal.

References
[1] M. Bebendorf. Hierarchical LU decomposition based preconditioners for BEM. Computing, 74:225–247,
2005.
[2] H. Cheng, W. Y. Crutchfield, Z. Gimbutas, L. Greengard, J. F. Ethridge, J. Huang, V. Rokhlin,
N. Yarvin, and J. Zhao. A wideband fast multipole method for the Helmholtz equation in three dimensions. Journal of Computational Physics, 216(1):300–325, 2006.
[3] E. Darve and P. Havé. Efficient fast multipole method for low-frequency scattering. Journal of Computational Physics, 197(1):341–363, 2004.
[4] 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., 16(1):1–17, March 1990.
[5] W. Fong and E. Darve. The black-box fast multipole method. Journal of Computational Physics, 228
(23):8712–8725, 2009.
[6] Z. Gimbutas and V. Rokhlin. A generalized fast multipole method for nonoscillatory kernels. SIAM
Journal on Scientific Computing, 24(3):796–817, March 2002.
[7] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. Journal of Computational
Physics, 73:325–348, 1987.
[8] L. Greengard and V. Rokhlin. A new version of the fast multipole method for the Laplace equation in
three dimensions. Acta Numerica, 6:229–269, 1997.
[9] Intel.
Intel
Math
Kernel
Library
(Intel
MKL)
10.3,
2012.
URL
http://software.intel.com/en-us/articles/intel-mkl/. [Online; accessed 28-August-2012].
[10] P. G. Martinsson and V. Rokhlin. An accelerated kernel-independent fast multipole method in one
dimension. SIAM Journal on Scientific Computing, 29(3):1160–1178, May 2007.
[11] J. C. Mason and D. C. Handscomb. Chebyshev Polynomials. Chapman & Hall/CRC, 2003.
[12] M. Messner, M. Schanz, and E. Darve. Fast directional multilevel summation for oscillatory kernels
based on chebyshev interpolation. Journal of Computational Physics, 231(4):1175 – 1196, 2012.
[13] V. Rokhlin. Diagonal forms of translation operators for the Helmholtz equation in three dimensions.
Applied and Computational Harmonic Analysis, 1:82–93, 1993.
[14] L. Ying, G. Biros, and D. Zorin. A kernel-independent adaptive fast multipole algorithm in two and
three dimensions. Journal of Computational Physics, 196(2):591 – 626, 2004.

23

