A C++ scientific computing library: linear algebra, ODE/PDE solvers, spectral methods, statistical mechanics, and gnuplot-based plotting – all under a single #include "numerics.hpp".

Example: Lorenz attractor

The Lorenz system is integrated with num::solve + num::RK45 (adaptive Dormand-Prince) and the trajectory is plotted with num::plt – the full program is ~30 lines (examples/lorenz.cpp):

#include "numerics.hpp"
 
int main() {
    const double sigma = 10.0, rho = 28.0, beta = 8.0 / 3.0;
 
    auto lorenz = [&](double, const num::Vector& s, num::Vector& ds) {
        ds[0] = sigma * (s[1] - s[0]);
        ds[1] = s[0] * (rho - s[2]) - s[1];
        ds[2] = s[0] * s[1] - beta * s[2];
    };
 
    num::Series xz;
    auto observer = [&](double, const num::Vector& y) {
        xz.emplace_back(y[0], y[2]);
    };
 
    num::solve(
        num::ODEProblem{lorenz, {1.0, 0.0, 0.0}, 0.0, 50.0},
        num::RK45{.rtol = 1e-8, .atol = 1e-10, .max_steps = 2000000},
        observer);
 
    num::plt::plot(xz);
    num::plt::title("Lorenz attractor (sigma=10, rho=28, beta=8/3)");
    num::plt::xlabel("x");
    num::plt::ylabel("z");
    num::plt::savefig("lorenz.png");
}

Lorenz attractor phase portrait (x vs z)

Example: 2D heat equation

2D Heat Equation

Implicit backward Euler via sparse CG. Grid2D carries geometry; ScalarField2D carries data. The sparse system A = I - coeff*L is assembled once and reused at every step:

num::Grid2D       grid{N, h};
num::SparseMatrix A      = num::pde::backward_euler_matrix(grid, coeff);
num::LinearSolver solver = num::make_cg_solver(A);
 
num::ScalarField2D u(grid, [=](double x, double y) {
    return num::gaussian2d(x, y, 0.5, 0.5, sigma);
});
num::solve(u, num::BackwardEuler{.solver=solver, .dt=dt, .nstep=nstep});

Initial Gaussian and diffused distribution at t=T_end.

Library Modules

kernel – Performance substrate

The inner-loop layer shared by all solvers. Three tiers:

Tier 1 (kernel::raw): always-inline dot/axpy/scale, BLAS-dispatched
Tier 2 (kernel::array, kernel::reduce): tag-dispatched ops with seq_t/par_t compile-time policy tags; zero runtime overhead
Tier 3 (kernel::subspace): LinearOp interface, mgs_orthogonalize, arnoldi_step — the shared inner loops for GMRES, Lanczos, and expv

core – Vectors, matrices, and high-level dispatch

num::Vector, num::Matrix, num::SparseMatrix, num::BandedMatrix
num::Backend – seq, blocked, simd, blas, omp, gpu (for matmul/matvec)
num::matmul, num::matmul_blocked, num::matmul_register_blocked, num::matmul_simd
num::matvec, num::matadd
num::dot, num::axpy, num::scale, num::norm
High-Performance Computing – cache blocking, SIMD, BLAS, backend selection
Parallel and Distributed Computing – OpenMP and CUDA backends

factorization – Direct linear solvers

num::lu, num::lu_solve, num::lu_det, num::lu_inv
num::qr, num::qr_solve
num::thomas – O(n) tridiagonal (Thomas algorithm)

solvers – Iterative solvers

num::cg, num::cg_matfree – conjugate gradient
num::gmres – restarted GMRES (Krylov)
num::jacobi, num::gauss_seidel – stationary iterative

eigen – Eigenvalue methods

num::power_iteration, num::inverse_iteration, num::rayleigh_iteration
num::lanczos – Krylov-Lanczos with Ritz extraction
num::eig_sym – dense symmetric eigensolver (Jacobi sweeps)

svd – Singular value decomposition

num::svd – full SVD
num::svd_truncated – randomized truncated SVD

spectral – Fourier transforms

Fast Fourier Transform
num::spectral::fft, num::spectral::ifft – complex DFT / inverse DFT
num::spectral::rfft, num::spectral::irfft – real-to-complex and inverse
num::spectral::FFTPlan – precomputed plan for repeated transforms
FFTBackend::seq (Cooley-Tukey radix-2), FFTBackend::fftw (FFTW3, optional)

analysis – Quadrature and root finding

Quadrature: num::trapz, num::simpson, num::gauss_legendre, num::adaptive_simpson, num::romberg
Root finding: num::bisection, num::newton, num::secant, num::brent

stats – Simulation observables

num::RunningStats – Welford online mean + variance, O(1) memory
num::Histogram – fixed-bin histogram; reweighting for WHAM analysis
num::autocorr_time – integrated autocorrelation time (Madras-Sokal windowing)

markov – Markov chain Monte Carlo

num::markov::metropolis_sweep – Metropolis sweep; computes exp(-beta*dE) internally
num::markov::metropolis_sweep_prob – caller-supplied acceptance probability (e.g. Boltzmann table)
num::markov::umbrella_sweep, num::markov::umbrella_sweep_prob – umbrella sampling with save/restore
num::markov::UmbrellaWindow, MetropolisStats, UmbrellaStats
num::markov::make_seeded_rng – hardware-entropy RNG seeding

fields – Grid geometry and field storage

num::Grid2D – geometry only: N, h, coordinate accessors, flat index
num::Grid3D – 3D analogue with nx, ny, nz
num::ScalarField2D – owns Grid2D + data vector; callable constructor fills from f(x,y)
num::ScalarField3D, num::VectorField3D – 3D field types with trilinear sampling

solve – Unified solver dispatcher

num::solve(problem, algorithm) – C++20 concept-dispatched entry point
num::ODEProblem{f, u0, t0, tf} – explicit ODE formulation
num::MCMCProblem{accept_prob, propose, n_sites} – MCMC formulation
Algorithm tags: num::Euler, num::RK4, num::RK45, num::BackwardEuler, num::Metropolis

sparse and banded – Structured sparse formats

num::SparseMatrix – CSR format; num::sparse_matvec
num::BandedMatrix – band-storage; num::banded_matvec, num::banded_solve
num::LinearSolver – universal callable: (rhs, x) -> SolverResult

meshfree – Particle-based spatial structures

num::CellList2D, num::CellList3D – linked-cell neighbour search
num::VerletList – Verlet neighbour list with skin radius

plot – Gnuplot pipe wrapper

num::Gnuplot – popen wrapper: operator<< for commands, send1d() for inline data
num::Series – std::vector<std::pair<double,double>> for (x, y) data
num::apply_siam_style(gp) – clean SIAM-style theme
num::save_png(gp, file, w, h) – redirect output to PNG
num::set_loglog(gp), num::set_logx(gp) – axis scaling helpers

Physics Simulations

Raylib-rendered simulations built on numerics. Each runs as a batch renderer: simulate all frames, export PNGs, composite with make_video.sh.

Monte Carlo

2D Ising Model

page_app_ising

Metropolis dynamics and umbrella-sampled nucleation on a 300x300 lattice. Reproduces Brendel et al. (2005) nucleation experiments.

Library feature	Role
`num::markov::metropolis_sweep_prob`	Sweep hot path with precomputed Boltzmann table
`num::markov::umbrella_sweep_prob`	Per-sweep rejection / save-restore for window sampling
`num::SparseMatrix` + `sparse_matvec`	Neighbour-sum matrix for energy observable (SIMD path)
`num::RunningStats`, `num::Histogram`	Online mean, variance, nucleus-size distribution per window
`num::newton`	Mean-field self-consistency equation solver

Fluid Dynamics

2D SPH Fluid Simulation

page_app_fluid

Weakly compressible SPH with heat transport, rigid bodies, and particle injection. Tait EOS, cubic-spline kernel, Morris viscosity Laplacian.

Library feature	Role
`num::CellList2D`	O(1) neighbour queries; Newton-3rd-law pair traversal
`num::Backend` dispatch	`seq` (pair loop) <-> `omp` (per-particle) selectable at runtime
`num::Vector`	Flat particle attribute arrays

3D SPH Fluid Simulation

page_app_fluid3d

3D WCSPH with opposing hose jets, heat transport, and a free-orbit camera whose orientation rotates the gravity vector.

Library feature	Role
`num::CellList3D`	3D counting-sort spatial hash; 13-stencil Newton-3rd-law pairs
`num::Backend` dispatch	`seq` (pair traversal) <-> `omp` (per-particle)
`num::Vector`	Flat particle attribute arrays

2D Incompressible Navier-Stokes

page_app_ns

Chorin projection on a staggered MAC grid; semi-Lagrangian advection; matrix-free CG pressure solve. Kelvin-Helmholtz double shear layer initial condition.

Library feature	Role
`num::cg_matfree`	Matrix-free CG for -grad^2p = r; warm-start from previous pressure field
`num::dot`	Inner products inside CG (dispatched to best available backend)
`num::Vector`	Velocity faces u, v; pressure p; intermediates u, v, rhs

Quantum Mechanics

2D Time-Dependent Schrodinger Equation

page_app_tdse

Strang operator splitting with Crank-Nicolson kinetic sweeps; Thomas algorithm for O(N) tridiagonal solves; Lanczos eigendecomposition; five interchangeable potentials.

Library feature	Role
`num::thomas`	O(N) complex tridiagonal solve per row/column (kinetic sweeps)
`num::lanczos`	Krylov subspace for lowest eigenmodes
`num::brent`	Bessel zero-finding for exact circular-well energies
`num::gauss_legendre`	Norm and energy quadrature

Electromagnetism

Electromagnetic Field Demo

page_app_em

DC current flow + magnetostatics on a 32^3 voxel grid. Four Poisson problems (one electric, three magnetic) solved with matrix-free CG; interactive magnetic dipole.

Library feature	Role
`num::cg_matfree`	All four Poisson solves (variable-coefficient electric + 3x magnetic)
`num::Grid3D`	3D scalar field storage for phi, A, J, B components
`num::Vector`	Flattened field arrays passed to the solver

Status Report

Run cmake --build build --target report to generate output/REPORT.md – a full snapshot of the current build:

Build environment – compiler, flags, detected backends (BLAS, FFTW3, OpenMP, CUDA, MPI)
Test summary – pass/fail counts per suite
Benchmark tables – throughput and timing for every module (linalg, FFT, banded, analysis)
Plots – PNG throughput curves (seq vs fftw, size scaling) in output/plots/

Sections whose backend was not found at configure time are left as placeholders, so the report always reflects exactly what is installed.

Performance Notes

Notes on hardware-aware implementation: cache blocking, register tiling, and SIMD vectorization applied to the matrix operations in this library.

Algorithm Theory

For derivations, error bounds, and proofs of the algorithms implemented here (LU/QR factorization, CG, GMRES, Lanczos, eigenvalue methods, quadrature, root finding, ODE integrators), see the companion notes:

Scientific Computing — adityadendukuri.github.io/cs111-scientific-computing-notes