la-stack

Fast, stack-allocated linear algebra for fixed dimensions in Rust.

This crate grew from the need to support delaunay with fast, stack-allocated linear algebra primitives and algorithms while keeping the API intentionally small and explicit.

📐 Introduction

la-stack provides a handful of const-generic, stack-backed building blocks:

• Vector<const D: usize> for fixed-length vectors ([f64; D] today)
• Matrix<const D: usize> for fixed-size square matrices ([[f64; D]; D] today)
• Lu<const D: usize> for LU factorization with partial pivoting (solve + det)
• Ldlt<const D: usize> for LDLT factorization without pivoting (solve + det; symmetric SPD/PSD)

✨ Design goals

• ✅ Copy types where possible
• ✅ Const-generic dimensions (no dynamic sizes)
• ✅ Explicit algorithms (LU, solve, determinant)
• ✅ No runtime dependencies (dev-dependencies are for contributors only)
• ✅ Stack storage only (no heap allocation in core types)
• ✅ unsafe forbidden

🚫 Anti-goals

• Bare-metal performance: see blas-src, lapack-src, openblas-src
• Comprehensive: use nalgebra if you need a full-featured library
• Large matrices/dimensions with parallelism: use faer if you need this

🔢 Scalar types

Today, the core types are implemented for f64. The intent is to support f32 and f64 (and f128 if/when Rust gains a stable primitive for it). Longer term, we may add optional arbitrary-precision support (e.g. via rug) depending on performance.

🚀 Quickstart

Add this to your Cargo.toml:

[dependencies]
la-stack = "0.1"

Solve a 5×5 system via LU:

use la_stack::prelude::*;

// This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
// A = J - I: zeros on diagonal, ones elsewhere.
let a = Matrix::<5>::from_rows([
    [0.0, 1.0, 1.0, 1.0, 1.0],
    [1.0, 0.0, 1.0, 1.0, 1.0],
    [1.0, 1.0, 0.0, 1.0, 1.0],
    [1.0, 1.0, 1.0, 0.0, 1.0],
    [1.0, 1.0, 1.0, 1.0, 0.0],
]);

let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]);

let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap();
let x = lu.solve_vec(b).unwrap().into_array();

// Floating-point rounding is expected; compare with a tolerance.
let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
for (x_i, e_i) in x.iter().zip(expected.iter()) {
    assert!((*x_i - *e_i).abs() <= 1e-12);
}

Compute a determinant for a symmetric SPD matrix via LDLT (no pivoting).

For symmetric positive-definite matrices, LDL^T is essentially a square-root-free form of the Cholesky decomposition (you can recover a Cholesky factor by absorbing sqrt(D) into L):

use la_stack::prelude::*;

// This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
let a = Matrix::<5>::from_rows([
    [1.0, 1.0, 0.0, 0.0, 0.0],
    [1.0, 2.0, 1.0, 0.0, 0.0],
    [0.0, 1.0, 2.0, 1.0, 0.0],
    [0.0, 0.0, 1.0, 2.0, 1.0],
    [0.0, 0.0, 0.0, 1.0, 2.0],
]);

let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det();
assert!((det - 1.0).abs() <= 1e-12);

🧩 API at a glance

Type	Storage	Purpose	Key methods
Vector<D>	[f64; D]	Fixed-length vector	new, zero, dot, norm2_sq
Matrix<D>	[[f64; D]; D]	Fixed-size square matrix	from_rows, zero, identity, lu, ldlt, det
Lu<D>	Matrix<D> + pivot array	Factorization for solves/det	solve_vec, det
Ldlt<D>	Matrix<D>	Factorization for symmetric SPD/PSD solves/det	solve_vec, det

Storage shown above reflects the current f64 implementation.

📋 Examples

The examples/ directory contains small, runnable programs:

just examples
# or:
cargo run --example solve_5x5
cargo run --example det_5x5

🤝 Contributing

A short contributor workflow:

cargo install just
just setup        # install/verify dev tools + sync Python deps
just check        # lint/validate (non-mutating)
just fix          # apply auto-fixes (mutating)
just ci           # lint + tests + examples + bench compile

For the full set of developer commands, see just --list and AGENTS.md.

📝 Citation

If you use this library in academic work, please cite it using CITATION.cff (or GitHub's "Cite this repository" feature). A Zenodo DOI will be added for tagged releases.

📚 References

For canonical references to LU / LDL^T algorithms used by this crate, see REFERENCES.md.

📊 Benchmarks (vs nalgebra/faer)

Raw data: docs/assets/bench/vs_linalg_lu_solve_median.csv

Summary (median time; lower is better). The “la-stack vs nalgebra/faer” columns show the % time reduction relative to each baseline (positive = la-stack faster):

D	la-stack median (ns)	nalgebra median (ns)	faer median (ns)	la-stack vs nalgebra	la-stack vs faer
2	2.044	18.266	164.197	+88.8%	+98.8%
3	13.465	23.723	214.231	+43.2%	+93.7%
4	27.774	53.689	238.476	+48.3%	+88.4%
5	46.982	71.070	301.806	+33.9%	+84.4%
8	138.664	177.992	388.146	+22.1%	+64.3%
16	629.219	589.141	915.520	-6.8%	+31.3%
32	2,669.149	2,484.327	2,937.819	-7.4%	+9.1%
64	16,673.839	14,833.982	12,528.617	-12.4%	-33.1%

📄 License

BSD 3-Clause License. See LICENSE.