Add QR iteration and cleanup

3 files changed, 149 insertions, 33 deletions

diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp
index 9209a30..218b35c 100644
--- a/experiments/hilbert_qr.cpp
+++ b/experiments/hilbert_qr.cpp
@@ -1,9 +1,3 @@
-// Experiment: QR methods on Hilbert matrices
-//
-// The Hilbert matrix H[i][j] = 1/(i+j+1) is the canonical ill-conditioned
-// dense matrix.  Its condition number grows roughly as (3.5 * e)^n / sqrt(n),
-// reaching ~10^13 at n=10 and ~10^18 at n=14.
-//
 // We compare classical GS, modified GS, and Householder QR on:
 //   - reconstruction error   ||A - QR||_F
 //   - orthogonality error    ||Q^T Q - I||_F
@@ -68,9 +62,6 @@ double orthogonality_error(const QRResult& qr) {
     return std::sqrt(err);
 }
 
-// ---------------------------------------------------------------------------
-// Timing
-// ---------------------------------------------------------------------------
 
 using Clock = std::chrono::high_resolution_clock;
 using Seconds = std::chrono::duration<double>;
@@ -88,10 +79,6 @@ double min_time(Fn fn, int trials = 5) {
     return best;
 }
 
-// ---------------------------------------------------------------------------
-// Run one method on one size, return {recon, ortho, time} or nullopt on failure
-// ---------------------------------------------------------------------------
-
 using QRFn = std::function<QRResult(const Matrix&)>;
 
 struct Result { double recon, ortho, time_s; };
@@ -110,10 +97,6 @@ std::optional<Result> measure(const Matrix& A, QRFn fn) {
     }
 }
 
-// ---------------------------------------------------------------------------
-// Pretty printing
-// ---------------------------------------------------------------------------
-
 void print_row(const std::string& method, std::optional<Result> r) {
     std::cout << std::left << std::setw(16) << method;
     if (!r) {
diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp
new file mode 100644
index 0000000..8ef0ac5
--- /dev/null
+++ b/experiments/matmul.cpp
@@ -0,0 +1,149 @@
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+
+#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) && \
+      defined(__ARM_NEON) && defined(__aarch64__) && \
+      defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
+#   define MATMUL_BACKEND "NEON"
+#else
+#   define MATMUL_BACKEND "scalar"
+#endif
+
+using linalg::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(
+            "naive_matmul: lhs.cols() != rhs.rows()");
+    }
+    const std::size_t m = lhs.rows();
+    const std::size_t n = rhs.cols();
+    const std::size_t k = lhs.cols();
+
+    Matrix result(m, n, 0.0);
+
+    const double* A = lhs.data();
+    const double* B = rhs.data();
+    double*       C = result.data();
+
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            for (std::size_t p = 0; p < k; ++p)
+                C[i * n + j] += A[i * k + p] * B[p * n + j];
+
+    return result;
+}
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+// Fill an n×n matrix with a deterministic pattern so the compiler cannot
+// optimise multiplications away.
+Matrix make_matrix(std::size_t n) {
+    Matrix M(n, n);
+    const double inv = 1.0 / static_cast<double>(n + 1);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            M(i, j) = static_cast<double>(i + j + 1) * inv;
+    return M;
+}
+
+template <typename Fn>
+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();
+        const double elapsed = Seconds(t1 - t0).count();
+        if (elapsed < best) best = elapsed;
+    }
+    return best;
+}
+
+double flops(std::size_t n) {
+    const double nd = static_cast<double>(n);
+    return 2.0 * nd * nd * nd;
+}
+
+
+int main() {
+    const std::vector<std::size_t> sizes = {
+        8, 16, 32, 64, 128, 256, 512
+    };
+
+    constexpr std::size_t small_threshold = 128;
+    constexpr int         trials_small    = 9;
+    constexpr int         trials_large    = 3;
+
+    std::cout << std::string(72, '*') << "\n";
+    std::cout << "  Matmul benchmark: naive (ijk) vs SIMD (" MATMUL_BACKEND
+                 ") + transpose\n";
+    std::cout << "  C = A * B,  A and B both n×n\n";
+    std::cout << std::string(72, '*') << "\n\n";
+
+    // Column header.
+    std::cout << std::left
+              << std::setw(6)  << "n"
+              << std::setw(14) << "naive ms"
+              << std::setw(14) << "naive GFLOP/s"
+              << std::setw(14) << "SIMD ms"
+              << std::setw(14) << "SIMD GFLOP/s"
+              << std::setw(10) << "speedup"
+              << "\n";
+    std::cout << std::string(72, '-') << "\n";
+
+    for (const std::size_t n : sizes) {
+        const Matrix A = make_matrix(n);
+        const Matrix B = make_matrix(n);
+
+        const int trials = (n <= small_threshold) ? trials_small : trials_large;
+
+        volatile double sink_naive = naive_matmul(A, B)(0, 0);
+        volatile double sink_simd  = (A * B)(0, 0);
+        (void)sink_naive;
+        (void)sink_simd;
+
+        const double t_naive = min_time_s([&] { (void)naive_matmul(A, B); }, trials);
+        const double t_simd  = min_time_s([&] { (void)(A * B); },           trials);
+
+        const double fp       = flops(n);
+        const double gf_naive = fp / t_naive / 1e9;
+        const double gf_simd  = fp / t_simd  / 1e9;
+        const double speedup  = t_naive / t_simd;
+
+        std::cout << std::left  << std::setw(6) << n
+                  << std::fixed << std::setprecision(3)
+                  << std::setw(14) << t_naive * 1e3
+                  << std::setprecision(2)
+                  << std::setw(14) << gf_naive
+                  << std::setprecision(3)
+                  << std::setw(14) << t_simd * 1e3
+                  << std::setprecision(2)
+                  << std::setw(14) << gf_simd
+                  << std::setprecision(2) << std::setw(10) << speedup
+                  << "x\n";
+    }
+
+    std::cout << "\n(each cell = minimum over "
+              << trials_small << " trials for n<=" << small_threshold
+              << ", " << trials_large << " trials for larger n)\n";
+
+    return 0;
+}
diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp
index 0e71dc1..57281d8 100644
--- a/experiments/pivoting_vs_no_pivoting.cpp
+++ b/experiments/pivoting_vs_no_pivoting.cpp
@@ -1,8 +1,3 @@
-// Experiment: partial pivoting vs no-pivot LU
-//
-// Demonstrates why partial pivoting is essential for numerical stability.
-// Run the binary and inspect the residuals printed to stdout.
-
 #include "lu.hpp"
 #include "matrix.hpp"
 #include "norms.hpp"
@@ -23,7 +18,6 @@ using linalg::Vector;
 
 // ---------------------------------------------------------------------------
 // Local no-pivot LU for comparison only.
-// This is intentionally naive — it is here to show what breaks without pivoting.
 // ---------------------------------------------------------------------------
 
 struct NoPivotLU {
@@ -55,8 +49,6 @@ NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) {
     return NoPivotLU{std::move(L), std::move(U), false, 0};
 }
 
-// Solve using a no-pivot LU (L unit lower triangular, U upper triangular).
-// If the factorization failed or U is numerically singular, returns nullopt.
 std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) {
     if (f.failed) return std::nullopt;
     try {
@@ -178,8 +170,6 @@ void exp_random(std::size_t n = 8) {
 }
 
 // 2. Badly row-scaled matrix
-//    Rows differ in magnitude by ~10^14.  Without pivoting, tiny early pivots
-//    amplify round-off; with pivoting, the large-row is selected first.
 void exp_badly_scaled() {
     const Matrix A{
         {1e-14, 1.0,   2.0  },
@@ -191,21 +181,15 @@ void exp_badly_scaled() {
 }
 
 // 3. Classic pathological example for no-pivot LU.
-//    With epsilon = 1e-15, no-pivot computes a huge multiplier (1/epsilon),
-//    which causes catastrophic cancellation in the updated rows.
-//    With pivoting, we swap first and the multiplier is bounded by 1.
 void exp_epsilon_pathology() {
     constexpr double eps = 1e-15;
     const Matrix A{{eps, 1.0}, {1.0, 2.0}};
-    //  True solution of [eps 1; 1 2] * x = [1+eps; 3] is x = [1; 1].
     const Vector b{1.0 + eps, 3.0};
     run_case("Epsilon pathology [[1e-15,1],[1,2]] (classic)", A, b);
     std::cout << "  Note: exact solution is x = [1, 1]\n";
 }
 
 // 4. Matrix where no-pivot LU diverges visibly on a 4x4 example.
-//    The first pivot is small (0.001) but rows below have entries ~1000.
-//    No pivot causes multipliers of magnitude 10^6, annihilating subdiagonal info.
 void exp_amplified_multiplier() {
     const Matrix A{
         {0.001, 1.0,   0.0,   0.0  },