aboutsummaryrefslogtreecommitdiff
path: root/experiments
diff options
context:
space:
mode:
Diffstat (limited to 'experiments')
-rw-r--r--experiments/hilbert_qr.cpp17
-rw-r--r--experiments/matmul.cpp149
-rw-r--r--experiments/pivoting_vs_no_pivoting.cpp16
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 },
generated by cgit v1.2.3 (git 2.51.0) at 2026-09-01 23:48:52 +0000