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 within Krylov solvers of boundary element methods. Instead, we propose using FMM for the volume-to-volume contribution of inhomogeneous Poisson-like problems, where 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 multigrid methods, it is capable of comparable algebraic convergence rates down to the truncation error of the discretized PDE, and it offers potentially 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.Comment: 17 pages, 9 figure

A multigrid scheme for 3D Monge-Amp\`ere equations

The elliptic Monge-Amp\`ere equation is a fully nonlinear partial differential equation which has been the focus of increasing attention from the scientific computing community. Fast three dimensional solvers are needed, for example in medical image registration but are not yet available. We build fast solvers for smooth solutions in three dimensions using a nonlinear full-approximation storage multigrid method. Starting from a second-order accurate centered finite difference approximation, we present a nonlinear Gauss-Seidel iterative method which has a mechanism for selecting the convex solution of the equation. The iterative method is used as an effective smoother, combined with the full-approximation storage multigrid method. Numerical experiments are provided to validate the accuracy of the finite difference scheme and illustrate the computational efficiency of the proposed multigrid solver.Comment: 18 pages, 1 figure, 7 tables, 41 references. Accepted by International Journal of Computer Mathematics (published online: 21 Nov 2016

A fast and memory-efficient spectral Galerkin scheme for distributed elliptic optimal control problems

Many scientific and engineering challenges can be formulated as optimization problems which are constrained by partial differential equations (PDEs). These include inverse problems, control problems, and design problems. As a major challenge, the associated optimization procedures are inherently large-scale. To ensure computational tractability, the design of efficient and robust iterative methods becomes imperative. To meet this challenge, this paper introduces a fast and memory-efficient preconditioned iterative scheme for a class of distributed optimal control problems governed by convection-diffusion-reaction (CDR) equations. As an alternative to low-order discretizations and Schur-complement block preconditioners, the scheme combines a high-order spectral Galerkin method with an efficient preconditioner tailored specifically for the CDR application. The preconditioner is matrix-free and can be applied within linear complexity where the proportionality constant is small. Numerical results demonstrate that the preconditioner is ideal in the sense that appropriate Krylov subspace methods converge within a low number of iterations, independently of the problem size and the Tikhonov regularization parameter

High-performance modeling acoustic and elastic waves using the Parallel Dichotomy Algorithm

A high-performance parallel algorithm is proposed for modeling the propagation of acoustic and elastic waves in inhomogeneous media. An initial boundary-value problem is replaced by a series of boundary-value problems for a constant elliptic operator and different right-hand sides via the integral Laguerre transform. It is proposed to solve difference equations by the conjugate gradient method for acoustic equations and by the GMRES$(k)$ method for modeling elastic waves. A preconditioning operator was the Laplace operator that is inverted using the variable separation method. The novelty of the proposed algorithm is using the Dichotomy Algorithm (Terekhov, 2010), which was designed for solving a series of tridiagonal systems of linear equations, in the context of the preconditioning operator inversion. Via considering analytical solutions, it is shown that modeling wave processes for long instants of time requires high-resolution meshes. The proposed parallel fine-mesh algorithm enabled to solve real application seismic problems in acceptable time and with high accuracy. By solving model problems, it is demonstrated that the considered parallel algorithm possesses high performance and efficiency over a wide range of the number of processors (from 2 to 8192).Comment: The formula (2) has been correcte

An Efficient Coarse Grid Projection Method for Quasigeostrophic Models of Large-Scale Ocean Circulation

This paper puts forth a coarse grid projection (CGP) multiscale method to accelerate computations of quasigeostrophic (QG) models for large scale ocean circulation. These models require solving an elliptic sub-problem at each time step, which takes the bulk of the computational time. The method we propose here is a modular approach that facilitates data transfer with simple interpolations and uses black-box solvers for solving the elliptic sub-problem and potential vorticity equations in the QG flow solvers. After solving the elliptic sub-problem on a coarsened grid, an interpolation scheme is used to obtain the fine data for subsequent time stepping on the full grid. The potential vorticity field is then updated on the fine grid with savings in computational time due to the reduced number of grid points for the elliptic solver. The method is applied to both single layer barotropic and two-layer stratified QG ocean models for mid-latitude oceanic basins in the beta plane, which are standard prototypes of more realistic ocean dynamics. The method is found to accelerate these computations while retaining the same level of accuracy in the fine-resolution field. A linear acceleration rate is obtained for all the cases we consider due to the efficient linear-cost fast Fourier transform based elliptic solver used. We expect the speed-up of the CGP method to increase dramatically for versions of the method that use other, suboptimal, elliptic solvers, which are generally quadratic cost. It is also demonstrated that numerical oscillations due to lower grid resolutions, in which the Munk scales are not resolved adequately, are effectively eliminated with CGP method.Comment: International Journal for Multiscale Computational Engineering, 2013. arXiv admin note: substantial text overlap with arXiv:1212.0140, arXiv:1212.0922, arXiv:1104.273

The WR-HSS iteration method for a system of linear differential equations and its applications to the unsteady discrete elliptic problem

We consider the numerical method for non-self-adjoint positive definite linear differential equations, and its application to the unsteady discrete elliptic problem, which is derived from spatial discretization of the unsteady elliptic problem with Dirichlet boundary condition. Based on the idea of the alternating direction implicit (ADI) iteration technique and the Hermitian/skew-Hermitian splitting (HSS), we establish a waveform relaxation (WR) iteration method for solving the non-self-adjoint positive definite linear differential equations, called the WR-HSS method. We analyze the convergence property of the WR-HSS method, and prove that the WR-HSS method is unconditionally convergent to the solution of the system of linear differential equations. In addition, we derive the upper bound of the contraction factor of the WR-HSS method in each iteration which is only dependent on the Hermitian part of the corresponding non-self-adjoint positive definite linear differential operator. Finally, the applications of the WR-HSS method to the unsteady discrete elliptic problem demonstrate its effectiveness and the correctness of the theoretical results.Comment: 30 pages, 5 figures, 13 table

A multidomain spectral method for solving elliptic equations

We present a new solver for coupled nonlinear elliptic partial differential equations (PDEs). The solver is based on pseudo-spectral collocation with domain decomposition and can handle one- to three-dimensional problems. It has three distinct features. First, the combined problem of solving the PDE, satisfying the boundary conditions, and matching between different subdomains is cast into one set of equations readily accessible to standard linear and nonlinear solvers. Second, touching as well as overlapping subdomains are supported; both rectangular blocks with Chebyshev basis functions as well as spherical shells with an expansion in spherical harmonics are implemented. Third, the code is very flexible: The domain decomposition as well as the distribution of collocation points in each domain can be chosen at run time, and the solver is easily adaptable to new PDEs. The code has been used to solve the equations of the initial value problem of general relativity and should be useful in many other problems. We compare the new method to finite difference codes and find it superior in both runtime and accuracy, at least for the smooth problems considered here.Comment: 31 pages, 8 figure

Spectral Methods for Numerical Relativity. The Initial Data Problem

Numerical relativity has traditionally been pursued via finite differencing. Here we explore pseudospectral collocation (PSC) as an alternative to finite differencing, focusing particularly on the solution of the Hamiltonian constraint (an elliptic partial differential equation) for a black hole spacetime with angular momentum and for a black hole spacetime superposed with gravitational radiation. In PSC, an approximate solution, generally expressed as a sum over a set of orthogonal basis functions (e.g., Chebyshev polynomials), is substituted into the exact system of equations and the residual minimized. For systems with analytic solutions the approximate solutions converge upon the exact solution exponentially as the number of basis functions is increased. Consequently, PSC has a high computational efficiency: for solutions of even modest accuracy we find that PSC is substantially more efficient, as measured by either execution time or memory required, than finite differencing; furthermore, these savings increase rapidly with increasing accuracy. The solution provided by PSC is an analytic function given everywhere; consequently, no interpolation operators need to be defined to determine the function values at intermediate points and no special arrangements need to be made to evaluate the solution or its derivatives on the boundaries. Since the practice of numerical relativity by finite differencing has been, and continues to be, hampered by both high computational resource demands and the difficulty of formulating acceptable finite difference alternatives to the analytic boundary conditions, PSC should be further pursued as an alternative way of formulating the computational problem of finding numerical solutions to the field equations of general relativity.Comment: 15 pages, 5 figures, revtex, submitted to PR

Matrix-equation-based strategies for convection-diffusion equations

We are interested in the numerical solution of nonsymmetric linear systems arising from the discretization of convection-diffusion partial differential equations with separable coefficients and dominant convection. Preconditioners based on the matrix equation formulation of the problem are proposed, which naturally approximate the original discretized problem. For certain types of convection coefficients, we show that the explicit solution of the matrix equation can effectively replace the linear system solution. Numerical experiments with data stemming from two and three dimensional problems are reported, illustrating the potential of the proposed methodology

An extrapolation cascadic multigrid method combined with a fourth order compact scheme for 3D poisson equation

In this paper, we develop an EXCMG method to solve the three-dimensional Poisson equation on rectangular domains by using the compact finite difference (FD) method with unequal meshsizes in different coordinate directions. The resulting linear system from compact FD discretization is solved by the conjugate gradient (CG) method with a relative residual stopping criterion. By combining the Richardson extrapolation and tri-quartic Lagrange interpolation for the numerical solutions from two-level of grids (current and previous grids), we are able to produce an extremely accurate approximation of the actual numerical solution on the next finer grid, which can greatly reduce the number of relaxation sweeps needed. Additionally, a simple method based on the midpoint extrapolation formula is used for the fourth-order FD solutions on two-level of grids to achieve sixth-order accuracy on the entire fine grid cheaply and directly. The gradient of the numerical solution can also be easily obtained through solving a series of tridiagonal linear systems resulting from the fourth-order compact FD discretizations. Numerical results show that our EXCMG method is much more efficient than the classical V-cycle and W-cycle multigrid methods. Moreover, only few CG iterations are required on the finest grid to achieve full fourth-order accuracy in both the $L^2$-norm and $L^{\infty}$-norm for the solution and its gradient when the exact solution belongs to $C^6$. Finally, numerical result shows that our EXCMG method is still effective when the exact solution has a lower regularity, which widens the scope of applicability of our EXCMG method.Comment: Accepted for publication in Journal of Scientific Computing. arXiv admin note: text overlap with arXiv:1506.0298