MTL5 — Matrix Template Library 5

C++20 header-only linear algebra library for mixed-precision algorithm design and optimization.

MTL5 is a modernized successor to MTL4, preserving the proven architecture for custom number types (posits, LNS, etc.) while leveraging C++20 features for cleaner, safer, and more maintainable code. Its distinguishing feature is a mixed-precision accumulator model — store narrow, accumulate wide, round out once — carried through the dense kernels, iterative solvers, and native sparse factorizations.

• Header-only — no build step; drop include/ on your path and #include <mtl/mtl.hpp>
• Zero Boost dependency — pure C++20 (concepts, constexpr, std::span, ranges)
• Mixed precision throughout — a shared accumulator_traits policy expresses the three precisions of a kernel (element, accumulate, result); e.g. mult<float>(A_bf16, B_bf16, C_bf16) accumulates in fp32 and stores bf16 once. The policy threads through the dense/sparse kernels, every Krylov solver, and the smoothers, plus dense LU iterative refinement for low-precision-factor + high-precision-residual solves
• Custom arithmetic types — designed for posits, LNS, and other Universal number types; the accumulator is a pluggable policy, so an exact quire super-accumulator drops in from the Universal pairing (MTL5 itself stays library-free)
• Dense & sparse — CSR/CSC, COO, ELLPACK, block-diagonal, plus dense row/column-major
• Complete BLAS surface — L1/L2/L3, generic over any type, auto-dispatching to external BLAS/LAPACK for dense float/double
• Sparse direct solvers — native Cholesky/LDLᵀ/LU/QR, supernodal LU/LDLᵀ, KLU, fill-reducing orderings, plus wrappers for SuiteSparse/SuperLU
• Iterative solvers — a full Krylov suite, preconditioners, smoothers, and multigrid
• Eigen & SVD — dense (symmetric + general with eigenvectors), matrix-free iterative (Lanczos/Arnoldi/power), and sparse shift-invert
• On-node threading — a dependency-free thread pool (no OpenMP/TBB); every substantial dense and sparse kernel is threaded — BLAS L1/2/3 (incl. multi-loop GEMM), dense factorizations (LU/Cholesky/QR), sparse direct solvers (level-scheduled), and element-wise sweeps — and bit-identical across thread counts, serial by default
• Optional SIMD — Google Highway-backed batch kernels with narrow→wide widening on load

• CMake 3.22+
• C++20 compiler (GCC 11+, Clang 14+, MSVC 2022+, Apple Clang)
• Catch2 v3 (fetched automatically for tests)

Portable across Linux (x64/ARM64), macOS (ARM64), and Windows (MSVC / Clang-CL).

# Clone
git clone https://github.com/stillwater-sc/mtl5.git
cd mtl5

# Build (uses dev preset: Debug, tests ON)
cmake --preset dev
cmake --build build -j$(nproc)

# Run tests
ctest --test-dir build

# Run an example
./build/examples/phase01_core_types/...

Because MTL5 is header-only, using it in your own project needs no build of MTL5 itself:

#include <mtl/mtl.hpp>

int main() {
    mtl::mat::dense2D<double> A(3, 3);
    mtl::vec::dense_vector<double> x(3), b(3);
    // ... fill A, b ...
    mtl::lu_apply(A, x, b);   // factor A and solve A*x = b (A modified in place)
}

• Dense — dense2D (row- or column-major), dense_vector, strided_vector_ref, unit_vector
• Sparse matrices — compressed2D (CSR/CSC), coordinate2D (COO), ell_matrix (ELLPACK), block_diagonal2D, identity2D, permutation_matrix, plus sparse_vector
• Expression templates & views — lazy element-wise expressions (CRTP), transposed / sub-matrix / row / column views
• Mathematical tensors — tensor<T, Rank, Dim> (stack-allocated, compile-time rank/dimension), symmetric/antisymmetric tensors, metric and index helpers
• NumPy-style arrays — array::ndarray<T, N, Order>: static rank, runtime shape/strides, owning or view, C- or F-order, slicing, broadcasting, and interop

• BLAS Level 1 — dot, dot_real, axpy, scal, and the norm family (two_norm/nrm2, one_norm, inf_norm, frobenius_norm)
• BLAS Level 2 — gemv, ger (rank-1 update), symv, trmv, trsv
• BLAS Level 3 — gemm/mult, trmm, trsm, symm, syrk, syr2k
• Factorizations — LU, QR, LQ, Cholesky, LDLᵀ (and Bunch–Kaufman ldlt_bk), SVD, Hessenberg reduction, Householder & Givens
• Eigenvalues — symmetric solver and general (Francis implicit double-shift QR); eigen returns eigenvalues and eigenvectors
• Solvers & inverses — triangular solves, inv, sparse_solve
• Mixed-precision iterative refinement — lu_iterative_refine<Working>(A, b, x) factors once in a low Working precision and corrects with a higher-precision residual (best-iterate return, patience, optional scale-and-round), plus normwise_backward_error as the quality/termination metric
• Structure & utility — trace, diagonal, kron (Kronecker product), projection, reorder, trans, conj, real/imag, product, sum, min/max, fill, random
• Transcendentals — the full element-wise set: trig, inverse trig, hyperbolic, exp/log/pow/sqrt/cbrt, erf/erfc, rounding
• Property predicates — is_symmetric/is_hermitian, is_spd/is_positive_definite, is_singular/is_invertible, determinant, condition_number/rcond, numerical_rank/nullity, spectral_radius, inertia, is_orthogonal/is_unitary/is_normal, and many structural/vector checks

• math::accumulator_traits<Acc, Value> — a cross-cutting policy expressing element (storage), accumulator (compute), and result (serialize) precisions; the accumulate→output conversion is fused into the final store. It abstracts a sum-of-products reduction so a kernel writes one loop regardless of how terms combine, and covers three configurations: plain acc += product (the default primary template), fused multiply-add via math::fma_accumulator<T> (one rounding per term, no intermediate product round), and a caller-supplied super-accumulator (e.g. an exact quire) for single-rounding dot products
• convert — standalone element-wise re-quantization (distinct from the fused epilogue)
• Threaded through the dense kernels (dot, gemm/mult, gemv, sum-of-squares norms), the sparse factorizations, every Krylov solver, and the stationary smoothers — each takes an optional Accumulator template parameter; default (Accumulator = void) is byte-identical to the plain path

Google Highway-backed batch, algorithm, and blocking layers with widening-on-load (load_widen): the micro-kernels promote narrow operands into wide accumulator registers, giving large speedups for mixed-precision GEMM and dot. Enabled with -DMTL5_WITH_HIGHWAY=ON.

A persistent detail::thread_pool built on the C++ standard concurrency runtime (no OpenMP/TBB) with parallel_for and parallel_reduce. Every substantial dense and sparse kernel is threaded and bit-identical across thread counts — the parallel result matches the serial one exactly (==), not just to a tolerance:

• Level 1/2/3 BLAS — axpy/scal, dot/nrm2, GEMV, and a BLIS-style blocked GEMM with multi-loop (2D jc × ic) parallelism that scales tall, wide, and square shapes
• Dense factorizations — LU, Cholesky, and Householder QR (parallel trailing/column updates and reflector application)
• Sparse direct solvers — the triangular solves of sparse Cholesky, LDLᵀ, LU, and their supernodal variants, parallelized by level scheduling (value-agnostic schedules that survive same-pattern in-place refactorization)
• Sparse SpMV and the element-wise expression sweeps behind y = a + b, C += A + B, etc.

Bit-identity is preserved by construction: reductions are recast to one-write-per-output form, the sparse solves replay the serial accumulation order, and each GEMM C-block gets the same FMAs in the same order regardless of grid shape. Threading is off by default (serial, zero overhead); MTL5_NUM_THREADS sizes the pool. Iterative and eigen solvers inherit the SpMV/L1 threading with no solver-code changes. See docs/design/parallelization-patterns-and-pitfalls.md for the design rationale.

• Krylov — cg, bicg, bicgstab, bicgstab_ell, cgs, gmres, idr_s, minres, qmr, tfqmr — all accumulator-aware (route their dot/mult through accumulator_traits)
• Preconditioners (itl::pc) — identity, diagonal, block_diagonal, ic_0, ildl, ilu_0, ilut, ssor
• Smoothers (itl::smoother) — jacobi; gauss_seidel with backward_gauss_seidel and symmetric_gauss_seidel (SGS); sor with backward_sor and symmetric_sor (SSOR) — forward/backward/symmetric sweep directions, all accumulator-aware
• Multigrid (itl::mg) — geometric multigrid with prolongation/restriction
• Iterative eigensolvers (itl::eigen) — power_iteration, lanczos, arnoldi, matrix-free through the LinearOperator concept

• Fill-reducing orderings — RCM, AMD, COLAMD, minimum-degree, Dulmage–Mendelsohn
• Symbolic analysis — elimination tree (O(nnz)), column elimination tree, postorder, supernode partitioning
• Numeric factorization — sparse_cholesky (LLᵀ), sparse_ldlt, sparse_lu (threshold partial pivoting), sparse_qr (least squares), native supernodal LU/LDLᵀ, and native KLU; all generic over the mixed-precision accumulator
• Refactorization — reuse a prior symbolic structure + pivot sequence to recompute same-pattern matrices ~2–3× faster (the SPICE-transient path)
• Iterative refinement — Universal-free, templated residual precision, scaled variant for narrow-exponent low-precision factors
• Sparse eigen — largest-magnitude Arnoldi and shift-invert for eigenpairs near a target

Auto-dispatching bindings that engage when the type qualifies (dense column-major float/double) and the library is present, otherwise falling back to the in-house path: BLAS (L1/L2/L3), LAPACK (factorizations, syev/geev eigensolvers), UMFPACK, SuperLU, KLU, CHOLMOD, SPQR. Any non-default accumulator forces the native kernel, since external BLAS cannot honor a custom accumulator.

• Matrix Market reader/writer, with transparent gzip (.mtx.gz) when built with zlib
• .el edge-list read/write
• From-first-principles PNG writer and spy sparsity-pattern visualization

• Parametric matrices — classic and random test matrices at arbitrary size with controllable conditioning: hilbert, frank, wilkinson, clement, companion, forsythe, kahan, lehmer, lotkin, minij, moler, pascal, rosser, vandermonde, laplacian, poisson, ones, magic, and randomized randorth, randspd, randsym, randsvd
• Range / spacing vectors — NumPy-style arange, linspace, logspace, geomspace returning dense_vector<T>
• Named test-matrix catalog (mtl::testsuite) — well-known reference matrices (SuiteSparse / textbook problems) with published condition numbers, so studies run on the same named problems: testsuite::by_name("bcsstk01"), testsuite::kappa(name), testsuite::names()

Threading is a runtime setting, not a build option: set the MTL5_NUM_THREADS environment variable (unset or 1 runs the serial paths).

include/mtl/
├── concepts/       # C++20 concepts (Scalar, Matrix, Vector, LinearOperator, ...)
├── tag/            # Compile-time tags (orientation, sparsity, shape, storage)
├── traits/         # Type traits and metafunctions
├── math/           # Algebraic identities + mixed-precision accumulator policy
├── detail/         # Internal: memory blocks, GEMM kernels/packing, thread pool
├── simd/           # Highway-backed batch/blocking SIMD layer
├── mat/            # Matrix types, expressions, views
├── vec/            # Vector types and expressions
├── tensor/         # Mathematical tensors (compile-time rank/dimension)
├── array/          # NumPy-style ndarray (slicing, broadcasting, interop)
├── operation/      # Free-function operations (BLAS, decompositions, eigen, predicates)
├── functor/        # Scalar and typed functors
├── recursion/      # Block-recursive infrastructure
├── generators/     # Test matrix generators
├── io/             # Matrix Market, edge-list, PNG, spy
├── itl/            # Iterative solvers, preconditioners, smoothers, multigrid, eigen
├── sparse/         # Sparse direct solvers: orderings, analysis, factorization
└── interface/      # Optional BLAS/LAPACK/UMFPACK/SuperLU/KLU/CHOLMOD/SPQR bindings

The examples/ directory contains a phased tour (phase01–phase15) from core types through iterative solvers, sparse assembly, decompositions, eigen/SVD, expression templates, I/O, and sparse direct solvers, plus applied demos (e.g. an unscented Kalman filter).

cmake --preset release
cmake --build build-release
cmake --install build-release --prefix /usr/local

Then in your project:

find_package(MTL5 REQUIRED)
target_link_libraries(myapp PRIVATE MTL5::mtl5)

Full documentation — architecture, algorithm write-ups (mixed-precision kernels, on-node threading, eigenvalues, measuring solver accuracy), and a Doxygen C++ API reference — is published from the docs/ tree to the project's GitHub Pages site.

MIT License — see LICENSE.

MTL5 builds on the foundational work of MTL4 by Peter Gottschling and the Simunova team.