Crate ferreus_bbfmm 

Black Box Fast Multipole Method (BBFMM)

This crate is a parallel implementation of the Black Box Fast Multipole Method in Rust.

BBFMM is a kernel-independent, hierarchical algorithm for rapidly evaluating all pairwise interactions in a collection of particles.

While originally developed for radial basis function (RBF) interpolation problems, ferreus_bbfmm has been generalised to support a broad range of FMM use cases where the kernel is smooth (i.e. non-oscillatory).

Features:

• 1D (binary tree), 2D (quadtree) and 3D (octree)
• Optimised low-rank M2L interactions that leverage symmetries and compression
• Both adaptive and non-adaptive tree structures
• Multiple right-hand sides

Example: Fast Matrix-Vector Product

use ferreus_bbfmm::{FmmTree, KernelFunction};
use faer::{Mat, RowRef};
use rand::{Rng, SeedableRng};
use rand::rngs::StdRng;

// Define a kernel that implements the KernelFunction trait
pub struct LinearRbfKernel;

impl KernelFunction for LinearRbfKernel {
    #[inline(always)]
    fn evaluate(&self, target: RowRef<f64>, source: RowRef<f64>) -> f64 {
        let mut dist = 0.0;
        for (t, s) in target.iter().zip(source.iter()) {
            let diff = t - s;
            dist += diff * diff;
        }
        dist.sqrt()
        
        -dist
    }
}
 
let kernel = LinearRbfKernel;
 
// Generate random source points in 3D
let num_points = 10000;
let dim = 3;
let mut rng = StdRng::seed_from_u64(42);
let num_rhs = 2;

let source_points = Mat::from_fn(num_points, dim, |_, _| rng.random_range(-1.0..1.0));
let weights = Mat::from_fn(num_points, num_rhs, |_, _| rng.random_range(0.0..1.0));

// Interpolation order defines the number of Chebyshev nodes in each dimension
// used in the far-field approximation
// A higher interpolation order is more accurate, but takes longer to compute
let interpolation_order = 7;

// Create an adaptive tree
let adaptive_tree = true;

// No need to store empty leaves for fast matrix-vector product
let sparse_tree = true;

// Create a new tree
let mut tree = FmmTree::new(
    source_points.clone(),
    interpolation_order,
    kernel,
    adaptive_tree,
    sparse_tree,
    None,
    None,
);

// Set the weights - this performs an upward pass through the tree
// and sets the multipole coefficients
tree.set_weights(&weights.as_ref());

// Evaluate at the source points
let target_points = source_points.clone();

// Perform a downward pass to set the local coefficients, then perform a leaf evaluation
tree.evaluate(&weights.as_ref(), &target_points);

println!("Evaluated values at source locations: {:?}", tree.target_values);

Example: RBF Evaluator

use ferreus_bbfmm::{FmmTree, FmmParams, M2LCompressionType, KernelFunction};
use faer::{Mat, RowRef};
use rand::{Rng, SeedableRng};
use rand::rngs::StdRng;

// Define a kernel that implements the KernelFunction trait
pub struct LinearRbfKernel;

impl KernelFunction for LinearRbfKernel {
    #[inline(always)]
    fn evaluate(&self, target: RowRef<f64>, source: RowRef<f64>) -> f64 {
        let mut dist = 0.0;
        for (t, s) in target.iter().zip(source.iter()) {
            let diff = t - s;
            dist += diff * diff;
        }
        dist.sqrt()
        
        -dist
    }
}
 
let kernel = LinearRbfKernel;
 
// Generate random source points in 3D
let num_points = 10000;
let dim = 3;
let mut rng = StdRng::seed_from_u64(42);
let num_rhs = 1;

let source_points = Mat::from_fn(num_points, dim, |_, _| rng.random_range(-1.0..1.0));
let weights = Mat::from_fn(num_points, num_rhs, |_, _| rng.random_range(0.0..1.0));

let interpolation_order = 7;

// Creating an adaptive tree for the evaluator uses less memory
let adaptive_tree = true;

// Store empty leaves for general RBF evaluation
let sparse_tree = false;

// For the evaluator we may wish to evaluate over a larger region than the source points cover
let extents = vec![-2.0, -2.0, -2.0, 2.0, 2.0, 2.0f64];

// Optionally define some tuning parameters
let params = FmmParams{
    max_points_per_cell: 256,
    compression_type: M2LCompressionType::ACA,
    epsilon: 10f64.powi(-(interpolation_order as i32)),
    eval_chunk_size: 1024,
};
 
// Create a new tree
let mut tree = FmmTree::new(
    source_points.clone(),
    interpolation_order,
    kernel,
    adaptive_tree,
    sparse_tree,
    Some(extents),
    Some(params),
);

// Set the weights - this performs an upward pass through the tree
// and sets the multipole coefficients
tree.set_weights(&weights.as_ref());

// For implicit modelling where a 'surface following' method of generating an isosurface
// is used, the evaluator may be called many times. In this case it's more efficient to
// perform a single downward pass to set all the local coefficients, then call the evaluator
// on the relevant leaves for each evaluation
tree.set_local_coefficients(&weights.as_ref());

// Create some arbritrary target points
let num_target_points = 100;
let target_points = Mat::from_fn(num_target_points, dim, |_, _| rng.random_range(-2.0..2.0));

// Perform a leaf evaluation
tree.evaluate_leaves(&weights.as_ref(), &target_points);

println!("Evaluated values at target locations: {:?}", tree.target_values);

// Create some more target points
let num_target_points = 1000;
let target_points = Mat::from_fn(num_target_points, dim, |_, _| rng.random_range(-2.0..2.0));

// Perform another leaf evaluation
tree.evaluate_leaves(&weights.as_ref(), &target_points);

println!("Evaluated values at target locations: {:?}", tree.target_values);

References

1. Fong, W., & Darve, E. (2009). The black-box fast multipole method. Journal of Computational Physics, 228(23), 8712–8725.

2. Messner, M., Bramas, B., Coulaud, O., & Darve, E. (2012). Optimized M2L kernels for the Chebyshev interpolation-based fast multipole method.

3. Pouransari, H., & Darve, E. (2015). Optimizing the adaptive fast multipole method for fractal sets. SIAM Journal on Scientific Computing, 37, A1040–A1066.

Structs

FmmParams

Optional parameters for tuning the FMM performance.

FmmTree

A Fast Multipole Method (FMM) tree that organises source points into a hierarchical spatial structure to accelerate kernel summation tasks.

Enums

FmmError

Errors that can occur during FMM tree operations.

M2LCompressionType

Supported M2L operator compression methods.

Traits

KernelFunction

Evaluates a kernel function between a target and source point.