Fast Multipole Preconditioners for Sparse Matrices Arising from Elliptic Equations

Among optimal hierarchical algorithms for the computational solution of elliptic problems, the Fast Multipole Method (FMM) stands out for its adaptability to emerging architectures, having high arithmetic intensity, tunable accuracy, and relaxable global synchronization requirements. We demonstrate that, beyond its traditional use as a solver in problems for which explicit free-space kernel representations are available, the FMM has applicability as a preconditioner in finite domain elliptic boundary value problems, by equipping it with boundary integral capability for satisfying conditions at finite boundaries and by wrapping it in a Krylov method for extensibility to more general operators. Here, we do not discuss the well developed applications of FMM to implement matrix-vector multiplications in within Krylov solvers of boundary element methods. Instead, we propose using FMM for the volume-to-volume contribution of inhomogeneous Poisson-like problems, and the boundary integral is a small part of the overall computation. Our method may be used to precondition sparse matrices arising from finite difference/element discretizations, and can handle a broader range of scientific applications. Compared with multilevel methods, it is capable of comparable algebraic convergence rates down to the truncation error of the discretized PDE, and it has superior multicore and distributed memory scalability properties on commodity architecture supercomputers. Compared with other methods exploiting the low rank character of off-diagonal blocks of the dense resolvent operator, FMM-preconditioned Krylov iteration may reduce the amount of communication because it is matrix-free and exploits the tree structure of FMM. We describe our tests in reproducible detail with freely available codes and outline directions for further extensibility.


Introduction
Elliptic PDEs arise in a vast number of applications in scientific computing.A significant class of these involve the Laplace operator, which appears not only in potential calculations but also in, for example, Stokes and Navier-Stokes problems [21,Chapters 5 and 7], electron density computations [43, Part II] and reaction-convection-diffusion equations [37,Part IV].Consequently, the rapid solution of PDEs involving the Laplace operator is of wide interest.
Although many successful numerical methods for such PDEs exist, changing computer architectures necessitate new paradigms for computing and the development of new algorithms.Computer architectures of the future will favor algorithms with high concurrency, high arithmetic intensity (Flop/Byte), asynchronicity, and data locality.This trend is manifest on GPUs and coprocessors, where some algorithms are accelerated much less than others on the class of architectures that can be extended to extreme scale.There is always a balance between algorithmic efficiency in a convergence sense, and how well an algorithm scales on parallel architectures.This balance is shifting towards increased parallelism, even at the cost of increasing computation.There is therefore a strong incentive to choose/modify the algorithms that promise the best use of future architectures.Since the processor frequency has plateaued for the last decade, Moore's law holds continued promise only for those who are willing to make algorithmic changes.
Among the scientific applications ripe for reconsideration, those governed by elliptic PDEs will be among the most challenging.A common solution strategy for such systems is to discretize the partial differential equations by low-order finite element, finite volume or finite difference methods and then solve the resulting large, sparse linear system.However, elliptic systems are global in nature, and this conflicts with the requirements of future architectures, e.g.high concurrency, asynchronicity, and data locality.The linear solver must enable the transfer of information from one end of the domain to the other, either through successive local communications (as in many iterative methods), or a direct global communication (as in direct solvers with global recurrences and Krylov methods with global reductions).In either case, avoiding synchronization and reducing communication are the main challenges.There has been considerable effort in this direction in the dense linear algebra community [18].The directed-acyclic-graph-based technology developed in such efforts could be combined with iterative algorithms of optimal complexity for solving elliptic PDEs at extreme scale.
Scalable algorithms for solving elliptic PDEs tend to have a hierarchical structure, as in multigrid methods [56], fast multipole methods (FMM) [30], and H-matrices [33].This structure is crucial, not only for achieving optimal arith-metic complexity, but also for minimizing data movement.For example, a 3-D FFT requires O( √ P ) communication for the transpose between pencil-shaped subdomains on P processes [16], whereas these hierarchical methods require O(log P ) communication [41].This O(log P ) communication complexity is likely to be optimal for elliptic problems, since an appropriately coarsened form of a local representation must somehow arrive to all other parts of the domain for the elliptic equation to converge.In other words, an elliptic problem for which the solution is desired everywhere cannot have a communication complexity of O (1).However, a disadvantage of these hierarchical methods is that for certain problems we see slow convergence, or even divergence.
Krylov subspace methods provide another popular alternative to direct methods for general operators, although we note that methods such as Chebyshev semi-iteration can require even less communication when information about the spectrum of the coefficient matrix is known [23,Section 10.1.5],[24].Among the best known Krylov methods are the conjugate gradient method [35], MINRES [46] and GMRES [50], although a multitude of Krylov solvers are available in popular scalable solver libraries.The great advantage of these solvers is their robustness -for any consistent linear system there exists a Krylov method that will converge, in exact arithmetic, for sufficiently many iterations.However, the convergence rate of Krylov methods typically deteriorates as the discretization of an elliptic PDE is refined.Mesh-independent convergence for Krylov methods applied to systems from elliptic PDEs can often be recovered by preconditioning.Among the best performing preconditioners are the optimal hierarchical methods or, for multiphysics problems such as Stokes and Navier-Stokes equations, block preconditioners with these methods as components.By combining these hierarchical methods and Krylov subspace solvers we get the benefits of both approaches and obtain a linear solver that is fast but robust.These hierarchical methods all have multiple parameters for controlling the precision of the solution and are able to trade-off accuracy for speed, which is a useful feature for a preconditioner.Furthermore, in analogy to geometric multigrid and algebraic multigrid, H 2 -matrices can be thought of as an algebraic generalization of what FMMs do geometrically.There are advantages and disadvantages to using algebraic and geometric methods, and both have their place as preconditioners.
There has been some recent work on algebraic multigrids (AMG) in anticipation of the future hardware constraints mentioned above.Gahvari et al. developed a performance model for AMG and tested it on various HPC systems -Intrepid, Jaguar, Hera, Zeus, and Atlas [22].They found that network distance and contention were both substantial performance bottlenecks for AMG.Adams presents a low-memory matrix-free full multigrid (FMG) with a full approximation storage (FAS) [1].He revives an idea from the 1970s [9], which processes the multigrid algorithm vertically, and improves data locality and asynchronicity.Baker et al. compared the scalability of different smoothershybrid Gauss-Seidel, l 1 Gauss-Seidel, and Chebyshev polynomial, and showed that l 1 Gauss-Seidel and Chebychev smoothers scale much better [3].There is continuous effort in the multigrid community to adapt the algorithm according to future hardware constraints, and it is likely that multigrid will evolve to remain competitive.
On the other hand, performing a hierarchical low rank approximation (HLRA) of the off-diagonal blocks of a matrix leads to a whole new variety of O(N ) solvers or preconditioners.HLRA based methods include FMM [30], H-matrices [33], hierarchically semi-separable matrices (multifrontal methods) [13], and recursive skeletonization [36].These techniques can be applied to a dense matrix or the fill-in during a sparse direct solve, thus enabling an O(N ) matrix-vector multiplication of a N × N dense matrix or an O(N ) direct solve of a N × N sparse matrix to within a specified accuracy.These HLRA based methods require a decaying kernel which yields a block low-rank structure.The distinguishing features of the variants come in the way the low rank approximation is constructed -rank-revealing LU [47], rank-revealing QR [31], pivoted QR [39], truncated SVD [26], randomized SVD [42], adaptive cross approximation [48], hybrid cross approximation [8], and Chebychev interpolation [19] are all possibilities.Multipole/local expansions in the FMM constitute another way to construct the low rank approximations.In fact, many of the original developers of FMM are now working on these algebraic variants [29].
Literature on the HLRA based methods mentioned above mainly focuses on the error convergence of the low rank approximation and there is little investigation of the parallel scalability or direct comparison against multigrid.An exception is the work by Grasedyck et al. [27], where their H-LU preconditioner is compared with BoomerAMG, Pardiso, MUMPS, UMFPACK, SuperLU, and Spooles.However, their executions are serial, and show that their H-matrix code is not yet competitive with these other highly optimized libraries.
In the present work, we consider the Laplace and Stokes equations and devise highly scalable preconditioners for these problems.Our Poisson preconditioner is based on a boundary element method in which matrix-vector multiplies are performed using FMM; the result is an O(N ) preconditioner that is scalable.For the Stokes problem, we apply a block diagonal preconditioner, in which our Poisson preconditioner is combined with a simple diagonal matrix.Such FMM based preconditioners were first proposed by Sambavaram et al. [51].Such methods lacked practical motivation when flops were expensive, since they turn a sparse matrix into a dense matrix of the same size before hierarchically grouping the off-diagonal blocks.But in a world of cheap flops, the notion of a "compute-bound preconditioner" sounds more attractive.In the present work, we perform scalability benchmarks and compare the time-to-solution with state-of-the-art multigrid methods such as BoomerAMG in a high performance computing environment.
The rest of the manuscript is organized as follows.In Section 2 we present the model problems and in Section 3 we give an overview of Krylov subspace methods and preconditioning.The basis of our preconditioner is a boundary element method that is discussed in Section 4 and the FMM, the essential kernel that makes our method efficient and scalable, is described in Section 5. Our numerical results in Section 6 examine the convergence rates of FMM and multigrid for small Poisson and Stokes problems.Then, in Section 7 we scale up the Poisson problem tests and perform strong scalability runs, where we compare the time-to-solution against BoomerAMG [34] on up to 1024 cores.Our conclusions are given in Section 8.

Model problems
In this section we introduce the Poisson and Stokes model problems we wish to solve and describe properties of the linear systems that result from their discretization.We focus on low-order finite elements but note that discretization by low-order finite difference or finite volume methods give linear systems with similar properties.

Poisson model problem
The model Poisson problem we wish to solve is of the form where Ω ∈ R d , d = 2, 3 is a bounded connected domain with piecewise smooth boundary Γ, f is a forcing term and g defines the Dirichlet boundary condition.
Discretization of (1) by finite elements or finite differences leads to a large, sparse linear system of the form where A ∈ R N ×N is the stiffness matrix and b ∈ R N contains the forcing and boundary data.The matrix A is symmetric positive definite and its eigenvalues depend on the mesh size, which we denote by h, as is typical of discretizations of elliptic PDEs.In particular, the condition number κ = λ max (A)/λ min (A),  [7], [21]: Discretizing (3) by a stabilized1 finite element or finite difference approximation leads to the symmetric saddle point system where A ∈ R N ×N is the vector-Laplacian, a block diagonal matrix with blocks equal to the stiffness matrix from (2), B ∈ R M ×N is the discrete divergence matrix, C ∈ R M ×M is the symmetric positive definite pressure mass matrix and f ∈ R N and g ∈ R M contain the Dirichlet boundary data.
The matrix A is symmetric indefinite and the presence of the stiffness matrix means that the condition number of A increases as the mesh is refined.However, as we will see in the next section, the key ingredient in a preconditioner for A that mitigates this mesh dependence is a good preconditioner for the Poisson problem.This allows us to use our preconditioner for the Poisson problem in this more complicated fluid dynamics problem as well.
3 Iterative solvers and preconditioning

Krylov Subspace Methods
Large, sparse systems of the form (2) are often solved by Krylov subspace methods.We focus here on two Krylov methods: the conjugate gradient method (CG) [35] for systems with symmetric positive definite coefficient matrices and MINRES [46] for systems with symmetric indefinite matrices.For implementation and convergence details, we refer the reader to the books by Greenbaum [28] and Saad [49].
The convergence of these Krylov subspace methods depends on the spectrum of the coefficient matrix which for the Poisson and Stokes problems, as well as other elliptic PDEs, deteriorates as the mesh is refined.This dependence is removed by preconditioning.In the case of the Poisson problem (2), we can conceptually think of solving the equivalent linear system P −1 Ax = P −1 b (left preconditioning), or AP −1 y = b, with P −1 y = x (right preconditioning), for some P −1 ∈ R N ×N , and analogously for the Stokes equations ( 4).However, when the coefficient matrix is symmetric, we would like to preserve this property when preconditioning; this can be achieved by using a symmetric positive definite preconditioner [21,Chapters 2 and 6].We also note that in practice we never need P −1 explicitly but only the action of this matrix on a vector.This enables us to use matrix-free approaches such as multigrid or the fast multipole method.
Many preconditioners for the Poisson problem reduce the number of iterations, with geometric and algebraic multigrid among the most effective strategies [21], [56].However, to achieve a lower time-to-solution than can by obtained for the original system, it is also necessary to choose a preconditioner that can be cheaply applied at each iteration.Both geometric and algebraic multigrid methods are O(N ), and therefore exhibit good performance on machines and problems for which computation is expensive.However, stresses arise in parallel applications as discussed in the introduction.
For Stokes problems we consider the block diagonal preconditioner where P A ∈ R N ×N and P S ∈ R M ×M are symmetric positive definite matrices.The advantage of this preconditioner is that there is no coupling between the blocks, so P is scalable provided the blocks P A and P S are.
Appropriate choices for P A and P S have been well studied and it is known that mesh-independent convergence of MINRES can be recovered when P A is spectrally equivalent to A in (4) and P S is spectrally equivalent to the pressure mass matrix Q ∈ R M ×M [11], [21,Chapter 6].These spectral equivalence requirements imply that the eigenvalues of P −1 A A and P −1 S Q are bounded in an interval on the positive real line independently of the mesh width h.
It typically suffices to use the diagonal of Q [21, Chapter 6], [61] or a few steps of Chebyshev semi-iteration [60] for P S .Moreover, the diagonal matrix is extremely parallelizable.Thus, the key to obtaining a good preconditioner for A is to approximate the vector Laplacian effectively.This is typically the most computationally intensive part of the preconditioning process, since in most cases M N .

The FMM-BEM preconditioner
In this paper we propose an alternative preconditioner for Poisson and Stokes problems that heavily utilizes the fast multipole method (FMM).The FMM is O(N ) with compute intensive inner kernels.It has a hierarchical data structure that allows asynchronous communication and execution.These features make the FMM a promising preconditioner for large scale problems on future computer architectures.We show that this preconditioner is effective at improving the convergence of these Krylov subspace methods, and is effectively parallelized on today's highly distributed architectures.
The FMM in its original form relies on free-space Green's functions and is able to solve problems with free-field boundary conditions.In Section 4 the FMM preconditioner is extended to Dirichlet, Neumann or Robin boundary conditions for arbitrary geometries by coupling it with a boundary element method (BEM).Our approach uses the FMM as a preconditioner inside a sparse matrix solver and the BEM solve is inside the preconditioner.Numerous previous studies use FMM for the matrix-vector multiplication inside the Krylov solver for the dense matrix arising from the boundary element discretization.In the present method we are calculating problems with non-zero sources in the volume, and the FMM is used to calculate the volume-to-volume contribution.This means we are performing the action of an N ×N dense matrix-vector multiplication, where N is the number of points in the volume (not the boundary).Additionally, as discussed in Section 4.3, it is possible to extend the boundary element method to problems with variable diffusion coefficients, particularly since low accuracy solves are often sufficient in preconditioning.
Figure 1 shows the flow of calculation of our FMM-BEM preconditioner within the conjugate gradient method; its role in other Krylov solvers is similar.The FMM is used to approximate the matrix-vector multiplication of A −1 within the preconditioner.The BEM solver adapts the FMM to finitely applied boundary conditions.During each step of the iteration, the u vector from the previous iteration is used to determine ∂u/∂n at the boundary from (8), then ( 9) is used to compute the new u in the domain Ω.
4 Boundary Element Method

Formulation
We use a standard Galerkin boundary element method [52] with volume contributions to solve the Poisson equation.A brief description of the formulation is given here.Applying Green's third identity to (1a) gives where G is the Green's function of the Laplace operator, ∂ ∂n is the derivative in the outward normal direction, and Γ is the boundary.Following the definition of the Green's function ∇ 2 G = −δ, the third term in (6) becomes Therefore, we may solve the constant coefficient inhomogeneous Poisson problem by solving the following set of equations As an example, consider the case where Dirichlet boundary conditions are prescribed on ∂Ω.The unknowns are ∂u/∂n on Γ and u in Ω\Γ, where (8) solves for the former and ( 9) can be used to determine the latter.For Neumann boundary conditions one can simply switch the two boundary integral terms in (8) and solve for u instead of ∂u/∂n.In either case, we obtain both u and ∂u/∂n at each point on the boundary, then calculate (9) to obtain u at the internal points.The last term in (9) takes up most of the calculation time since it is a volume integral for every point in the volume, whereas other terms are either for every point on the boundary or are boundary integrals.

Discretization
The integrals in equations ( 8) and ( 9) are discretized in a similar fashion to finite element methods.In the following description of the discretization process, we will use the term on the left hand side in (8) as an example.The first step is to break the global integral into a discrete sum of piecewise local integrals over each element where N Γ is the number of boundary nodes.These piecewise integrals are performed by using quadratures over the basis functions [52].In the present case, we use constant elements so there are no nodal points at the corners of the square domain for the tests in Sections 6 and 7.By applying this discretization technique to all terms in (8) we obtain , where N Ω is the number of internal nodes.All values on the right hand side are known, and ∂u/∂n at the boundary is determined by solving the linear system.Similarly, we apply the discretization to (9) to have .
At this point, all values on the right hand side are known so one can perform three matrix-vector multiplications to obtain u at the internal nodes, and the solution to the original Poisson equation (1a).The third term on the right hand side involves an N Ω × N Ω matrix, and is the dominant part of the computational load.This matrix-vector multiplication can be approximated in O(N ) time by using the FMM described in Section 5. We also use the FMM to accelerate all other matrix-vector multiplications.

Variable coefficient problems
A natural question that arises is how to extend the boundary element method, which is the basis of our preconditioner, to problems with variable diffusion coefficients, i.e., problems of the type where x ∈ R d , d = 2, 3 and a ∈ R.
Several strategies for extending boundary element methods to problems with variable diffusion coefficients have been proposed (see, for example, the thesis of Brunton [10,Chapter 3]).Additionally, in this preconditioner setting we may not need to capture the variation in the diffusion coefficient to a high degree of accuracy; for a similar discussion in the context of additive Schwarz preconditioners see, for example, Graham et al. [25].
Although analytic fundamental solutions can sometimes be found for problems with variable diffusion (see, e.g., Cheng [14] and Clements [15]), in most cases numerical techniques are employed.One popular method is to introduce a number of subdomains, on each of which the diffusion coefficient is approximated by a constant function [40,57].
A second option is to split the differential operator into a part for which a fundamental solution exists and another which becomes part of the source term.Specifically, starting from (11), a similar approach to that described in Banerjee [6] and Cheng [14] leads to where again G is the standard fundamental solution for the Laplace operator, i.e., not the fundamental solution for (11).We can then proceed as described above for (6).
5 Fast Multipole Method

Introduction to FMM
The last term in Eq. ( 9) when discretized, has the form where i = 1, 2, ..., N Ω .If we calculate this equation directly, it will require O(N 2 ) operations.In Figure 2, we show a schematic of how the fast multipole method is able to calculate this in O(N ) operations.Figures 2(a) and 2(b) show how the source particles (red) interact with the target particles (blue) for the direct method and FMM, respectively.In the direct method, all source particles interact with all target particles directly.In the FMM, the source particles are first converted to multipole expansions using the P2M (particle to multipole) kernel.Figure 2(c) shows the corresponding geometric view of the hierarchical domain decomposition of the particle distribution.Then, multipole expansions are aggregated into larger groups using the M2M (multipole to multipole) kernel.Following that, the multipole expansions are translated to local expansions between well-separated cells using the M2L (multipole to local) kernel.Both Figures 2(b) and 2(c) show that the larger cells interact if they are significantly far away, and smaller cells may interact with slightly closer cells.The direct neighbors between the smallest cells are calculated using the P2P (particle to particle) kernel, which is equivalent to the direct method between a selected group of particles.Then, the local expansions of the larger cells are translated to smaller cells using the L2L (local to local) kernel.Finally, the local expansions at the smallest cells are translated into the potential on each particle using the L2P (local to particle) kernel.The mathematical formulae for these kernels will be given in Section 5.3.
In order to perform the FMM calculation mentioned above, one must first de-compose the domain in a hierarchical manner.It is common to use an octree in 3-D and quad-tree in 2-D, where the domain is split by its geometrical centerline.The splitting is performed recursively until the number of particles per cell reaches a prescribed threshold.The splitting is usually performed adaptively, so that the densely populated areas result in a deeper branching of the tree.A common requirement in FMMs is that these cells must be isotropic (cubes or squares and not rectangles), since they are used as units for measuring the well-separatedness as shown in Figure 2(c) during the M2L interaction.However, our FMM does not use the size of cells to measure the distance between them and allows the cells to be of any shape as long as they can be hierarchically grouped into a tree structure.Once the tree structure is constructed, it is trivial to find parent-child relationships between the cells/particles.This relation is all that is necessary for performing P2M, M2M, L2L, and L2P kernels.However, for the M2L and P2P kernels one must identify a group of well-separated and neighboring cells, respectively.We will describe an efficient method for finding well-separated cells in the following subsection.

Dual Tree Traversal
The simplest method for finding well-separated pairs of cells in the FMM is to "loop over all target cells and find their parent's neighbor's children that are non-neighbors", as shown by Greengard and Rokhlin [30].A scheme that permits the interaction of cells at different levels for an adaptive tree was introduced by Carrier et al. [12].This scheme is used in many modern FMM codes, and is sometimes called the UVWX-list [41].Another scheme to find well-separated pair of cells is to "simultaneously traverse the target and source tree while applying a multipole acceptance criterion", as shown by Warren and Salmon [59].Teng [55] showed that this dual tree traversal can produce interaction pairs that are almost identical to the adaptive interaction list by Carrier et al. [12].A concise explanation and optimized implementation of the dual tree traversal is provided by Dehnen [17].
The dual tree traversal has many favorable properties compared to the explicit construction of interaction lists.First of all, the definition of well-separatedness can be defined quite flexibly.For example, if one were to construct explicit interaction lists by extending the definition of neighbors from 3 × 3 × 3 to 5 × 5 × 5 using the traditional scheme, the M2L list size will increase rapidly from 6 3 − 3 3 = 189 to 10 3 − 5 3 = 875 in 3-D, which is never faster for any number of expansions.On the other hand, the dual tree traversal can adjust the definition of neighbors much more flexibly and the equivalent interaction list always has a spherical shape.(We say "equivalent interaction list" because there is no explicit interaction list construction in the dual tree traversal.) The cells no longer need to be cubic, since the cells themselves are not used to measure the proximity of cells.The cells can be any shape or size -even something like a hierarchical K-means.Of course, the explicit interaction list construction can be modified to include more flexibility, too [32].However, the resulting code becomes much more complicated than the dual tree traversal, which is literally a few lines of code. 2 This simplicity is a large advantage on its own.Furthermore, the parallel version of the dual tree traversal simply traverses the local tree for the target with the local essential tree [58] for the sources, so the serial dual tree traversal code can be used once the local essential tree is assembled.
A possible (but unlikely) limitation of dual tree traversals is the loss of explicit parallelism -it has no loops.It would not be possible to simply use an OpenMP "parallel for" directive to parallelize the dual tree traversal.In contrast, the traditional schemes always have an outer loop over the target cells, which can be easily parallelized and dynamically load balanced with OpenMP directives.However, this is not an issue since task based parallelization tools such as Intel Thread Building Blocks (TBB) can be used to parallelize the dual tree traversal.With the help of these tools, tasks are spawned as the tree is traversed and dispatched to idle threads dynamically.In doing so, we not only assure load-balance but also data-locality, so it may actually end up being a superior solution than parallelizing "for loops" with OpenMP, especially on NUMA architectures.
Considering the advantages mentioned above, we have decided to use the dual tree traversal in our current work.This allows us to perform low accuracy optimizations by adjusting the multipole acceptance criterion without increasing the order of expansions too much, which is the secret to our speed [62].These low accuracy optimizations can give the FMM a performance boost when used as a preconditioner.

Multipole Expansions
For the 2-D Laplace equation, the free space Green's function has the form where r ij = |x i − x j | is the distance between point i and point j.Eq. ( 12) can be written as where (z) represents the real part of z. Figure 3 shows the decomposition of vector x ij into five parts,  [12].We denote the nth order multipole expansion coefficient at x as M n (x), and the nth order local expansion coefficient as L n (x), where n = 0, 1, ..., p − 1 for a pth order truncation of the series.
(1) P2M from particle at x j to multipole expansion at x µ , (2) M2M from multipole expansion at x µ to multipole expansion at x M , (3) M2L from multipole expansion at x M to local expansion at x Λ , (4) L2L from local expansion at x Λ to local expansion at x λ , (5) L2P from local expansion at x λ to particle at x i , For the P2M, M2M, and M2L kernels, the first term requires special treatment.
The expansions are truncated at order p, so the accuracy of the FMM can be controlled by adjusting p.When recurrence relations are used to calculate the powers of z and the combinations they can be calculated at the cost of one multiplication per inner loop (k loop) iteration.In our implementation, we do not construct any matrices during the calculation of these kernels.The P2P kernel is vectorized with the use of SIMD intrinsics, and the log() function is calculated using a polynomial fit for log 2 (x)/(x − 1) using SIMD.

Numerical results
In this section we demonstrate the potential of the FMM-based preconditioner by applying it to a number of test problems and comparing it with standard preconditioners.The primary aim is to assess the effectiveness of the preconditioner at reducing the number of Krylov subspace iterations that are required for convergence to a given tolerance.Additionally, we seek to ascertain whether mesh independence is achieved.We defer reporting on performance to Section 7. Accordingly, we choose problems that are small enough to enable solution by Matlab.
Our Poisson problems are all two dimensional and include examples with homogeneous and inhomogeneous Dirichlet boundary conditions.We additionally present a two-dimensional Stokes flow problem and show that, as predicted, combining the FMM-based Poisson preconditioner with a block diagonal matrix gives an effective preconditioner for the saddle point problem (4).
Throughout, our stopping criterion is a decrease in the relative residual of six orders of magnitude.If such a decrease is not achieved after maxit iterations the computations are terminated; this is denoted by '-' in the tables.This maximum number of iterations is stated for each problem below.For all problems and preconditioners the initial iterate is zero at entries corresponding to points in the domain interior and equal to b at points corresponding to the  1 Preconditioned CG iterations for the relative residual to reduce by six orders of magnitude for the problem with −∇ 2 u = 1 and homogeneous boundary conditions.
boundary; this ensures that the residual vector is always zero at the boundary.Careful implementation of the preconditioned Krylov subspace method would improve efficiency since the zero elements in the residual make certain computations redundant; we do not consider this here.

The Poisson equation
We first test our preconditioner on three two-dimensional Poisson problems with a constant diffusion coefficient on the domain on [−1, 1] 2 .We discretize the problems by Q 1 finite elements using IFISS [20], [54], with default settings.Our fast multipole preconditioner is compared with the incomplete Cholesky (IC) factorization [44] with zero fill implemented in Matlab and the algebraic multigrid (AMG) and geometric multigrid (GMG) methods in IFISS.Within the GMG preconditioner we select point-damped Jacobi as a smoother instead of the default ILU, which is less amenable to parallelization.Otherwise, default settings for both multigrid methods are used.For all preconditioners, maxit = 20 and we apply preconditioned conjugate gradients.
Our first example is the first reference problem in Elman et al. [21, Section 1.1] for which Table 1 lists the preconditioned CG iterations for each preconditioner applied.The FMM preconditioner, as well as GMG and AMG appear to give meshindependent convergence, although the incomplete Cholesky factorization does not.The slight increase in the number of FMM iterations on the finest grid is peculiar to this problem, as is clear from subsequent examples, and with the exception of this grid the AMG and FMM preconditioners are equally competitive.
In Table 2 we plot the eigenvalues of the FMM preconditioned stiffness matrix for h = 2 −4 , 2 −5 and 2 −6 .It is clear that the smallest eigenvalue of A decreases h λ min (A) λ max (A) κ(A) λ min (P −1 A) λ max (P −1 A) κ(P −1 A)  2 Smallest (λ min ) and largest (λ max ) eigenvalues and condition number (κ) of the stiffness matrix A and FMM-preconditioned matrix P −1 A for the problem with −∇ 2 u = 1 and homogeneous boundary conditions.
as the mesh is refined; this is particularly problematic for Krylov subspace methods, since small eigenvalues can significantly hamper convergence.However, the eigenvalues of the FMM-preconditioned matrix are bounded away from the origin in a small interval that does not increase in size as the mesh is refined.This hints at spectral equivalence between the FMM-based preconditioner and the stiffness matrix.The condition number appears to be bounded, which is unsurprising given the mesh-independent convergence observed.
Our second example is the third reference problem from Elman et al.
From Table 3 we find that, similarly to the previous problem, the FMM preconditioner and both multigrid preconditioners are mesh independent but the Cholesky preconditioner is not.The FMM preconditioner is also competitive with the multigrid methods, requiring only one more CG iteration.Thus, on systems on which applying the FMM preconditioner is significantly faster than applying the multigrid preconditioners, we will achieve a faster time-tosolution with the former.We note that the eigenvalues and condition numbers obtained for the FMM preconditioned stiffness matrix are the same as those computed for the previous example.
The final problem we consider in this section is the Poisson problem with solution u(x, y) = x 2 + y 2 on [−1, 1] 2 , which has forcing term f ≡ −4 in the domain and inhomogeneous Dirichlet boundary conditions.The convergence results for this problem, given in Table 4, are similar to those for the previous problems.They show that the FMM preconditioner gives mesh independent convergence and is competitive with AMG and GMG.We also obtain the same eigenvalue results as for the previous examples.3 Preconditioned CG iterations for the relative residual to reduce by six orders of magnitude for the problem with −∇ 2 u = 0 and inhomogeneous boundary conditions.4 Preconditioned CG iterations for the relative residual to reduce by six orders of magnitude for the problem with −∇ 2 u = −4 and inhomogeneous boundary conditions.

Effect of FMM precision on convergence
For the results shown above, the FMM precision was set to preserve six significant digits.However, the FMM can be accelerated further by trading precision for speed.Since we are using the FMM as a preconditioner, the accuracy requirements are somewhat lower than that of general applications of FMM.Although this balance between the accuracy and speed of FMM is a critical factor for evaluating the usefulness of FMM as a preconditioner, the relation between the FMM precision and convergence rate has not been studied previously.
In Figure 4 the relative residual at each CG iteration is plotted against the number of iterations for FMM, AMG, GMG, and IC.The problem is the same as in Table 1.Three cases of FMM are used with six, four, and two significant digits of accuracy, respectively.The = 10 −6 case corresponds to the condition for the tests in Tables 1-4.Decreasing the FMM accuracy to four digits has little effect during the first few iterations, but slows down the convergence near the end.Decreasing the FMM accuracy further to two digits slows down the convergence significantly, but is still much better than the Increasing the precision of the FMM past six digits did not result in any noticeable improvement because truncation error begins to dominate.We are preconditioning a matrix resulting from a FEM discretization by using a integral equation with Green's function kernels.Each has its own error, below which algebraic error need not be reduced.We show in Figure 5 the convergence of spatial discretization error for the FEM and BEM approaches.We use the same reference problem as in Table 1, which has an analytical solution.The discretization error is measured by taking the relative L 2 norm of the difference between the analytical solution and the individual numerical solutions.We see that the FEM is second order and BEM is first order.The five different values of ∆x correspond to h = {2 −4 , 2 −5 , 2 −6 , 2 −7 , 2 −8 }, which were used in the previous experiments.For the current range of grid spacing, the discrepancies between the FEM and BEM truncation error is in the range of 10 −3 to 10 −4 .ments in Matlab using IFISS with default settings.As described in Section 3, by combining a stiffness matrix preconditioner P A with the diagonal of the pressure mass matrix P S , an effective preconditioner (5) for the saddle point system (4) is obtained.Here, we are interested in using the FMM preconditioner for P A , and we compare its performance with AMG and GMG.We do not consider the incomplete Cholesky factorization of A because of its poor performance on the stiffness matrix (see Tables 1, 3 and 4).We set maxit = 50 and apply preconditioned MINRES to the saddle point system.

Stokes problem
As for the Poisson problem, the FMM-based preconditioner provides a meshindependent preconditioner that is comparable to algebraic and geometric multigrid.Although three or four more iterations are required by the FMM preconditioner than the AMG preconditioner, if each iteration is faster the time-to-solution may well be lower.

Performance analysis
In this section we evaluate the performance of the FMM-based preconditioner by comparing its time-to-solution to an algebraic multigrid code Boomer-AMG.We have implemented our FMM-preconditioner into PETSc [4] , [5] via PetIGA [45].PetIGA is a software layer that sits on top of PETSc that facili- The time-to-solution is plotted against the problem size N in Figure 6.Since we are using PETSc, it is trivial to change the preconditioner to AMG by passing the option "--pc type hypre" during runtime.Therefore, the timeto-solution of BoomerAMG is shown as a reference in the same figure.For BoomerAMG we compared different relaxation, coarsening, and interpolation methods and found that "-pc hypre boomeramg relax type all backward-SOR/Jacobi -pc hypre boomeramg coarsen type modifiedRuge-Stueben -pc hypre bommeramg interp type classical" gives the best performance.
Both FMM and AMG runs are serial, where we used a single MPI process and a single thread.A majority of the time goes into the setup of the preconditioner "PCSetUp" and the actual preconditioning "PCApply", so only these events are shown in the legend.The "PCSetUp" is called only once for the entire run, while "PCApply" is called every iteration.For the present runs, both FMM and AMG required six iterations for the relative residual to drop six digits, so all runs are calling "PCApply" six times.The order of expansion for the FMM is set to p = 6 and θ = 0.4, which gives about six significant digits of accuracy.With this accuracy for the FMM, we are still able to converge in six iterations.The P2P kernel in the FMM code is performed in single precision using SIMD intrinsics, but this does not prevent us from reaching the required accuracy of six significant digits because we use Kahan's summation technique [38] for the reduction.
By taking a closer look at Figure 6, one can see that both the FMM and AMG show O(N ) asymptotic behavior.The FMM seems to have a slower preconditioning time, but a much faster setup time compared to AMG.The FMM also has a constant overhead which becomes evident when N is small.In summary, the time-to-solution of the FMM is approximately an order of magnitude larger than that of AMG for the serial runs.This is consistent with our intuition that FMM is not the preconditioner of choice for solving small problems on a single core.We will show in the following section that the FMM becomes competitive when scalability comes into the picture.

Parallel results
Using the same Poisson problem, we now compare the performance of FMM and AMG for parallel runs on Stampede.We also compare with a sparse direct solver MUMPS by invoking at runtime "-ksp type preonly -pc type lu -pc factor mat solver package mumps".
The strong scaling of FMM, AMG, and MUMPS are shown in Figure 7.We use the largest grid size in the previous runs N = 4096 2 .Stampede has 16 cores per node so all runs first parallelize among these cores and then among the nodes after the 16 cores are filled.The FMM strong scales quite well up to 1024 cores, while the parallel efficiency of AMG starts to decrease after 128 cores.The sparse direct solver has a much larger time-to-solution even on a single core, and is much less scalable than the other two hierarchical preconditioners.For this particular Poisson problem on this particular machine using this particular FMM code we see an advantage over BoomerAMG past 512 cores.

Comparison with other approaches
We mentioned in the introduction that the FMM is one of many hierarchical low rank approximation methods.In the present study we chose to use FMM instead of the other algebraic variants simply because one of the authors had already developed a highly optimized FMM code.Here, we justify this choice by comparing with the most recent results for H-matrices, hierarchically semiseparable matrices, and recursive skeletonization methods.We cannot perform a direct comparison in the same environment since some of these codes are not publicly available.Furthermore, most of these papers fail to mention the type of processor they used, so it is difficult to match their environment for our FMM runs.However, a rough comparison against the timings reported in their publication may be useful to some readers.We selected only the latest publications to compare with so that the processor speed would not be too different from ours.
As a point of reference, our FMM can solve the Poisson problem with N = 1024 2 mesh points in a little less than 3 seconds on a single core of a Xeon E5-2680, 2.7 GHz CPU.Note that this is the total time to solve the linear system, including the initial setup time and all iterations.Ambikasaran and Darve [2] solve a linear system based on Gaussian radial basis function interpolation, where the same N = 1024 2 system takes approximately 3000 seconds in total to assemble and solve.If we correct for the difference in precision by noticing that rank 30 takes about three times longer than rank 15 (which is closer to our accuracy), the equivalent time would be 1000 seconds.
Ho and Greengard [36] show an example of their recursive skeletonization technique by using it as a generalized FMM for the Laplace kernel.The N = 131, 072 2-D volume case with precision = 10 −9 takes 340 seconds for the matrix-compression, and the matrix-vector multiplication time is negligible compared to that.Correcting for the problem size will increase the calculation time by 1, 048, 576/131, 072 = 8, and correcting for the precision should roughly decrease the calculation time by three, similar to the previous case.Therefore, the estimated time-to-solution for an equivalent problem size and precision should be around 900 seconds.
Schmitz and Ying [53] demonstrate the performance of their multifrontal method with hierarchical matrices by solving a Poisson equation on a uniform 2-D mesh with zero Dirichlet boundary conditions.They show results for N = 1024 2 and r = 10 −6 in their table, which is directly comparable to our results.For this case, their setup time was 71 seconds and the solve time was 2 seconds.This is more than an order of magnitude faster than the two previous cases, but more than an order of magnitude slower than our current results.

Extension to 3-D
The results presented in this article are only for two dimensional problems.A natural question that arises is whether the extension to 3-D is straightforward, and whether FMM will still be competitive as a preconditioner or not.Our results showed that a dominant part of the calculation time for the FMM preconditioner is the "PCApply" stage, which is the dual tree traversal for calculation of M2L and P2P kernels.For 3-D kernels, the M2L operation is much more complicated so the calculation time of the FMM will increase, even for the same number of unknowns N .
Figure 8 shows the calculation time of our 2-D FMM and 3-D FMM, both for the Laplace kernel with four significant digits of accuracy on a single core of a Xeon E5-2680, 2.7 GHz CPU.The problem size N varies from 10 5 to 10 7 .We see that the 3-D FMM is about an order of magnitude slower than the 2-D FMM for the same problem size.The curse of dimensionality exists for FMMs as well, but it is not clear how competitive our 3-D FMM preconditioner will be until we perform a direct comparison against other preconditioners in 3-D.
The logical next step is to actually perform these 3-D comparisons.

Conclusions
The Fast Multipole Method, originally developed as a free-standing solver, can be effectively combined with Krylov iteration as a scalable and highly performant preconditioner for traditional low-order finite discretizations of elliptic boundary value problems.In model problems it performs similarly to algebraic multigrid in convergence rate, while excelling in scalings where AMG becomes memory bandwidth-bound locally and/or synchronization-bound globally.Additional algorithmic development and additional testing of implementations on emerging architectures are necessary to more fully define the niche in which FMM is the preconditioner of choice.Two of the extensions required relative to our current results are to coefficient variability and in the order of discretization of the boundary integral technique.No preconditioner considered in isolation can address the fundamental architectural challenges of Krylov methods for sparse linear systems, which are being simultaneously adapted to a less synchronization tolerant computational environments, but it is important to address the bottlenecks of preconditioning this most popular class of solvers by making a wide variety of tunable preconditioners available and better integrating them into the overall solver.Fast multipole-based preconditioners are demonstrably ready to play an important role in the migration of sparse iterative solvers to the exascale.

Fig. 1 .
Fig. 1.Flow chart of the FMM-BEM preconditioner within the conjugate gradient method.

Fig. 2 .
Fig. 2. Schematic of Fast Multipole Method.(a) shows the interactions for a O(N 2 ) direct method.(b) shows the interactions for the O(N ) FMM, describing the type of interaction between elements in the tree data structure.(c) shows the same FMM kernels as in (b), but from a geometric point of view of the hierarchical domain decomposition.

Fig. 3 .
Fig.3.Decomposition of the distance vector x ij = x i − x j into five parts, that correspond to the five stages P2M, M2M, M2L, L2L, and L2P in the FMM.

Fig. 4 .
Fig. 4. Convergence rate of the FMM preconditioner with different precision, plotted along with algebraic multigrid, geometric multigrid, and incomplete Cholesky preconditioners.The represents the precision of the FMM, where = 10 −6 corresponds to six significant digits of accuracy.incompleteCholesky.

FinallyFig. 5 .
Fig. 5. Convergence of spatial discretization error for the FEM and BEM.The relative L 2 norm of the difference between the analytical solution is plotted against the grid spacing ∆x.

Fig. 8 .
Fig. 8. Calculation time of 2-D and 3-D FMM for the same problem size.
where λ and Λ are the center of local expansions and µ and M are the center of multipole expansions.The lower case is used for the smaller cells and upper case is used for the larger cells.When assuming the relation |x ΛM | > |x iλ + x λΛ | + |x M µ + x µj | the following FMM approximations are valid