Changed: Refactors benchmarks to compare against other crates

Renames and updates the benchmarking suite to compare against
multiple linear algebra crates, including faer, in addition to
nalgebra, providing a broader performance context.
Adds faer implementations and benchmarks.
Updates data and plots.

Author: acgetchell

Cargo.lock

@@ -17,12 +17,13 @@ keywords = ["linear-algebra", "geometry", "const-generics"]
 [dev-dependencies]
 approx = "0.5.1"
 criterion = { version = "0.8.1", features = ["html_reports"] }
+faer = { version = "0.23.2", default-features = false, features = ["std", "linalg"] }
 nalgebra = "0.34.1"
 pastey = "0.2.0"
 proptest = "1.9.0"
 
 [[bench]]
-name = "vs_nalgebra"
+name = "vs_linalg"
 harness = false
 
 [lints.rust]

Cargo.toml

@@ -34,8 +34,9 @@ while keeping the API intentionally small and explicit.
 
 ## 🚫 Anti-goals
 
-- Comprehensive: use [`nalgebra`](https://crates.io/crates/nalgebra) if you need a full-featured library
 - Bare-metal performance: see [`blas-src`](https://crates.io/crates/blas-src), [`lapack-src`](https://crates.io/crates/lapack-src), [`openblas-src`](https://crates.io/crates/openblas-src)
+- Comprehensive: use [`nalgebra`](https://crates.io/crates/nalgebra) if you need a full-featured library
+- Large matrices/dimensions with parallelism: use [`faer`](https://crates.io/crates/faer) if you need this
 
 ## 🔢 Scalar types
 
@@ -112,25 +113,25 @@ just commit-check # lint + all tests + examples
 
 For the full set of developer commands, see `just --list` and `WARP.md`.
 
-## 📊 Benchmarks (vs nalgebra)
+## 📊 Benchmarks (vs nalgebra/faer)
 
-![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_nalgebra_lu_solve_median.svg)
+![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_linalg_lu_solve_median.svg)
 
-Raw data: [docs/assets/bench/vs_nalgebra_lu_solve_median.csv](docs/assets/bench/vs_nalgebra_lu_solve_median.csv)
+Raw data: [docs/assets/bench/vs_linalg_lu_solve_median.csv](docs/assets/bench/vs_linalg_lu_solve_median.csv)
 
-Summary (median time; lower is better). “la-stack vs nalgebra” is the % time reduction relative to nalgebra (positive = la-stack faster):
+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):
 
 <!-- BENCH_TABLE:lu_solve:median:new:BEGIN -->
-| D | la-stack median (ns) | nalgebra median (ns) | la-stack vs nalgebra |
-|---:|--------------------:|--------------------:|---------------------:|
-| 2 | 2.125 | 19.172 | +88.9% |
-| 3 | 13.562 | 24.082 | +43.7% |
-| 4 | 28.365 | 55.434 | +48.8% |
-| 5 | 48.567 | 76.793 | +36.8% |
-| 8 | 141.935 | 182.628 | +22.3% |
-| 16 | 642.935 | 605.115 | -6.3% |
-| 32 | 2,761.816 | 2,505.691 | -10.2% |
-| 64 | 17,009.208 | 14,696.410 | -15.7% |
+| D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | la-stack vs nalgebra | la-stack vs faer |
+|---:|--------------------:|--------------------:|----------------:|---------------------:|----------------:|
+| 2 | 2.065 | 18.375 | 160.418 | +88.8% | +98.7% |
+| 3 | 13.457 | 23.377 | 198.440 | +42.4% | +93.2% |
+| 4 | 27.750 | 54.267 | 228.744 | +48.9% | +87.9% |
+| 5 | 46.317 | 73.840 | 291.623 | +37.3% | +84.1% |
+| 8 | 138.183 | 177.982 | 389.006 | +22.4% | +64.5% |
+| 16 | 629.427 | 591.505 | 893.672 | -6.4% | +29.6% |
+| 32 | 2,688.216 | 2,503.157 | 2,908.436 | -7.4% | +7.6% |
+| 64 | 16,771.962 | 14,860.016 | 12,485.424 | -12.9% | -34.3% |
 <!-- BENCH_TABLE:lu_solve:median:new:END -->
 
 ## 📄 License

README.md

@@ -35,7 +35,7 @@ When making changes in this repo, prioritize (in order):
   - `src/lu.rs`: `Lu<const D: usize>` factorization with partial pivoting (`solve_vec`, `det`)
 - A minimal `justfile` exists for common workflows (see `just --list`).
 - The public API re-exports these items from `src/lib.rs`.
-- Dev-only benchmarks live in `benches/vs_nalgebra.rs` (Criterion + nalgebra comparison).
+- Dev-only benchmarks live in `benches/vs_linalg.rs` (Criterion + nalgebra/faer comparison).
 
 ## Publishing note
 

WARP.md

@@ -1,16 +1,58 @@
-//! Benchmark comparison between la-stack and nalgebra.
+//! Benchmark comparison between la-stack and other Rust linear algebra crates.
 //!
 //! Goal: like-for-like comparisons of the operations la-stack supports across several
 //! fixed dimensions.
 //!
 //! Notes:
-//! - Determinant is benchmarked via LU on both sides (nalgebra uses closed-forms for 1×1/2×2/3×3).
-//! - Matrix infinity norm is the maximum absolute row sum on both sides.
+//! - Determinant is benchmarked via LU on all sides (nalgebra uses closed-forms for 1×1/2×2/3×3).
+//! - Matrix infinity norm is the maximum absolute row sum on all sides.
 
 use criterion::Criterion;
+use faer::linalg::solvers::Solve;
+use faer::perm::PermRef;
 use pastey::paste;
 use std::hint::black_box;
 
+fn faer_perm_sign(p: PermRef<'_, usize>) -> f64 {
+    // Sign(det(P)) for a permutation matrix P is +1 for even permutations, -1 for odd.
+    // Parity can be computed from the number of cycles:
+    //   sign = (-1)^(n - cycles)
+    let (forward, _inverse) = p.arrays();
+    let n = forward.len();
+
+    let mut seen = vec![false; n];
+    let mut cycles = 0usize;
+
+    for start in 0..n {
+        if seen[start] {
+            continue;
+        }
+        cycles += 1;
+
+        let mut i = start;
+        while !seen[i] {
+            seen[i] = true;
+            i = forward[i];
+        }
+    }
+
+    if (n - cycles).is_multiple_of(2) {
+        1.0
+    } else {
+        -1.0
+    }
+}
+
+fn faer_det_from_partial_piv_lu(lu: &faer::linalg::solvers::PartialPivLu<f64>) -> f64 {
+    // For PA = LU with unit-lower L, det(A) = det(P) * det(U).
+    let u = lu.U();
+    let mut det = 1.0;
+    for i in 0..u.nrows() {
+        det *= u[(i, i)];
+    }
+    det * faer_perm_sign(lu.P())
+}
+
 #[inline]
 #[allow(clippy::cast_precision_loss)] // D, r, c are small integers, precision loss is not an issue.
 fn matrix_entry<const D: usize>(r: usize, c: usize) -> f64 {
@@ -81,7 +123,7 @@ fn nalgebra_inf_norm<const D: usize>(m: &nalgebra::SMatrix<f64, D, D>) -> f64 {
     max_row_sum
 }
 
-macro_rules! gen_vs_nalgebra_benches_for_dim {
+macro_rules! gen_vs_linalg_benches_for_dim {
     ($c:expr, $d:literal) => {
         paste! {{
             // Isolate each dimension's inputs to keep types and captures clean.
@@ -96,11 +138,17 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                 let nv1 = nalgebra::SVector::<f64, $d>::from_fn(|i, _| vector_entry(i, 0.0));
                 let nv2 = nalgebra::SVector::<f64, $d>::from_fn(|i, _| vector_entry(i, 1.0));
 
+                let fa = faer::Mat::<f64>::from_fn($d, $d, |r, c| matrix_entry::<$d>(r, c));
+                let frhs = faer::Mat::<f64>::from_fn($d, 1, |i, _| vector_entry(i, 0.0));
+                let fv1 = faer::Mat::<f64>::from_fn($d, 1, |i, _| vector_entry(i, 0.0));
+                let fv2 = faer::Mat::<f64>::from_fn($d, 1, |i, _| vector_entry(i, 1.0));
+
                 // Precompute LU once for solve-only / det-only benchmarks.
                 let a_lu = a
                     .lu(la_stack::DEFAULT_PIVOT_TOL)
                     .expect("matrix should be non-singular");
                 let na_lu = na.clone().lu();
+                let fa_lu = fa.partial_piv_lu();
 
                 let mut [<group_d $d>] = ($c).benchmark_group(concat!("d", stringify!($d)));
 
@@ -123,6 +171,14 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_det_via_lu", |bencher| {
+                    bencher.iter(|| {
+                        let lu = black_box(&fa).partial_piv_lu();
+                        let det = faer_det_from_partial_piv_lu(&lu);
+                        black_box(det);
+                    });
+                });
+
                 // === LU factorization ===
                 [<group_d $d>].bench_function("la_stack_lu", |bencher| {
                     bencher.iter(|| {
@@ -140,6 +196,13 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_lu", |bencher| {
+                    bencher.iter(|| {
+                        let lu = black_box(&fa).partial_piv_lu();
+                        black_box(lu);
+                    });
+                });
+
                 // === LU solve (factor + solve) ===
                 [<group_d $d>].bench_function("la_stack_lu_solve", |bencher| {
                     bencher.iter(|| {
@@ -163,6 +226,14 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_lu_solve", |bencher| {
+                    bencher.iter(|| {
+                        let lu = black_box(&fa).partial_piv_lu();
+                        let x = lu.solve(black_box(&frhs));
+                        black_box(x);
+                    });
+                });
+
                 // === Solve using a precomputed LU ===
                 [<group_d $d>].bench_function("la_stack_solve_from_lu", |bencher| {
                     bencher.iter(|| {
@@ -182,6 +253,13 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_solve_from_lu", |bencher| {
+                    bencher.iter(|| {
+                        let x = fa_lu.solve(black_box(&frhs));
+                        black_box(x);
+                    });
+                });
+
                 // === Determinant from a precomputed LU ===
                 [<group_d $d>].bench_function("la_stack_det_from_lu", |bencher| {
                     bencher.iter(|| {
@@ -197,6 +275,13 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_det_from_lu", |bencher| {
+                    bencher.iter(|| {
+                        let det = faer_det_from_partial_piv_lu(&fa_lu);
+                        black_box(det);
+                    });
+                });
+
                 // === Vector dot product ===
                 [<group_d $d>].bench_function("la_stack_dot", |bencher| {
                     bencher.iter(|| {
@@ -212,6 +297,18 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_dot", |bencher| {
+                    bencher.iter(|| {
+                        let mut sum = 0.0;
+                        let a = black_box(&fv1);
+                        let b = black_box(&fv2);
+                        for i in 0..$d {
+                            sum += a[(i, 0)] * b[(i, 0)];
+                        }
+                        black_box(sum);
+                    });
+                });
+
                 // === Vector norm squared ===
                 [<group_d $d>].bench_function("la_stack_norm2_sq", |bencher| {
                     bencher.iter(|| {
@@ -227,6 +324,18 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_norm2_sq", |bencher| {
+                    bencher.iter(|| {
+                        let mut sum = 0.0;
+                        let v = black_box(&fv1);
+                        for i in 0..$d {
+                            let x = v[(i, 0)];
+                            sum += x * x;
+                        }
+                        black_box(sum);
+                    });
+                });
+
                 // === Matrix infinity norm (max absolute row sum) ===
                 [<group_d $d>].bench_function("la_stack_inf_norm", |bencher| {
                     bencher.iter(|| {
@@ -242,6 +351,25 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
                     });
                 });
 
+                [<group_d $d>].bench_function("faer_inf_norm", |bencher| {
+                    bencher.iter(|| {
+                        let m = black_box(&fa);
+                        let mut max_row_sum = 0.0;
+
+                        for r in 0..$d {
+                            let mut row_sum = 0.0;
+                            for c in 0..$d {
+                                row_sum += m[(r, c)].abs();
+                            }
+                            if row_sum > max_row_sum {
+                                max_row_sum = row_sum;
+                            }
+                        }
+
+                        black_box(max_row_sum);
+                    });
+                });
+
                 [<group_d $d>].finish();
             }
         }}
@@ -251,15 +379,15 @@ macro_rules! gen_vs_nalgebra_benches_for_dim {
 fn main() {
     let mut c = Criterion::default().configure_from_args();
 
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 2);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 3);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 4);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 5);
+    gen_vs_linalg_benches_for_dim!(&mut c, 2);
+    gen_vs_linalg_benches_for_dim!(&mut c, 3);
+    gen_vs_linalg_benches_for_dim!(&mut c, 4);
+    gen_vs_linalg_benches_for_dim!(&mut c, 5);
 
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 8);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 16);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 32);
-    gen_vs_nalgebra_benches_for_dim!(&mut c, 64);
+    gen_vs_linalg_benches_for_dim!(&mut c, 8);
+    gen_vs_linalg_benches_for_dim!(&mut c, 16);
+    gen_vs_linalg_benches_for_dim!(&mut c, 32);
+    gen_vs_linalg_benches_for_dim!(&mut c, 64);
 
     c.final_summary();
 }

benches/vs_nalgebra.rs benches/vs_linalg.rsbenches/vs_nalgebra.rs renamed to benches/vs_linalg.rs

@@ -5,6 +5,7 @@
   "words": [
     "acgetchell",
     "blas",
+    "capsys",
     "Clippy",
     "clippy",
     "codacy",
@@ -16,6 +17,8 @@
     "f128",
     "f32",
     "f64",
+    "faer",
+    "frhs",
     "generics",
     "Getchell",
     "gnuplot",
@@ -24,6 +27,7 @@
     "keepends",
     "laerror",
     "lapack",
+    "linalg",
     "linespoints",
     "logscale",
     "lu",
@@ -38,6 +42,7 @@
     "nonfinite",
     "noplot",
     "nrhs",
+    "nrows",
     "openblas",
     "pastey",
     "patchlevel",

cspell.json

@@ -0,0 +1,9 @@
+D,la_stack,la_lo,la_hi,nalgebra,na_lo,na_hi,faer,fa_lo,fa_hi
+2,2.064783337073873,2.057689704502174,2.068888221867889,18.37497049638361,18.331137225709913,18.44905068359342,160.41828167015206,159.771806584937,161.23095742945685
+3,13.456732290936305,13.447517145467188,13.49236442840063,23.37662094350666,23.266710970814643,23.477406862396062,198.4402652520574,197.38808795625323,199.0284125157543
+4,27.75044852748178,27.599588184050052,27.785226320734317,54.266619121715074,54.15885008672231,54.36197439817226,228.74410991553628,227.56318439122288,230.055154703485
+5,46.31708492242545,46.21936017779927,46.502000712392494,73.83977618264294,73.57429654893484,74.07010889514254,291.6230709464669,290.6976886943525,292.8926041338858
+8,138.18343752426807,137.71388876876108,138.5654231809486,177.9820086523135,177.7122781945933,178.24448790485158,389.0063341187339,387.97176861931916,389.80891746278604
+16,629.4267672715603,626.9209242618742,638.6291503313327,591.5050977290701,590.798294714918,592.3790737673589,893.6719342969343,891.1105343477684,898.4784384384384
+32,2688.2164628623186,2684.227603602204,2691.7876890359166,2503.1570307509674,2500.668872475772,2506.5409929308526,2908.435890151515,2905.222288277013,2916.704963235294
+64,16771.9616225278,16754.37058346066,16817.485164835165,14860.015805946792,14773.547535211268,14956.477776495036,12485.423971036585,12471.988315217392,12501.272443181819