Example usage 1

Suppose we want to evaluate the following sum:

\[f(y_i) = \sum_{j=1}^N k(y_i - x_j) \mu_j\]

where \(k\) evaluates the kernel for 3D electrostatic interactions:

\[k(y_i-x_j) = \frac{1}{4\pi \|y_i-x_j\|}\]

In this setting, we can think of the points \(x_j\) as sources, each with strength \(\mu_j\). The points \(y_i\) are targets. The goal is to evaluate \(f(y_i)\) for many target points \(y_i\) efficiently. We will call the number of target points \(M\) and the number of source points \(N\).

First, we need to define the kernel as a function handle:

function evals = kern(src_info, targ_info)
    % evaluate electrostatic kernel between N source points src_info.r and
    % M target points targ_info.r
    %
    % Output shape is (M, N)

    % Shape (M, N)
    rx = src_info.r(1, :) - targ_info.r(1, :).';
    ry = src_info.r(2, :) - targ_info.r(2, :).';
    rz = src_info.r(3, :) - targ_info.r(3, :).';

    dist = sqrt(rx.^2 + ry.^2 + rz.^2);
    k_evals = 1 ./ dist;
    evals = k_evals / (4*pi);
 end

Next, we need to define the source and target points. For simplicity, we will just use random points in a box:

N = 1000; % number of sources
M = 3000; % number of targets
src_info = struct();
targ_info = struct();
src_info.r = rand(3, N);
targ_info.r = rand(3, M);
mu = rand(N, 1);

To evaluate the sum, we first need to construct a regular grid that will be used for the FFT. We need to choose the grid spacing carefully so that sources spread on the regular grid can accurately represent the far-field kernel evaluations. get_grid() handles this for us:

tol = 1e-6;
[grid_info, proxy_info] = get_grid(@kern, src_info, targ_info, tol);

Next, we call get_spread() to compute the matrices which spread the sources and targets to the regular grid.

[A_spread_src, srt_info_src] = get_spread(@kern, [], src_info, ...
                                             grid_info, proxy_info);
[A_spread_targ, srt_info_targ] = get_spread(@kern, [], targ_info, ...
                                             grid_info, proxy_info);

Note that the second argument was left empty as we do not wish to take any derivatives of k.

The matrices A_spread_src and A_spread_targ take care of the far-field interactions, but we need to correct for near-field interactions which must be computed exactly. We do this by calling get_addsub():

A_addsub = get_addsub(@kern, [], grid_info, proxy_info,  ...
                      srt_info_src, srt_info_targ, ...
                      A_spread_src, A_spread_targ);

Finally, we can evaluate the sum by calling get_kernhat() and pcfft_apply():

kern_hat = get_kernhat(@kern, grid_info);
f = pcfft_apply(mu, A_spread_src, A_spread_targ, ...
                A_addsub, kern_hat);

Notice that we didn’t use mu until the final line of this example. If the source points, target points, and kernel are fixed, we can think of all of the steps before the call to pcfft_apply() as a precomputation step. Once the precomputation is done, we can apply the same operator to different source strengths \(\mu\) very quickly by calling pcfft_apply() again with the new source strengths and the same matrices.