numerics

C++23 numerical routines for dense and structured linear algebra, ODE/PDE integration, Fourier transforms, quadrature, statistics, Monte Carlo sampling, and particle geometry.

The package is a solver suite. Operators provide a common argument form for Krylov and action routines, but the library also includes direct solvers, eigenvalue routines, SVD, FFTs, finite-difference routines, stochastic kernels, and optional BLAS/LAPACK/LAPACKE, FFTW3, OpenMP, CUDA, and MPI backends.

Examples

• PDE Examples heat equation, backward Euler, ADI diffusion
• Poisson Solver Example Dirichlet Poisson solve by DST-I
• FFT Examples real FFT, reusable plans, backend selection
• Operator Examples and Concepts dense, sparse, and matrix-free Krylov solves
• Field Examples 3D scalar/vector fields and Poisson utilities
• Stencil Operator Examples matrix-free stencil operators
• SPH Kernels SPH density, pressure, and viscosity kernels
• PBCLattice2D periodic lattice indexing
• Connected Components component labeling on grid data
• Boltzmann Acceptance Table Metropolis acceptance tables

Implementation Notes

• Cache Blocking Implementation Note cache-blocked dense matmul
• Register Blocking Implementation Note register-blocked dense matmul
• SIMD Kernel Implementation Note SIMD dense kernels and FFT backend selection
• Parallel, GPU, and MPI Implementation Note OpenMP, CUDA, MPI, and backend boundaries
• Performance and Report Workflow benchmark/report workflow and kernel files
• Application Example: Projection Step matrix-free projection solve pattern

Dense Linear Algebra

#include <numerics.hpp>
 
num::Matrix A(2, 2, 0.0);
A(0, 0) = 4.0; A(0, 1) = 1.0;
A(1, 0) = 1.0; A(1, 1) = 3.0;
 
num::Vector b{1.0, 2.0};
auto F = num::lu(A);
 
num::Vector x;
num::lu_solve(F, b, x);

Iterative Solvers

num::operators::DenseOp Aop(A);
num::Vector x(2, 0.0);
 
num::SolverResult info = num::cg(Aop, b, x, 1e-10, 100);

Matrix-Free Operators

auto Lop = num::operators::make_op(
    [N](const num::Vector& u, num::Vector& Lu) {
        apply_negative_laplacian(u, Lu, N);
    },
    N * N);
 
num::SolverResult info = num::gmres(Lop, rhs, sol);

ODE Integration

auto rhs = [](double, const num::Vector& y, num::Vector& dy) {
    dy[0] = y[1];
    dy[1] = -y[0];
};
 
auto result = num::solve(
    num::ODEProblem{rhs, {1.0, 0.0}, 0.0, 20.0},
    num::RK45{.rtol = 1e-8, .atol = 1e-10});

Build

cmake -B build -DNUMERICS_BUILD_TESTS=ON
cmake --build build -j$(nproc)
ctest --test-dir build --output-on-failure

CMake Use

find_package(numerics REQUIRED)
target_link_libraries(my_program PRIVATE numerics::numerics)