Compression methods

Inti.jl wraps several external libraries providing acceleration routines for integral operators. In general, acceleration routines have the signature assemble_*(iop, args...; kwargs...), and take an IntegralOperator as a first argument. They return a new object that represents a compressed version of the operator. The following methods are available:

• assemble_matrix: create a dense Matrix representation of the integral operator. Not really a compression method, but useful for debugging and small problems.
• assemble_hmatrix: assemble a hierarchical matrix representation of the operator using the HMatrices library.
• assemble_fmm: return a LinearMap object that represents the operator using the fast multipole method. This method is powered by the FMM2D, FMMLIB2D and FMM3D libraries, and is only available for certain kernels.

To illustrate the use of compression methods, we will use the following problem as an example. Note that for such a small problem, compression methods are not likely not necessary, but they are useful for larger problems.

using Inti
using LinearAlgebra
# define the quadrature
geo = Inti.GeometricEntity("ellipsoid")
Ω = Inti.Domain(geo)
Γ = Inti.boundary(Ω)
Q = Inti.Quadrature(Γ; meshsize = 0.4, qorder = 5)
# create the operator
op = Inti.Helmholtz(; dim = 3, k = 2π)
K = Inti.SingleLayerKernel(op)
Sop = Inti.IntegralOperator(K, Q, Q)
x = rand(eltype(Sop), length(Q))
rtol = 1e-8

In what follows we compress Sop using the different methods available.

Dense matrix

assemble_matrix(iop::IntegralOperator; threads = true)

Typically used for small problems, the dense matrix representation converts the IntegralOperator into a Matrix object. The underlying type of the Matrix is determined by the eltype of the IntegralOperator, and depends on the inferred type of the kernel. Here is how assemble_matrix can be used:

Smat = Inti.assemble_matrix(Sop; threads=true)
er = norm(Sop * x - Smat * x, Inf) / norm(Sop * x, Inf)
println("Forward map error: $er")

Forward map error: 6.037452551574307e-15

Since the returned object is plain Julia Matrix, it can be used with any of the linear algebra routines available in Julia (e.g. \, lu, qr, *, etc.)

Hierarchical matrix

assemble_hmatrix(iop[; atol, rank, rtol, eta])

The hierarchical matrix representation is a compressed representation of the underlying operator; as such, it takes a tolerance parameter that determines the relative error of the compression. Here is an example of how to use the assemble_hmatrix method to compress the previous problem:

using HMatrices
Shmat = Inti.assemble_hmatrix(Sop; rtol = 1e-8)
er = norm(Smat * x - Shmat * x, Inf) / norm(Smat * x, Inf)
println("Forward map error: $er")

Forward map error: 2.0400023872607663e-8

Note that HMatrices are said to be kernel-independent, meaning that they efficiently compress a wide range of integral operators provided they satisfy a certain asymptotic smoothness criterion (see e.g. [3, 4]).

The HMatrix object can be used to solve linear systems, both iteratively through e.g. GMRES, or directly using an LU factorization.

Fast multipole method

assemble_fmm(iop; atol)

The fast multipole method (FMM) is an acceleration technique based on an analytic multipole expansion of the kernel in the integral operator [5, 6]. It provides a very memory-efficient and fast way to evaluate certain types of integral operators. Here is how assemble_fmm can be used:

using FMM3D
Sfmm = Inti.assemble_fmm(Sop; rtol = 1e-8)
er = norm(Sop * x - Sfmm * x, Inf) / norm(Sop * x, Inf)
println("Forward map error: $er")

Forward map error: 1.0535441704352412e-10

Tips on choosing a compression method

The choice of compression method depends on the problem at hand, as well as on the available hardware. Here is a rough guide on how to choose a compression:

1. For small problems (say less than 5k degrees of freedom), use the dense matrix representation. It is the simplest and most straightforward method, and does not require any additional packages. It is also the most accurate since it does not introduce any approximation errors.
2. If the integral operator is supported by the assemble_fmm, and if an iterative solver is acceptable, use it. The FMM is a very efficient method for certain types of kernels, and can handle problems with up to a few million degrees of freedom on a laptop.
3. If the kernel is not supported by assemble_fmm, if iterative solvers are not an option, or if the system needs solution for many right-hand sides, use the assemble_hmatrix method. It is a very general method that can handle a wide range of kernels, and although assembling the HMatrix can be time and memory consuming (the complexity is still log-linear in the DOFs for many kernels of interest, but the constants can be large), the resulting HMatrix object is very efficient to use. For example, the forward map is usually significantly faster than the one obtained through assemble_fmm.