direct_sum.hpp File Reference

The dense direct sum of a pairwise interaction: for every target, the sum over every source. More...

Namespaces
namespace	mundy

Functions
template<class ExecSpace, class Interaction, class Accumulate>
void	mundy::direct_sum (const ExecSpace &space, const size_t num_targets, const size_t num_sources, const Interaction &interaction, const Accumulate &accumulate)
	accumulate(t, sum over s of interaction(t, s)) for every target t, split the best way for space.

Detailed Description

direct_sum is the exact O(targets x sources) evaluation that a fast multipole method approximates, for any dense long-range kernel (Stokeslet, RPY, double layer, Coulomb, ...) on one process. The caller supplies two functions: interaction(t, s), the contribution of source s to target t as a scalar, a mundy::Vector, or a mundy::Matrix; and accumulate(t, sum), which receives each target's total exactly once and decides where it goes.

How the work is split depends on the execution space:

• On the host, each thread owns a panel of 4 targets and sweeps the sources once for the whole panel, so each source is loaded once per panel and the arithmetic vectorizes across the panel's targets. Each target adds its sources in order 0, 1, 2, ..., so the result does not depend on the number of threads, up to FMA contraction.
• On a device, each target gets a warp: its 32 vector lanes split the sources and one reduction combines them. Splitting every target's sum keeps a GPU busy from about a thousand targets, where one thread per target needs about a hundred thousand. The result is the same from run to run, but it adds in a different order than the host, so the two agree to rounding. With fewer sources than lanes, lanes idle; with many targets and only a few sources (under about 100), one thread per target, is faster.

To vectorize on the host, keep interaction free of branches: select with ?: rather than branch around a division or a square root, e.g. near ? 0 : rsqrt(near ? 1 : r2) instead of if (near) return 0. MundyMath builds with -fno-math-errno, without which every sqrt keeps a branch that blocks vectorization.