aboutsummaryrefslogtreecommitdiff
path: root/experiments/blas_comparison.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'experiments/blas_comparison.cpp')
-rw-r--r--experiments/blas_comparison.cpp1127
1 files changed, 1127 insertions, 0 deletions
diff --git a/experiments/blas_comparison.cpp b/experiments/blas_comparison.cpp
new file mode 100644
index 0000000..f6972e9
--- /dev/null
+++ b/experiments/blas_comparison.cpp
@@ -0,0 +1,1127 @@
+// 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;
+}
generated by cgit v1.2.3 (git 2.51.0) at 2026-09-01 21:58:50 +0000