Module refactor

4 files changed, 24 insertions, 1202 deletions

diff --git a/experiments/blas_comparison.cpp b/experiments/blas_comparison.cpp
deleted file mode 100644
index f6972e9..0000000
--- a/experiments/blas_comparison.cpp
+++ /dev/null
@@ -1,1127 +0,0 @@
-// blas_comparison.cpp
-//
-// Compares every operation in this linalg library against the corresponding
-// BLAS / LAPACK reference routine for correctness and performance.
-//
-// Sections:
-//   §1  Level 1 BLAS  : ddot, dnrm2
-//   §2  Level 2 BLAS  : dgemv  (matrix–vector multiply, square and rectangular)
-//   §3  Level 3 BLAS  : dgemm  (matrix–matrix multiply)
-//   §4  Triangular    : dtrsv  (forward and backward substitution)
-//   §5  LU solve      : dgesv  (single RHS and multiple RHS)
-//   §6  QR            : dgeqrf + dorgqr (vs all three of this library's QR methods)
-//   §7  Ill-conditioned: LU solve and QR on Hilbert matrices
-//
-// Accuracy metric  : compare output to BLAS/LAPACK (or known exact solution).
-// Performance metric: minimum wall-clock time over several trials.
-//
-// "this/blas" ratio < 1 means this library's implementation is faster.
-//
-// Note on LAPACK timing: calls to lapack_lu_solve() and lapack_qr() include
-// to/from column-major conversion overhead because this library's storage is row-major.
-// This is the real cost of calling LAPACK from a row-major library.
-
-#ifdef __APPLE__
-#  include <Accelerate/Accelerate.h>
-   using lapack_int_t = __CLPK_integer;
-#else
-#  include <cblas.h>
-   extern "C" {
-       void dgesv_(int*, int*, double*, int*, int*, double*, int*, int*);
-       void dgeqrf_(int*, int*, double*, int*, double*, double*, int*, int*);
-       void dorgqr_(int*, int*, int*, double*, int*, double*, double*, int*, int*);
-   }
-   using lapack_int_t = int;
-#endif
-
-#include "lu.hpp"
-#include "matrix.hpp"
-#include "norms.hpp"
-#include "qr.hpp"
-#include "triangular_solve.hpp"
-#include "vector.hpp"
-
-#include <chrono>
-#include <cmath>
-#include <cstddef>
-#include <functional>
-#include <iomanip>
-#include <iostream>
-#include <optional>
-#include <random>
-#include <stdexcept>
-#include <string>
-#include <vector>
-
-using linalg::Matrix;
-using linalg::Vector;
-using Clock   = std::chrono::high_resolution_clock;
-using Seconds = std::chrono::duration<double>;
-
-// Volatile sink prevents the compiler from eliminating timed computations.
-static volatile double g_sink = 0.0;
-
-// ===========================================================================
-// Random data (fixed seed for reproducibility)
-// ===========================================================================
-
-static std::mt19937_64 rng(0xDEADBEEF42ULL);
-
-static double rand_dbl(double lo = -1.0, double hi = 1.0) {
-    return std::uniform_real_distribution<double>(lo, hi)(rng);
-}
-
-static Vector random_vec(std::size_t n) {
-    Vector v(n);
-    for (std::size_t i = 0; i < n; ++i) v[i] = rand_dbl();
-    return v;
-}
-
-static Matrix random_mat(std::size_t rows, std::size_t cols) {
-    Matrix M(rows, cols);
-    for (std::size_t i = 0; i < rows; ++i)
-        for (std::size_t j = 0; j < cols; ++j)
-            M(i, j) = rand_dbl();
-    return M;
-}
-
-// Lower-triangular with diagonal entries in [1, 2] (well-conditioned).
-static Matrix random_lower(std::size_t n) {
-    Matrix L(n, n, 0.0);
-    for (std::size_t i = 0; i < n; ++i) {
-        for (std::size_t j = 0; j < i; ++j) L(i, j) = rand_dbl();
-        L(i, i) = 1.0 + rand_dbl(0.0, 1.0);
-    }
-    return L;
-}
-
-// Upper-triangular with diagonal entries in [1, 2].
-static Matrix random_upper(std::size_t n) {
-    Matrix U(n, n, 0.0);
-    for (std::size_t i = 0; i < n; ++i) {
-        U(i, i) = 1.0 + rand_dbl(0.0, 1.0);
-        for (std::size_t j = i + 1; j < n; ++j) U(i, j) = rand_dbl();
-    }
-    return U;
-}
-
-// Hilbert matrix H[i][j] = 1/(i+j+1).
-static Matrix hilbert(std::size_t n) {
-    Matrix H(n, n);
-    for (std::size_t i = 0; i < n; ++i)
-        for (std::size_t j = 0; j < n; ++j)
-            H(i, j) = 1.0 / static_cast<double>(i + j + 1);
-    return H;
-}
-
-// ===========================================================================
-// Error metrics
-// ===========================================================================
-
-static double vec_l2(const Vector& v) {
-    double s = 0.0;
-    for (std::size_t i = 0; i < v.size(); ++i) s += v[i] * v[i];
-    return std::sqrt(s);
-}
-
-static double vec_diff_l2(const Vector& a, const Vector& b) {
-    double s = 0.0;
-    for (std::size_t i = 0; i < a.size(); ++i) {
-        const double d = a[i] - b[i];
-        s += d * d;
-    }
-    return std::sqrt(s);
-}
-
-static double frob_diff(const Matrix& A, const Matrix& B) {
-    double s = 0.0;
-    for (std::size_t i = 0; i < A.rows(); ++i)
-        for (std::size_t j = 0; j < A.cols(); ++j) {
-            const double d = A(i, j) - B(i, j);
-            s += d * d;
-        }
-    return std::sqrt(s);
-}
-
-static double qr_recon_err(const Matrix& A, const linalg::QRResult& qr) {
-    return frob_diff(A, qr.Q * qr.R);
-}
-
-static double qr_ortho_err(const linalg::QRResult& qr) {
-    const Matrix& Q = qr.Q;
-    const std::size_t n = Q.cols();
-    const Matrix QtQ = linalg::transpose(Q) * Q;
-    double s = 0.0;
-    for (std::size_t i = 0; i < n; ++i)
-        for (std::size_t j = 0; j < n; ++j) {
-            const double d = QtQ(i, j) - (i == j ? 1.0 : 0.0);
-            s += d * d;
-        }
-    return std::sqrt(s);
-}
-
-// ===========================================================================
-// Column-major conversion (this library's Matrix is row-major; LAPACK expects col-major)
-// ===========================================================================
-
-static std::vector<double> to_col_major(const Matrix& A) {
-    const std::size_t m = A.rows(), n = A.cols();
-    std::vector<double> buf(m * n);
-    for (std::size_t i = 0; i < m; ++i)
-        for (std::size_t j = 0; j < n; ++j)
-            buf[j * m + i] = A(i, j);
-    return buf;
-}
-
-static Matrix from_col_major(const std::vector<double>& buf,
-                              std::size_t m, std::size_t n) {
-    Matrix A(m, n);
-    for (std::size_t i = 0; i < m; ++i)
-        for (std::size_t j = 0; j < n; ++j)
-            A(i, j) = buf[j * m + i];
-    return A;
-}
-
-// ===========================================================================
-// Timing
-// ===========================================================================
-
-template <typename Fn>
-static double min_time_s(Fn fn, int trials) {
-    double best = 1e30;
-    for (int t = 0; t < trials; ++t) {
-        const auto t0 = Clock::now();
-        fn();
-        const auto t1 = Clock::now();
-        best = std::min(best, Seconds(t1 - t0).count());
-    }
-    return best;
-}
-
-// ===========================================================================
-// LAPACK wrappers
-// ===========================================================================
-
-// Solve A*x = b using LAPACK dgesv_ (LU with partial pivoting).
-// Includes to/from column-major conversion.
-static Vector lapack_lu_solve(const Matrix& A, const Vector& b) {
-    const std::size_t n = A.rows();
-    lapack_int_t ni   = static_cast<lapack_int_t>(n);
-    lapack_int_t nrhs = 1;
-    lapack_int_t lda  = ni;
-    lapack_int_t ldb  = ni;
-    lapack_int_t info = 0;
-
-    std::vector<double>       a_cm(to_col_major(A));
-    std::vector<double>       b_cm(b.data(), b.data() + n);
-    std::vector<lapack_int_t> ipiv(n);
-
-    dgesv_(&ni, &nrhs, a_cm.data(), &lda, ipiv.data(),
-           b_cm.data(), &ldb, &info);
-
-    if (info != 0) throw std::runtime_error("dgesv_ failed (info=" +
-                                             std::to_string(info) + ")");
-    Vector x(n);
-    for (std::size_t i = 0; i < n; ++i) x[i] = b_cm[i];
-    return x;
-}
-
-// Solve A*X = B using LAPACK dgesv_ with multiple RHS columns.
-// Returns solution matrix X (n x nrhs), stored row-major.
-static Matrix lapack_lu_solve_multi(const Matrix& A, const Matrix& B) {
-    const std::size_t n    = A.rows();
-    const std::size_t nrhs = B.cols();
-    lapack_int_t ni    = static_cast<lapack_int_t>(n);
-    lapack_int_t nrhsi = static_cast<lapack_int_t>(nrhs);
-    lapack_int_t lda   = ni;
-    lapack_int_t ldb   = ni;
-    lapack_int_t info  = 0;
-
-    std::vector<double>       a_cm(to_col_major(A));
-    // B stored col-major for LAPACK: each RHS is a column
-    std::vector<double>       b_cm(to_col_major(B));
-    std::vector<lapack_int_t> ipiv(n);
-
-    dgesv_(&ni, &nrhsi, a_cm.data(), &lda, ipiv.data(),
-           b_cm.data(), &ldb, &info);
-
-    if (info != 0) throw std::runtime_error("dgesv_ (multi) failed");
-    return from_col_major(b_cm, n, nrhs);
-}
-
-// QR factorization via LAPACK dgeqrf_ + dorgqr_.
-// Returns thin QR (m×n Q, n×n R), including col-major conversion overhead.
-static std::optional<linalg::QRResult> lapack_qr(const Matrix& A) {
-    const std::size_t m = A.rows(), n = A.cols();
-    lapack_int_t mi  = static_cast<lapack_int_t>(m);
-    lapack_int_t ni  = static_cast<lapack_int_t>(n);
-    lapack_int_t ki  = ni;   // number of reflectors = n for square/tall A
-    lapack_int_t lda = mi;   // column-major leading dimension
-    lapack_int_t info = 0;
-
-    std::vector<double> a_cm(to_col_major(A));
-    std::vector<double> tau(n);
-
-    // --- dgeqrf: workspace query then factorize ---
-    {
-        lapack_int_t lwork = -1;
-        double wq = 0.0;
-        dgeqrf_(&mi, &ni, a_cm.data(), &lda, tau.data(), &wq, &lwork, &info);
-        if (info != 0) return std::nullopt;
-        lwork = static_cast<lapack_int_t>(wq);
-        std::vector<double> work(static_cast<std::size_t>(lwork));
-        dgeqrf_(&mi, &ni, a_cm.data(), &lda, tau.data(),
-                work.data(), &lwork, &info);
-        if (info != 0) return std::nullopt;
-    }
-
-    // Extract R from the upper triangle of a_cm *before* dorgqr overwrites it.
-    Matrix R(n, n, 0.0);
-    for (std::size_t j = 0; j < n; ++j)
-        for (std::size_t i = 0; i <= j; ++i)
-            R(i, j) = a_cm[j * m + i];
-
-    // --- dorgqr: workspace query then form explicit Q ---
-    {
-        lapack_int_t lwork = -1;
-        double wq = 0.0;
-        dorgqr_(&mi, &ni, &ki, a_cm.data(), &lda, tau.data(),
-                &wq, &lwork, &info);
-        if (info != 0) return std::nullopt;
-        lwork = static_cast<lapack_int_t>(wq);
-        std::vector<double> work(static_cast<std::size_t>(lwork));
-        dorgqr_(&mi, &ni, &ki, a_cm.data(), &lda, tau.data(),
-                work.data(), &lwork, &info);
-        if (info != 0) return std::nullopt;
-    }
-
-    Matrix Q = from_col_major(a_cm, m, n);
-    return linalg::QRResult{Q, R};
-}
-
-// ===========================================================================
-// Formatting helpers
-// ===========================================================================
-
-static void separator(char c = '=', int w = 78) {
-    std::cout << std::string(static_cast<std::size_t>(w), c) << "\n";
-}
-
-static void ratio_col(double r) {
-    // Print ratio with a directional note.
-    std::cout << std::fixed << std::setprecision(2)
-              << std::setw(10) << r
-              << (r < 1.0 ? " (faster)\n" : " (slower)\n");
-}
-
-// ===========================================================================
-// §1  Level 1 BLAS — ddot and dnrm2
-// ===========================================================================
-
-static void section_level1() {
-    std::cout << "\n"; separator();
-    std::cout << "  §1  Level 1 BLAS — dot product (ddot) and L2 norm (dnrm2)\n";
-    separator();
-    std::cout << "\n";
-
-    const std::vector<std::size_t> sizes = {64, 256, 1024, 4096, 16384, 65536};
-    constexpr int trials = 30;
-
-    // ---- ddot ----
-    std::cout << "  cblas_ddot  vs  linalg::dot\n\n";
-    std::cout << std::left
-              << std::setw(10) << "n"
-              << std::setw(18) << "|this - blas|"
-              << std::setw(14) << "blas µs"
-              << std::setw(14) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 66);
-
-    for (std::size_t n : sizes) {
-        const Vector x  = random_vec(n);
-        const Vector y  = random_vec(n);
-        const int    ni = static_cast<int>(n);
-
-        const double dot_blas = cblas_ddot(ni, x.data(), 1, y.data(), 1);
-        const double dot_ours = linalg::dot(x, y);
-
-        const double t_blas = min_time_s([&]{
-            g_sink += cblas_ddot(ni, x.data(), 1, y.data(), 1);
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            g_sink += linalg::dot(x, y);
-        }, trials);
-
-        std::cout << std::left  << std::setw(10) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(18) << std::abs(dot_ours - dot_blas)
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // ---- dnrm2 ----
-    std::cout << "\n  cblas_dnrm2  vs  linalg::norm2\n\n";
-    std::cout << std::left
-              << std::setw(10) << "n"
-              << std::setw(18) << "|this - blas|"
-              << std::setw(14) << "blas µs"
-              << std::setw(14) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 66);
-
-    for (std::size_t n : sizes) {
-        const Vector x  = random_vec(n);
-        const int    ni = static_cast<int>(n);
-
-        const double nrm_blas = cblas_dnrm2(ni, x.data(), 1);
-        const double nrm_ours = linalg::norm2(x);
-
-        const double t_blas = min_time_s([&]{
-            g_sink += cblas_dnrm2(ni, x.data(), 1);
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            g_sink += linalg::norm2(x);
-        }, trials);
-
-        std::cout << std::left  << std::setw(10) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(18) << std::abs(nrm_ours - nrm_blas)
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-}
-
-// ===========================================================================
-// §2  Level 2 BLAS — dgemv  (y = A x)
-// ===========================================================================
-
-static void section_dgemv() {
-    std::cout << "\n"; separator();
-    std::cout << "  §2  Level 2 BLAS — matrix–vector multiply (dgemv)\n";
-    separator();
-    std::cout << "\n";
-
-    constexpr int trials = 20;
-
-    auto run_dgemv = [&](const std::vector<std::size_t>& row_sizes,
-                         const std::vector<std::size_t>& col_sizes,
-                         const std::string& label) {
-        std::cout << "  " << label << "\n\n";
-        std::cout << std::left
-                  << std::setw(8)  << "rows"
-                  << std::setw(8)  << "cols"
-                  << std::setw(20) << "||y_this - y_blas||"
-                  << std::setw(14) << "blas µs"
-                  << std::setw(14) << "this µs"
-                  << std::setw(10) << "this/blas"
-                  << "\n";
-        separator('-', 74);
-
-        for (std::size_t i = 0; i < row_sizes.size(); ++i) {
-            const std::size_t m  = row_sizes[i];
-            const std::size_t k  = col_sizes[i];
-            const Matrix A       = random_mat(m, k);
-            const Vector x       = random_vec(k);
-            const int    mi      = static_cast<int>(m);
-            const int    ki      = static_cast<int>(k);
-
-            Vector y_blas(m, 0.0);
-            cblas_dgemv(CblasRowMajor, CblasNoTrans, mi, ki,
-                        1.0, A.data(), ki, x.data(), 1,
-                        0.0, y_blas.data(), 1);
-            const Vector y_ours = A * x;
-            const double err    = vec_diff_l2(y_ours, y_blas);
-
-            const double t_blas = min_time_s([&]{
-                Vector tmp(m, 0.0);
-                cblas_dgemv(CblasRowMajor, CblasNoTrans, mi, ki,
-                            1.0, A.data(), ki, x.data(), 1,
-                            0.0, tmp.data(), 1);
-                g_sink += tmp[0];
-            }, trials);
-            const double t_ours = min_time_s([&]{
-                Vector r = A * x;
-                g_sink += r[0];
-            }, trials);
-
-            std::cout << std::left  << std::setw(8) << m
-                      << std::setw(8) << k
-                      << std::scientific << std::setprecision(2)
-                      << std::setw(20) << err
-                      << std::fixed << std::setprecision(3)
-                      << std::setw(14) << t_blas * 1e6
-                      << std::setw(14) << t_ours * 1e6;
-            ratio_col(t_ours / t_blas);
-        }
-        std::cout << "\n";
-    };
-
-    // Square matrices
-    run_dgemv({8, 32, 64, 128, 256, 512, 1024},
-              {8, 32, 64, 128, 256, 512, 1024},
-              "Square  y = A*x,  A is n×n");
-
-    // Tall matrices (more rows than cols)
-    run_dgemv({256, 512, 1024, 2048},
-              { 32,  64,  128,  256},
-              "Tall  y = A*x,  A is m×k  (m >> k)");
-
-    // Wide matrices (more cols than rows)
-    run_dgemv({ 32,  64,  128,  256},
-              {256, 512, 1024, 2048},
-              "Wide  y = A*x,  A is m×k  (m << k)");
-}
-
-// ===========================================================================
-// §3  Level 3 BLAS — dgemm  (C = A B)
-// ===========================================================================
-
-static void section_dgemm() {
-    std::cout << "\n"; separator();
-    std::cout << "  §3  Level 3 BLAS — matrix–matrix multiply (dgemm)\n";
-    separator();
-    std::cout << "\n";
-
-    const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512};
-    constexpr std::size_t thresh_small   = 128;
-    constexpr int         trials_small   = 10;
-    constexpr int         trials_large   = 3;
-
-    std::cout << "  C = A*B,  all matrices n×n\n\n";
-    std::cout << std::left
-              << std::setw(8)  << "n"
-              << std::setw(20) << "||C_this - C_blas||_F"
-              << std::setw(12) << "blas ms"
-              << std::setw(12) << "this ms"
-              << std::setw(14) << "GFLOP/s blas"
-              << std::setw(14) << "GFLOP/s this"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 90);
-
-    for (std::size_t n : sizes) {
-        const Matrix A  = random_mat(n, n);
-        const Matrix B  = random_mat(n, n);
-        const int    ni = static_cast<int>(n);
-        const int    trials = (n <= thresh_small) ? trials_small : trials_large;
-
-        Matrix C_blas(n, n, 0.0);
-        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                    ni, ni, ni, 1.0, A.data(), ni, B.data(), ni,
-                    0.0, C_blas.data(), ni);
-
-        const Matrix C_ours = A * B;
-        const double err    = frob_diff(C_ours, C_blas);
-
-        const double t_blas = min_time_s([&]{
-            Matrix tmp(n, n, 0.0);
-            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                        ni, ni, ni, 1.0, A.data(), ni, B.data(), ni,
-                        0.0, tmp.data(), ni);
-            g_sink += tmp(0, 0);
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            Matrix r = A * B;
-            g_sink += r(0, 0);
-        }, trials);
-
-        const double fp_ops  = 2.0 * static_cast<double>(n)
-                                   * static_cast<double>(n)
-                                   * static_cast<double>(n);
-        const double gf_blas = fp_ops / t_blas / 1e9;
-        const double gf_ours = fp_ops / t_ours / 1e9;
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(12) << t_blas * 1e3
-                  << std::setw(12) << t_ours * 1e3
-                  << std::setprecision(2)
-                  << std::setw(14) << gf_blas
-                  << std::setw(14) << gf_ours;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // Non-square: C (m×n) = A (m×k) * B (k×n)
-    std::cout << "\n  Non-square  C = A*B,  shapes (m×k) * (k×n) -> m×n\n\n";
-    std::cout << std::left
-              << std::setw(8)  << "m"
-              << std::setw(8)  << "k"
-              << std::setw(8)  << "n"
-              << std::setw(20) << "||C_this - C_blas||_F"
-              << std::setw(14) << "blas µs"
-              << std::setw(14) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 82);
-
-    const std::vector<std::array<std::size_t,3>> shapes = {
-        {64, 32,  128},
-        {128, 64, 256},
-        {256, 128, 64},
-        {512, 32,  256},
-    };
-
-    for (const auto& [m, k, nc] : shapes) {
-        const Matrix A  = random_mat(m, k);
-        const Matrix B  = random_mat(k, nc);
-        const int    mi = static_cast<int>(m);
-        const int    ki = static_cast<int>(k);
-        const int    ni = static_cast<int>(nc);
-
-        Matrix C_blas(m, nc, 0.0);
-        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                    mi, ni, ki, 1.0, A.data(), ki, B.data(), ni,
-                    0.0, C_blas.data(), ni);
-
-        const Matrix C_ours = A * B;
-        const double err    = frob_diff(C_ours, C_blas);
-
-        const double t_blas = min_time_s([&]{
-            Matrix tmp(m, nc, 0.0);
-            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
-                        mi, ni, ki, 1.0, A.data(), ki, B.data(), ni,
-                        0.0, tmp.data(), ni);
-            g_sink += tmp(0, 0);
-        }, 10);
-        const double t_ours = min_time_s([&]{
-            Matrix r = A * B;
-            g_sink += r(0, 0);
-        }, 10);
-
-        std::cout << std::left  << std::setw(8) << m
-                  << std::setw(8) << k
-                  << std::setw(8) << nc
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // A^T * B (transposed LHS)
-    std::cout << "\n  Transposed  C = A^T * B,  A is k×m, B is k×n -> m×n\n"
-              << "  (BLAS uses CblasTrans; this library calls linalg::transpose(A) * B)\n\n";
-    std::cout << std::left
-              << std::setw(8)  << "k"
-              << std::setw(8)  << "m"
-              << std::setw(8)  << "n"
-              << std::setw(20) << "||C_this - C_blas||_F"
-              << std::setw(14) << "blas µs"
-              << std::setw(14) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 82);
-
-    for (std::size_t sz : {64UL, 128UL, 256UL}) {
-        const std::size_t k  = sz;
-        const std::size_t mm = sz / 2;
-        const std::size_t nc = sz;
-        const Matrix A  = random_mat(k, mm);   // k × m
-        const Matrix B  = random_mat(k, nc);   // k × n
-        const int    ki = static_cast<int>(k);
-        const int    mi = static_cast<int>(mm);
-        const int    ni = static_cast<int>(nc);
-
-        // BLAS: C = A^T * B  using CblasTrans for A
-        Matrix C_blas(mm, nc, 0.0);
-        cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
-                    mi, ni, ki, 1.0, A.data(), mi, B.data(), ni,
-                    0.0, C_blas.data(), ni);
-
-        const Matrix C_ours = linalg::transpose(A) * B;
-        const double err    = frob_diff(C_ours, C_blas);
-
-        const double t_blas = min_time_s([&]{
-            Matrix tmp(mm, nc, 0.0);
-            cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
-                        mi, ni, ki, 1.0, A.data(), mi, B.data(), ni,
-                        0.0, tmp.data(), ni);
-            g_sink += tmp(0, 0);
-        }, 10);
-        const double t_ours = min_time_s([&]{
-            Matrix r = linalg::transpose(A) * B;
-            g_sink += r(0, 0);
-        }, 10);
-
-        std::cout << std::left  << std::setw(8) << k
-                  << std::setw(8) << mm
-                  << std::setw(8) << nc
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-}
-
-// ===========================================================================
-// §4  Triangular solve — dtrsv vs forward/backward_substitution
-// ===========================================================================
-
-static void section_dtrsv() {
-    std::cout << "\n"; separator();
-    std::cout << "  §4  Triangular solve — dtrsv vs forward/backward_substitution\n";
-    separator();
-    std::cout << "\n";
-
-    const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512, 1024};
-    constexpr int trials = 20;
-
-    auto print_header = [] {
-        std::cout << std::left
-                  << std::setw(8)  << "n"
-                  << std::setw(20) << "||x_this - x_blas||"
-                  << std::setw(14) << "blas µs"
-                  << std::setw(14) << "this µs"
-                  << std::setw(10) << "this/blas"
-                  << "\n";
-        separator('-', 66);
-    };
-
-    // ---- Forward substitution: Lx = b ----
-    std::cout << "  Forward substitution  Lx = b  (L lower triangular, non-unit diagonal)\n\n";
-    print_header();
-
-    for (std::size_t n : sizes) {
-        const Matrix L  = random_lower(n);
-        const Vector b  = random_vec(n);
-        const int    ni = static_cast<int>(n);
-
-        Vector x_blas = b;
-        cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit,
-                    ni, L.data(), ni, x_blas.data(), 1);
-        const Vector x_ours = linalg::forward_substitution(L, b);
-        const double err    = vec_diff_l2(x_ours, x_blas);
-
-        const double t_blas = min_time_s([&]{
-            Vector tmp = b;
-            cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit,
-                        ni, L.data(), ni, tmp.data(), 1);
-            g_sink += tmp[0];
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            Vector r = linalg::forward_substitution(L, b);
-            g_sink += r[0];
-        }, trials);
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // ---- Backward substitution: Ux = b ----
-    std::cout << "\n  Backward substitution  Ux = b  (U upper triangular, non-unit diagonal)\n\n";
-    print_header();
-
-    for (std::size_t n : sizes) {
-        const Matrix U  = random_upper(n);
-        const Vector b  = random_vec(n);
-        const int    ni = static_cast<int>(n);
-
-        Vector x_blas = b;
-        cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit,
-                    ni, U.data(), ni, x_blas.data(), 1);
-        const Vector x_ours = linalg::backward_substitution(U, b);
-        const double err    = vec_diff_l2(x_ours, x_blas);
-
-        const double t_blas = min_time_s([&]{
-            Vector tmp = b;
-            cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit,
-                        ni, U.data(), ni, tmp.data(), 1);
-            g_sink += tmp[0];
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            Vector r = linalg::backward_substitution(U, b);
-            g_sink += r[0];
-        }, trials);
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // ---- Unit-diagonal forward substitution ----
-    std::cout << "\n  Forward substitution  Lx = b  (unit diagonal)\n\n";
-    print_header();
-
-    for (std::size_t n : sizes) {
-        // Build unit lower triangular
-        Matrix L = random_lower(n);
-        for (std::size_t i = 0; i < n; ++i) L(i, i) = 1.0;
-        const Vector b  = random_vec(n);
-        const int    ni = static_cast<int>(n);
-
-        Vector x_blas = b;
-        cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasUnit,
-                    ni, L.data(), ni, x_blas.data(), 1);
-        const Vector x_ours = linalg::forward_substitution(L, b,
-                                  /*singular_tolerance=*/1e-12,
-                                  /*unit_diagonal=*/true);
-        const double err = vec_diff_l2(x_ours, x_blas);
-
-        const double t_blas = min_time_s([&]{
-            Vector tmp = b;
-            cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasUnit,
-                        ni, L.data(), ni, tmp.data(), 1);
-            g_sink += tmp[0];
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            Vector r = linalg::forward_substitution(L, b, 1e-12, true);
-            g_sink += r[0];
-        }, trials);
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(20) << err
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(14) << t_blas * 1e6
-                  << std::setw(14) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-}
-
-// ===========================================================================
-// §5  LU solve — dgesv vs lu_factor + lu_solve
-// ===========================================================================
-
-static void section_lu_solve() {
-    std::cout << "\n"; separator();
-    std::cout << "  §5  LU solve — dgesv vs lu_factor + lu_solve\n";
-    separator();
-    std::cout << "\n";
-
-    // ---- Single RHS ----
-    std::cout << "  Single RHS  Ax = b\n\n";
-    std::cout << std::left
-              << std::setw(8)  << "n"
-              << std::setw(18) << "||x_this-x_blas||"
-              << std::setw(18) << "||res_this||"
-              << std::setw(18) << "||res_blas||"
-              << std::setw(12) << "blas µs"
-              << std::setw(12) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 96);
-
-    const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512};
-    constexpr std::size_t thresh = 128;
-    constexpr int ts = 20, tl = 5;
-
-    for (std::size_t n : sizes) {
-        const Matrix A  = random_mat(n, n);
-        const Vector b  = random_vec(n);
-        const int trials = (n <= thresh) ? ts : tl;
-
-        const Vector x_blas = lapack_lu_solve(A, b);
-        const auto   lu     = linalg::lu_factor(A);
-        const Vector x_ours = linalg::lu_solve(lu, b);
-
-        const Vector res_ours = (A * x_ours) - b;
-        const Vector res_blas = (A * x_blas) - b;
-
-        const double t_blas = min_time_s([&]{
-            Vector r = lapack_lu_solve(A, b);
-            g_sink += r[0];
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            auto lu2 = linalg::lu_factor(A);
-            Vector r = linalg::lu_solve(lu2, b);
-            g_sink += r[0];
-        }, trials);
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(18) << vec_diff_l2(x_ours, x_blas)
-                  << std::setw(18) << vec_l2(res_ours)
-                  << std::setw(18) << vec_l2(res_blas)
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(12) << t_blas * 1e6
-                  << std::setw(12) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-
-    // ---- Multiple RHS ----
-    // LAPACK dgesv handles multiple RHS in one shot.
-    // Our lu_solve only handles one vector at a time; we loop.
-    std::cout << "\n  Multiple RHS  AX = B  (nrhs = 8)\n"
-              << "  BLAS: one dgesv call.  Ours: lu_factor once, lu_solve 8 times.\n\n";
-    std::cout << std::left
-              << std::setw(8)  << "n"
-              << std::setw(22) << "||X_this - X_blas||_F"
-              << std::setw(18) << "||res_this||_F"
-              << std::setw(18) << "||res_blas||_F"
-              << std::setw(12) << "blas µs"
-              << std::setw(12) << "this µs"
-              << std::setw(10) << "this/blas"
-              << "\n";
-    separator('-', 100);
-
-    constexpr std::size_t nrhs = 8;
-
-    for (std::size_t n : sizes) {
-        const Matrix A  = random_mat(n, n);
-        const Matrix B  = random_mat(n, nrhs);
-        const int trials = (n <= thresh) ? ts : tl;
-
-        // LAPACK (single call, multiple RHS)
-        const Matrix X_blas = lapack_lu_solve_multi(A, B);
-
-        // Ours: factor once, solve per column
-        const auto lu = linalg::lu_factor(A);
-        Matrix X_ours(n, nrhs);
-        for (std::size_t j = 0; j < nrhs; ++j) {
-            Vector col_b(n);
-            for (std::size_t i = 0; i < n; ++i) col_b[i] = B(i, j);
-            const Vector col_x = linalg::lu_solve(lu, col_b);
-            for (std::size_t i = 0; i < n; ++i) X_ours(i, j) = col_x[i];
-        }
-
-        const double sol_err  = frob_diff(X_ours, X_blas);
-        const double res_ours = frob_diff(A * X_ours, B);
-        const double res_blas = frob_diff(A * X_blas, B);
-
-        const double t_blas = min_time_s([&]{
-            Matrix r = lapack_lu_solve_multi(A, B);
-            g_sink += r(0, 0);
-        }, trials);
-        const double t_ours = min_time_s([&]{
-            auto lu2 = linalg::lu_factor(A);
-            for (std::size_t j = 0; j < nrhs; ++j) {
-                Vector col_b(n);
-                for (std::size_t i = 0; i < n; ++i) col_b[i] = B(i, j);
-                Vector col_x = linalg::lu_solve(lu2, col_b);
-                g_sink += col_x[0];
-            }
-        }, trials);
-
-        std::cout << std::left  << std::setw(8) << n
-                  << std::scientific << std::setprecision(2)
-                  << std::setw(22) << sol_err
-                  << std::setw(18) << res_ours
-                  << std::setw(18) << res_blas
-                  << std::fixed << std::setprecision(3)
-                  << std::setw(12) << t_blas * 1e6
-                  << std::setw(12) << t_ours * 1e6;
-        ratio_col(t_ours / t_blas);
-    }
-    std::cout << "  Note: blas timing includes column-major conversion.\n";
-}
-
-// ===========================================================================
-// §6  QR factorization — dgeqrf+dorgqr vs this library's three methods
-// ===========================================================================
-
-static void section_qr() {
-    std::cout << "\n"; separator();
-    std::cout << "  §6  QR factorization — dgeqrf+dorgqr vs this library's three methods\n";
-    separator();
-    std::cout << "\n";
-    std::cout << "  Metrics per method:\n"
-              << "    ||A - QR||_F  : reconstruction error\n"
-              << "    ||Q^TQ - I||_F: orthogonality loss\n"
-              << "    time µs       : minimum wall-clock time\n"
-              << "  (lapack timing includes to/from col-major conversion)\n\n";
-
-    using OurFn = std::function<linalg::QRResult(const Matrix&)>;
-
-    const std::vector<std::pair<std::string, OurFn>> methods = {
-        {"lapack",       [](const Matrix& A) -> linalg::QRResult {
-                             auto r = lapack_qr(A);
-                             if (!r) throw std::runtime_error("lapack_qr failed");
-                             return *r;
-                         }},
-        {"classical_gs", [](const Matrix& A){ return linalg::qr_classical_gs(A); }},
-        {"modified_gs",  [](const Matrix& A){ return linalg::qr_modified_gs(A); }},
-        {"householder",  [](const Matrix& A){ return linalg::qr_householder(A); }},
-    };
-
-    auto run_qr_block = [&](const std::vector<std::size_t>& row_vec,
-                             const std::vector<std::size_t>& col_vec,
-                             const std::string& label) {
-        for (std::size_t idx = 0; idx < row_vec.size(); ++idx) {
-            const std::size_t m = row_vec[idx];
-            const std::size_t n = col_vec[idx];
-            const Matrix A = random_mat(m, n);
-            const int trials = (n <= 64) ? 15 : (n <= 128 ? 8 : 4);
-
-            separator('-', 78);
-            std::cout << "  " << label << "  m=" << m << "  n=" << n << "\n\n";
-            std::cout << std::left
-                      << std::setw(16) << "method"
-                      << std::setw(16) << "||A-QR||_F"
-                      << std::setw(16) << "||QtQ-I||_F"
-                      << std::setw(12) << "time µs"
-                      << "\n";
-            separator('-', 60);
-
-            for (const auto& [name, fn] : methods) {
-                try {
-                    const linalg::QRResult qr = fn(A);
-                    const double re = qr_recon_err(A, qr);
-                    const double oe = qr_ortho_err(qr);
-                    const double t  = min_time_s([&]{ fn(A); }, trials);
-
-                    std::cout << std::left << std::setw(16) << name
-                              << std::scientific << std::setprecision(2)
-                              << std::setw(16) << re
-                              << std::setw(16) << oe
-                              << std::fixed << std::setprecision(2)
-                              << std::setw(12) << t * 1e6
-                              << "\n";
-                } catch (const std::exception& e) {
-                    std::cout << std::left << std::setw(16) << name
-                              << "  FAILED: " << e.what() << "\n";
-                }
-            }
-            std::cout << "\n";
-        }
-    };
-
-    run_qr_block({8, 32, 64, 128, 256}, {8, 32, 64, 128, 256},
-                 "Square random");
-
-    run_qr_block({128, 256, 512, 256}, {32, 64, 64, 128},
-                 "Tall rectangular (m > n)");
-}
-
-// ===========================================================================
-// §7  Ill-conditioned accuracy — Hilbert matrices
-// ===========================================================================
-
-static void section_ill_conditioned() {
-    std::cout << "\n"; separator();
-    std::cout << "  §7  Accuracy on ill-conditioned systems (Hilbert matrices)\n";
-    separator();
-    std::cout << "\n";
-    std::cout << "  H[i][j] = 1/(i+j+1).  Condition number grows ~exponentially.\n"
-              << "  LU: true solution x* = ones  (b = H * ones).\n"
-              << "  QR: reconstruction and orthogonality errors.\n\n";
-
-    const std::vector<std::size_t> sizes = {4, 6, 8, 10, 12, 14};
-
-    // ---- LU solve ----
-    std::cout << "  LU solve on Hilbert matrices\n\n";
-    std::cout << std::left
-              << std::setw(6)  << "n"
-              << std::setw(20) << "||res_this||"
-              << std::setw(20) << "||res_blas||"
-              << std::setw(20) << "||x_this - x*||"
-              << std::setw(20) << "||x_blas - x*||"
-              << "\n";
-    separator('-', 86);
-
-    for (std::size_t n : sizes) {
-        const Matrix H    = hilbert(n);
-        const Vector ones(n, 1.0);
-        const Vector b    = H * ones;
-
-        try {
-            const auto   lu     = linalg::lu_factor(H);
-            const Vector x_ours = linalg::lu_solve(lu, b);
-            const Vector x_blas = lapack_lu_solve(H, b);
-
-            const Vector res_ours = (H * x_ours) - b;
-            const Vector res_blas = (H * x_blas) - b;
-            const Vector err_ours = x_ours - ones;
-            const Vector err_blas = x_blas - ones;
-
-            std::cout << std::left  << std::setw(6) << n
-                      << std::scientific << std::setprecision(2)
-                      << std::setw(20) << vec_l2(res_ours)
-                      << std::setw(20) << vec_l2(res_blas)
-                      << std::setw(20) << vec_l2(err_ours)
-                      << std::setw(20) << vec_l2(err_blas)
-                      << "\n";
-        } catch (const std::exception& e) {
-            std::cout << std::setw(6) << n
-                      << "  FAILED: " << e.what() << "\n";
-        }
-    }
-
-    // ---- QR on Hilbert matrices ----
-    std::cout << "\n  QR factorization on Hilbert matrices\n\n";
-    std::cout << std::left
-              << std::setw(6)  << "n"
-              << std::setw(16) << "method"
-              << std::setw(18) << "||A-QR||_F"
-              << std::setw(18) << "||QtQ-I||_F"
-              << "\n";
-    separator('-', 58);
-
-    const std::vector<std::size_t> qr_sizes = {4, 6, 8, 10, 12};
-
-    using OurFn2 = std::function<linalg::QRResult(const Matrix&)>;
-    const std::vector<std::pair<std::string, OurFn2>> methods2 = {
-        {"lapack",       [](const Matrix& A) -> linalg::QRResult {
-                             auto r = lapack_qr(A);
-                             if (!r) throw std::runtime_error("failed");
-                             return *r;
-                         }},
-        {"classical_gs", [](const Matrix& A){ return linalg::qr_classical_gs(A); }},
-        {"modified_gs",  [](const Matrix& A){ return linalg::qr_modified_gs(A); }},
-        {"householder",  [](const Matrix& A){ return linalg::qr_householder(A); }},
-    };
-
-    for (std::size_t n : qr_sizes) {
-        const Matrix H = hilbert(n);
-        bool first = true;
-        for (const auto& [name, fn] : methods2) {
-            try {
-                const linalg::QRResult qr = fn(H);
-                const double re = qr_recon_err(H, qr);
-                const double oe = qr_ortho_err(qr);
-                std::cout << std::left
-                          << std::setw(6)  << (first ? std::to_string(n) : "")
-                          << std::setw(16) << name
-                          << std::scientific << std::setprecision(2)
-                          << std::setw(18) << re
-                          << std::setw(18) << oe
-                          << "\n";
-            } catch (const std::exception& e) {
-                std::cout << std::setw(6)  << (first ? std::to_string(n) : "")
-                          << std::setw(16) << name
-                          << "  FAILED: " << e.what() << "\n";
-            }
-            first = false;
-        }
-        std::cout << "\n";
-    }
-}
-
-// ===========================================================================
-// main
-// ===========================================================================
-
-int main() {
-    separator('*');
-    std::cout << "  BLAS / LAPACK vs linalg";
-    separator('*');
-
-    section_level1();
-    section_dgemv();
-    section_dgemm();
-    section_dtrsv();
-    section_lu_solve();
-    section_qr();
-    section_ill_conditioned();
-
-    return 0;
-}
diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp
index 849bc8d..3a52194 100644
--- a/experiments/hilbert_qr.cpp
+++ b/experiments/hilbert_qr.cpp
@@ -1,20 +1,8 @@
-#include "matrix.hpp"
-#include "qr.hpp"
+import linalgebra;
+import std;
 
-#include <chrono>
-#include <cmath>
-#include <cstddef>
-#include <functional>
-#include <iomanip>
-#include <iostream>
-#include <optional>
-#include <stdexcept>
-#include <string>
-
-using linalg::Matrix;
-using linalg::QRResult;
-
-// --- Matrix construction ---
+using linalgebra::Matrix;
+using linalgebra::QRResult;
 
 Matrix hilbert(std::size_t n) {
     Matrix H(n, n);
@@ -24,8 +12,6 @@ Matrix hilbert(std::size_t n) {
     return H;
 }
 
-// --- Metrics ---
-
 double reconstruction_error(const Matrix& A, const QRResult& qr) {
     const std::size_t m = A.rows();
     const std::size_t n = A.cols();
@@ -53,7 +39,6 @@ double orthogonality_error(const QRResult& qr) {
     return std::sqrt(err);
 }
 
-
 using Clock = std::chrono::high_resolution_clock;
 using Seconds = std::chrono::duration<double>;
 
@@ -98,7 +83,6 @@ void print_row(const std::string& method, std::optional<Result> r) {
               << std::setw(10) << r->time_s * 1e6 << " µs\n";
 }
 
-
 int main() {
     std::cout << std::string(70, '*') << "\n";
     std::cout << "  Hilbert QR Experiment: comparing GS variants and Householder\n";
@@ -124,11 +108,11 @@ int main() {
         std::cout << std::string(70, ' ') << "\n";
 
         print_row("classical_gs",
-            measure(H, [](const Matrix& A) { return linalg::qr_classical_gs(A); }));
+            measure(H, [](const Matrix& A) { return linalgebra::qr_classical_gs(A); }));
         print_row("modified_gs",
-            measure(H, [](const Matrix& A) { return linalg::qr_modified_gs(A); }));
+            measure(H, [](const Matrix& A) { return linalgebra::qr_modified_gs(A); }));
         print_row("householder",
-            measure(H, [](const Matrix& A) { return linalg::qr_householder(A); }));
+            measure(H, [](const Matrix& A) { return linalgebra::qr_householder(A); }));
     }
 
     return 0;
diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp
index ed08c3c..a33b501 100644
--- a/experiments/matmul.cpp
+++ b/experiments/matmul.cpp
@@ -1,19 +1,7 @@
-#include "linalg_error.hpp"
-#include "matrix.hpp"
+import linalgebra;
+import std;
 
-#include <chrono>
-#include <cstddef>
-#include <iomanip>
-#include <iostream>
-#include <vector>
-
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__)
-#   define MATMUL_BACKEND "AVX512"
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__)
-#   define MATMUL_BACKEND "AVX2"
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__)
-#   define MATMUL_BACKEND "AVX"
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \
+#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \
       defined(__ARM_NEON) && defined(__aarch64__) && \
       defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
 #   define MATMUL_BACKEND "NEON"
@@ -21,14 +9,14 @@
 #   define MATMUL_BACKEND "scalar"
 #endif
 
-using linalg::Matrix;
+using linalgebra::Matrix;
 using Clock   = std::chrono::high_resolution_clock;
 using Seconds = std::chrono::duration<double>;
 
 
 Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) {
     if (lhs.cols() != rhs.rows()) {
-        throw linalg::DimensionMismatchError(
+        throw linalgebra::DimensionMismatchError(
             "naive_matmul: lhs.cols() != rhs.rows()");
     }
     const std::size_t m = lhs.rows();
@@ -49,8 +37,6 @@ Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) {
     return result;
 }
 
-// --- Helpers ---
-
 Matrix make_matrix(std::size_t n) {
     Matrix M(n, n);
     const double inv = 1.0 / static_cast<double>(n + 1);
diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp
index 6478d86..2f043cb 100644
--- a/experiments/pivoting_vs_no_pivoting.cpp
+++ b/experiments/pivoting_vs_no_pivoting.cpp
@@ -1,27 +1,13 @@
-#include "lu.hpp"
-#include "matrix.hpp"
-#include "norms.hpp"
-#include "triangular_solve.hpp"
-#include "vector.hpp"
+import linalgebra;
+import std;
 
-#include <cmath>
-#include <cstddef>
-#include <iomanip>
-#include <iostream>
-#include <optional>
-#include <random>
-#include <string>
-#include <vector>
-
-using linalg::Matrix;
-using linalg::Vector;
-
-// --- Local no-pivot LU for comparison only. ---
+using linalgebra::Matrix;
+using linalgebra::Vector;
 
 struct NoPivotLU {
     Matrix L;
     Matrix U;
-    bool failed = false;   
+    bool failed = false;
     std::size_t fail_step = 0;
 };
 
@@ -50,20 +36,18 @@ NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) {
 std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) {
     if (f.failed) return std::nullopt;
     try {
-        const Vector y = linalg::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true);
-        return linalg::backward_substitution(f.U, y);
+        const Vector y = linalgebra::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true);
+        return linalgebra::backward_substitution(f.U, y);
     } catch (...) {
         return std::nullopt;
     }
 }
 
-// --- Metrics ---
-
 double solve_residual(const Matrix& A, const Vector& x, const Vector& b) {
-    return linalg::norm2(A * x - b);
+    return linalgebra::norm2(A * x - b);
 }
 
-double reconstruction_error(const Matrix& A, const linalg::LUResult& lu) {
+double reconstruction_error(const Matrix& A, const linalgebra::LUResult& lu) {
     const std::size_t n = A.rows();
     Matrix PA(n, n);
     for (std::size_t i = 0; i < n; ++i)
@@ -79,8 +63,6 @@ double reconstruction_error(const Matrix& A, const linalg::LUResult& lu) {
     return std::sqrt(err);
 }
 
-// --- Reporting ---
-
 void print_header(const std::string& title) {
     std::cout << "\n" << std::string(60, '=') << "\n";
     std::cout << "  " << title << "\n";
@@ -95,8 +77,8 @@ void print_header(const std::string& title) {
 
 void report_pivoted(const Matrix& A, const Vector& b) {
     try {
-        const linalg::LUResult lu = linalg::lu_factor(A);
-        const Vector x = linalg::lu_solve(lu, b);
+        const linalgebra::LUResult lu = linalgebra::lu_factor(A);
+        const Vector x = linalgebra::lu_solve(lu, b);
         std::cout << std::left << std::setw(22) << "Pivoted LU"
                   << std::setw(20) << std::scientific << std::setprecision(3)
                   << solve_residual(A, x, b)
@@ -143,8 +125,6 @@ void run_case(const std::string& label, const Matrix& A, const Vector& b) {
     report_no_pivot(A, b);
 }
 
-// --- Experiment cases ---
-
 void exp_random(std::size_t n = 8) {
     std::mt19937 rng(42);
     std::uniform_real_distribution<double> dist(-5.0, 5.0);
@@ -165,7 +145,7 @@ void exp_badly_scaled() {
         {1.0,   3.0,   4.0  },
         {2.0,   5.0,   7.0  }
     };
-    const Vector b{1e-14 + 3.0, 8.0, 14.0};  // b = A * [1, 1, 1]
+    const Vector b{1e-14 + 3.0, 8.0, 14.0};
     run_case("Badly scaled (row norms differ by 10^14)", A, b);
 }
 
@@ -199,7 +179,6 @@ void exp_permutation() {
     run_case("Multiple row swaps required (zeros in pivot positions)", A, b);
 }
 
-
 int main() {
     std::cout << std::string(60, '*') << "\n";
     std::cout << "  Pivoting vs No-Pivoting LU Experiment\n";