diff options
| -rw-r--r-- | CMakeLists.txt | 22 | ||||
| -rw-r--r-- | experiments/blas_comparison.cpp | 1127 | ||||
| -rw-r--r-- | experiments/hilbert_qr.cpp | 22 | ||||
| -rw-r--r-- | experiments/matmul.cpp | 15 | ||||
| -rw-r--r-- | experiments/pivoting_vs_no_pivoting.cpp | 27 | ||||
| -rw-r--r-- | include/lu.hpp | 21 | ||||
| -rw-r--r-- | src/lu.cpp | 8 | ||||
| -rw-r--r-- | src/qr.cpp | 52 | ||||
| -rw-r--r-- | src/qr_iteration.cpp | 108 |
9 files changed, 1186 insertions, 216 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 1f40b0c..e345571 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -64,7 +64,27 @@ if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS) add_executable(matmul experiments/matmul.cpp) target_link_libraries(matmul PRIVATE linear_algebra::core) -endif() + +# find_package(BLAS) +# find_package(LAPACK) +# if(BLAS_FOUND AND LAPACK_FOUND) +# add_executable(blas_comparison experiments/blas_comparison.cpp) +# target_link_libraries(blas_comparison PRIVATE +# linear_algebra::core +# ${BLAS_LIBRARIES} +# ${LAPACK_LIBRARIES} +# ) +# if(DEFINED BLAS_INCLUDE_DIRS) +# target_include_directories(blas_comparison PRIVATE ${BLAS_INCLUDE_DIRS}) +# endif() +# if(APPLE) +# # Suppress deprecation warnings from Accelerate headers on recent macOS. +# target_compile_options(blas_comparison PRIVATE -Wno-deprecated-declarations) +# endif() +# else() +# message(STATUS "BLAS/LAPACK not found — skipping blas_comparison experiment") +# endif() +# endif() if(LINEAR_ALGEBRA_BUILD_TESTS) include(FetchContent) 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; +} diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp index f30d2ac..849bc8d 100644 --- a/experiments/hilbert_qr.cpp +++ b/experiments/hilbert_qr.cpp @@ -1,8 +1,3 @@ -// We compare classical GS, modified GS, and Householder QR on: -// - reconstruction error ||A - QR||_F -// - orthogonality error ||Q^T Q - I||_F -// - wall-clock time (minimum over several trials) - #include "matrix.hpp" #include "qr.hpp" @@ -19,9 +14,7 @@ using linalg::Matrix; using linalg::QRResult; -// --------------------------------------------------------------------------- -// Matrix construction -// --------------------------------------------------------------------------- +// --- Matrix construction --- Matrix hilbert(std::size_t n) { Matrix H(n, n); @@ -31,9 +24,7 @@ Matrix hilbert(std::size_t n) { return H; } -// --------------------------------------------------------------------------- -// Metrics -// --------------------------------------------------------------------------- +// --- Metrics --- double reconstruction_error(const Matrix& A, const QRResult& qr) { const std::size_t m = A.rows(); @@ -84,11 +75,9 @@ struct Result { double recon, ortho, time_s; }; std::optional<Result> measure(const Matrix& A, QRFn fn) { try { - // Run once to get metrics. const QRResult qr = fn(A); const double re = reconstruction_error(A, qr); const double oe = orthogonality_error(qr); - // Time over multiple trials. const double t = min_time([&] { fn(A); }); return Result{re, oe, t}; } catch (const std::exception&) { @@ -109,18 +98,15 @@ void print_row(const std::string& method, std::optional<Result> r) { << std::setw(10) << r->time_s * 1e6 << " µs\n"; } -// --------------------------------------------------------------------------- -// main -// --------------------------------------------------------------------------- int main() { std::cout << std::string(70, '*') << "\n"; std::cout << " Hilbert QR Experiment: comparing GS variants and Householder\n"; std::cout << std::string(70, '*') << "\n\n"; std::cout << - "H[i][j] = 1/(i+j+1). Condition number grows ~exponentially with n.\n" + "H[i][j] = 1/(i+j+1). Condition number grows ~exponentially with n.\n" "Orthogonality loss in classical GS tracks condition number directly.\n" - "Modified GS recovers ~half the lost digits. Householder is unaffected.\n\n"; + "Modified GS recovers ~half the lost digits. Householder is unaffected.\n\n"; const std::size_t sizes[] = {2, 3, 4, 5, 6, 7, 8, 10, 12}; diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp index 8ef0ac5..ed08c3c 100644 --- a/experiments/matmul.cpp +++ b/experiments/matmul.cpp @@ -49,12 +49,8 @@ Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) { return result; } -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- +// --- 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); @@ -98,7 +94,6 @@ int main() { 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" @@ -121,12 +116,12 @@ int main() { (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 t_simd = min_time_s([&] { (void)(A * B); }, trials); - const double fp = flops(n); + 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; + 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) diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp index 57281d8..6478d86 100644 --- a/experiments/pivoting_vs_no_pivoting.cpp +++ b/experiments/pivoting_vs_no_pivoting.cpp @@ -16,14 +16,12 @@ using linalg::Matrix; using linalg::Vector; -// --------------------------------------------------------------------------- -// Local no-pivot LU for comparison only. -// --------------------------------------------------------------------------- +// --- Local no-pivot LU for comparison only. --- struct NoPivotLU { Matrix L; Matrix U; - bool failed = false; // true if a zero pivot was encountered + bool failed = false; std::size_t fail_step = 0; }; @@ -59,9 +57,7 @@ std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) { } } -// --------------------------------------------------------------------------- -// Metrics -// --------------------------------------------------------------------------- +// --- Metrics --- double solve_residual(const Matrix& A, const Vector& x, const Vector& b) { return linalg::norm2(A * x - b); @@ -83,9 +79,7 @@ double reconstruction_error(const Matrix& A, const linalg::LUResult& lu) { return std::sqrt(err); } -// --------------------------------------------------------------------------- -// Reporting -// --------------------------------------------------------------------------- +// --- Reporting --- void print_header(const std::string& title) { std::cout << "\n" << std::string(60, '=') << "\n"; @@ -127,7 +121,6 @@ void report_no_pivot(const Matrix& A, const Vector& b) { << "FAILED during solve (singular U)\n"; return; } - // Compute reconstruction error without perm (no-pivot uses A directly). const Matrix LU_prod = f.L * f.U; double rec_err = 0.0; for (std::size_t i = 0; i < A.rows(); ++i) @@ -150,11 +143,8 @@ void run_case(const std::string& label, const Matrix& A, const Vector& b) { report_no_pivot(A, b); } -// --------------------------------------------------------------------------- -// Experiment cases -// --------------------------------------------------------------------------- +// --- Experiment cases --- -// 1. Random well-conditioned matrix void exp_random(std::size_t n = 8) { std::mt19937 rng(42); std::uniform_real_distribution<double> dist(-5.0, 5.0); @@ -169,7 +159,6 @@ void exp_random(std::size_t n = 8) { run_case("Random 8x8 (well-conditioned)", A, b); } -// 2. Badly row-scaled matrix void exp_badly_scaled() { const Matrix A{ {1e-14, 1.0, 2.0 }, @@ -180,7 +169,6 @@ void exp_badly_scaled() { run_case("Badly scaled (row norms differ by 10^14)", A, b); } -// 3. Classic pathological example for no-pivot LU. void exp_epsilon_pathology() { constexpr double eps = 1e-15; const Matrix A{{eps, 1.0}, {1.0, 2.0}}; @@ -189,7 +177,6 @@ void exp_epsilon_pathology() { std::cout << " Note: exact solution is x = [1, 1]\n"; } -// 4. Matrix where no-pivot LU diverges visibly on a 4x4 example. void exp_amplified_multiplier() { const Matrix A{ {0.001, 1.0, 0.0, 0.0 }, @@ -202,7 +189,6 @@ void exp_amplified_multiplier() { std::cout << " Note: exact solution is x = [1, 2, 3, 4]\n"; } -// 5. Matrix requiring multiple row swaps (permutation is non-trivial). void exp_permutation() { const Matrix A{ {0.0, 0.0, 3.0}, @@ -213,9 +199,6 @@ void exp_permutation() { run_case("Multiple row swaps required (zeros in pivot positions)", A, b); } -// --------------------------------------------------------------------------- -// main -// --------------------------------------------------------------------------- int main() { std::cout << std::string(60, '*') << "\n"; diff --git a/include/lu.hpp b/include/lu.hpp index 6717733..e7fe939 100644 --- a/include/lu.hpp +++ b/include/lu.hpp @@ -8,18 +8,6 @@ namespace linalg { -// Result of LU factorization with partial pivoting. -// -// The factorization satisfies PA = LU, where: -// P is the permutation matrix encoded by `perm` -// L is unit lower triangular (L[i][i] == 1) -// U is upper triangular -// -// `perm[i]` = index of the original row that ended up at position i. -// Applying P to a vector b means: (Pb)[i] = b[perm[i]]. -// -// `sign` is the sign of the permutation: +1 if an even number of row -// swaps were made, -1 if odd. Useful for computing det(A) = sign * prod(diag(U)). struct LUResult { Matrix L; Matrix U; @@ -27,17 +15,8 @@ struct LUResult { int sign; }; -// Compute the LU factorization of A with partial pivoting. -// -// Throws DimensionMismatchError if A is not square. -// Throws SingularMatrixError if A is (numerically) singular, i.e. any -// pivot is smaller in magnitude than `singular_tolerance`. LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12); -// Solve Ax = b given a precomputed LU factorization. -// -// Applies the stored permutation, then forward / backward substitution. -// Throws DimensionMismatchError if b.size() != lu.L.rows(). Vector lu_solve(const LUResult& lu, const Vector& b); } // namespace linalg @@ -19,10 +19,8 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { const std::size_t n = A.rows(); - // Working copy: elimination is performed in-place here. Matrix work = A; - // L starts as identity; multipliers fill the strict lower triangle. Matrix L = Matrix::zeros(n, n); for (std::size_t i = 0; i < n; ++i) { L(i, i) = 1.0; @@ -30,7 +28,6 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { Matrix U = Matrix::zeros(n, n); - // perm[i] = original row index now at position i. std::vector<std::size_t> perm(n); std::iota(perm.begin(), perm.end(), std::size_t{0}); int sign = 1; @@ -48,11 +45,9 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { } if (pivot_row != k) { - // Swap rows in the working matrix. for (std::size_t j = 0; j < n; ++j) { std::swap(work(k, j), work(pivot_row, j)); } - // Swap already-computed multipliers in L (columns 0 .. k-1). for (std::size_t j = 0; j < k; ++j) { std::swap(L(k, j), L(pivot_row, j)); } @@ -95,17 +90,14 @@ Vector lu_solve(const LUResult& lu, const Vector& b) { throw DimensionMismatchError(oss.str()); } - // Step 1: apply permutation P. (Pb)[i] = b[perm[i]] Vector pb(n); for (std::size_t i = 0; i < n; ++i) { pb[i] = b[lu.perm[i]]; } - // Step 2: forward substitution Ly = Pb (L has unit diagonal) const Vector y = forward_substitution(lu.L, pb, /*singular_tolerance=*/1e-14, /*unit_diagonal=*/true); - // Step 3: backward substitution Ux = y return backward_substitution(lu.U, y); } @@ -17,7 +17,6 @@ void require_tall(const Matrix& A, const char* name) { } } -// ||column j of M||_2 double col_norm(const Matrix& M, std::size_t j) { double s = 0.0; for (std::size_t i = 0; i < M.rows(); ++i) { @@ -26,7 +25,6 @@ double col_norm(const Matrix& M, std::size_t j) { return std::sqrt(s); } -// dot product of column j of M with column k of N (same number of rows) double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) { double s = 0.0; for (std::size_t i = 0; i < M.rows(); ++i) { @@ -37,16 +35,7 @@ double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) { } // namespace -// --------------------------------------------------------------------------- -// Classical Gram-Schmidt -// --------------------------------------------------------------------------- -// -// For column j: -// R[i][j] = <a_j, q_i> for i < j -// v = a_j - sum_i R[i][j] * q_i -// R[j][j] = ||v|| -// q_j = v / R[j][j] -// +// --- Gram-Schmidt --- QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { require_tall(A, "qr_classical_gs"); @@ -57,10 +46,8 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { Matrix R = Matrix::zeros(n, n); for (std::size_t j = 0; j < n; ++j) { - // Start with column j of A. for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j); - // Project out existing basis vectors using the *original* A column. for (std::size_t k = 0; k < j; ++k) { R(k, j) = col_dot(A, j, Q, k); // <a_j, q_k> for (std::size_t i = 0; i < m; ++i) { @@ -82,20 +69,7 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { return QRResult{std::move(Q), std::move(R)}; } -// --------------------------------------------------------------------------- -// Modified Gram-Schmidt -// --------------------------------------------------------------------------- -// -// For column j: -// v = a_j -// For k = 0 .. j-1: -// R[k][j] = <v, q_k> -// v = v - R[k][j] * q_k -// R[j][j] = ||v|| -// q_j = v / R[j][j] -// -// Each subtraction uses the already-updated v, so round-off is re-corrected -// at every sub-step rather than compounding into one subtraction. +// --- Modified Gram-Schmidt --- QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) { require_tall(A, "qr_modified_gs"); @@ -129,22 +103,7 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) { return QRResult{std::move(Q), std::move(R)}; } -// --------------------------------------------------------------------------- -// Householder QR -// --------------------------------------------------------------------------- -// -// At step k, build a Householder reflector H_k that maps R[k:, k] to -// -sign(R[k,k]) * ||R[k:,k]|| * e_1. -// -// H = I - (2 / (u^T u)) * u * u^T -// where u = x + sign(x_0) * ||x|| * e_1 (sign chosen to avoid cancellation) -// -// H is never formed explicitly. It is applied via the rank-1 update: -// M[k:, :] -= u * (2/(u^T u) * (u^T M[k:, :])) -// -// After n reflections, the working copy of A has become R (upper triangular). -// Q is accumulated by applying each H_k to an identity matrix from the left. -// The thin Q (m x n) is the first n columns of the full m x m orthogonal Q. +// --- Householder QR --- QRResult qr_householder(const Matrix& A) { require_tall(A, "qr_householder"); @@ -154,13 +113,11 @@ QRResult qr_householder(const Matrix& A) { // Will become R. Matrix work = A; - // Q accumulated as full m x m orthogonal matrix; trim to m x n. Matrix Q_full = Matrix::identity(m); for (std::size_t k = 0; k < n; ++k) { const std::size_t p = m - k; // length of the subvector - // Build Householder vector u from the subcolumn work[k:, k]. std::vector<double> u(p); for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k); @@ -172,7 +129,6 @@ QRResult qr_householder(const Matrix& A) { if (x_norm == 0.0) continue; - // sigma = sign(u[0]) * ||x|| const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm; u[0] += sigma; @@ -200,13 +156,11 @@ QRResult qr_householder(const Matrix& A) { } } - // Thin Q: first n columns of Q_full^T Matrix Q(m, n); for (std::size_t i = 0; i < m; ++i) for (std::size_t j = 0; j < n; ++j) Q(i, j) = Q_full(j, i); - // Thin R: first n rows of work Matrix R(n, n); for (std::size_t i = 0; i < n; ++i) for (std::size_t j = 0; j < n; ++j) diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp index 3e79d29..e476781 100644 --- a/src/qr_iteration.cpp +++ b/src/qr_iteration.cpp @@ -17,17 +17,11 @@ namespace linalg { namespace { -// --------------------------------------------------------------------------- -// Internal helpers -// --------------------------------------------------------------------------- // Frobenius norm of the strict lower triangle of an n×n matrix. // This is the standard convergence diagnostic for QR iteration: as A_k // approaches the real Schur form, all entries below the main diagonal // (excluding 2×2 block sub-diagonals) tend to zero. -// -// ||lower(A)||_F = sqrt( sum_{i > j} A(i,j)^2 ) -// // Ref: T&B §28; used as the convergence criterion in Algorithm 28.1. double lower_triangle_norm(const Matrix& A) { const std::size_t n = A.rows(); @@ -54,8 +48,8 @@ double lower_triangle_norm(const Matrix& A) { void extract_eigenvalues(const Matrix& T, double tol, Vector& real_out, Vector& imag_out) { const std::size_t n = T.rows(); - std::size_t out = 0; // next write position in real_out / imag_out - std::size_t i = 0; // current scan position in T + std::size_t out = 0; + std::size_t i = 0; while (i < n) { const bool is_last = (i + 1 == n); @@ -74,15 +68,15 @@ void extract_eigenvalues(const Matrix& T, double tol, // Characteristic polynomial: lambda^2 - (a+d)*lambda + (ad - bc) = 0. // Discriminant: (a-d)^2 + 4*b*c. // Ref: GVL §7.4.1. - const double a = T(i, i); - const double b = T(i, i + 1); - const double c = T(i + 1, i); - const double d = T(i + 1, i + 1); - const double tr = a + d; + const double a = T(i, i); + const double b = T(i, i + 1); + const double c = T(i + 1, i); + const double d = T(i + 1, i + 1); + const double tr = a + d; const double disc = (a - d) * (a - d) + 4.0 * b * c; if (disc >= 0.0) { - // Real eigenvalues — unusual in converged real Schur form, but + // Real eigenvalues unusual in converged real Schur form, but // handled robustly in case the block didn't fully split. const double sq = std::sqrt(disc); real_out[out] = 0.5 * (tr + sq); @@ -106,7 +100,6 @@ void extract_eigenvalues(const Matrix& T, double tol, assert(out == n); } -// Verify that A is square; throw DimensionMismatchError otherwise. void require_square(const Matrix& A, const char* fname) { if (A.rows() != A.cols()) { std::ostringstream oss; @@ -118,9 +111,7 @@ void require_square(const Matrix& A, const char* fname) { } // namespace -// --------------------------------------------------------------------------- -// Stage 1: Unshifted QR iteration -// --------------------------------------------------------------------------- +// --- Unshifted QR iteration --- // // Each step performs an orthogonal similarity transformation: // A_{k-1} = Q_k R_k (Householder QR; backward-stable) @@ -143,14 +134,9 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, require_square(A, "eigenvalues_unshifted"); const std::size_t n = A.rows(); - // Threshold for classifying a sub-diagonal entry as "zero" when reading - // eigenvalues out of the converged Schur form. Using the same value as - // the convergence tolerance is appropriate; we only reach extraction once - // ||lower(A_k)||_F < opts.tolerance. Ref: GVL §7.4.1. const double extract_tol = opts.tolerance; QRIterationResult result; - // Pre-size eigenvalue Vectors; they are always length n. result.eigenvalues_real = Vector(n, 0.0); result.eigenvalues_imag = Vector(n, 0.0); @@ -159,17 +145,14 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, static_cast<std::size_t>(opts.max_iterations)); } - // Handle the trivial 1×1 case immediately. if (n == 1) { result.eigenvalues_real[0] = A(0, 0); return result; } - // Working copy; becomes the quasi-upper-triangular Schur form A_k. Matrix Ak = A; for (int k = 0; k < opts.max_iterations; ++k) { - // --- QR step --- // Factor A_{k-1} = Q R using backward-stable Householder reflections. const QRResult qr = qr_householder(Ak); @@ -192,9 +175,6 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, } } - // Maximum iterations reached without convergence — fail loudly. - // Possible causes: eigenvalues too close in magnitude, or complex - // eigenvalue pairs that require a double shift (see Stage 2). std::ostringstream oss; oss << "eigenvalues_unshifted: did not converge in " << opts.max_iterations << " iterations " @@ -204,9 +184,7 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, throw NonConvergenceError(oss.str()); } -// --------------------------------------------------------------------------- -// Stage 2: Wilkinson-shifted QR iteration -// --------------------------------------------------------------------------- +// --- Wilkinson-shifted QR iteration --- // // The Wilkinson shift is the eigenvalue of the bottom-right 2×2 block // | a b | @@ -223,27 +201,14 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, namespace { -// Wilkinson shift: eigenvalue of the symmetric 2×2 trailing block -// | a b | -// | b d | -// that is closest to d. Only the subdiagonal entry b = A(n-1, n-2) is used -// for both off-diagonal positions; this treats the block as symmetric -// regardless of the actual superdiagonal, which is the standard convention -// (T&B Lecture 29, eq. 29.5; GVL §7.4.2). -// -// Numerically stable form avoids cancellation when |δ| >> b: -// σ = d − sign(δ) · b² / (|δ| + hypot(δ, b)) -// Discriminant δ² + b² is always ≥ 0, so no complex-shift fallback is needed. double wilkinson_shift(const Matrix& A) { const std::size_t n = A.rows(); - const double a = A(n - 2, n - 2); - const double b = A(n - 1, n - 2); // subdiagonal entry only - const double d = A(n - 1, n - 1); + const double a = A(n - 2, n - 2); + const double b = A(n - 1, n - 2); // subdiagonal entry only + const double d = A(n - 1, n - 1); const double delta = 0.5 * (a - d); - // denom = |δ| + sqrt(δ² + b²) = |δ| + hypot(δ, b) const double denom = std::abs(delta) + std::hypot(delta, b); if (denom == 0.0) return d; - // sign(δ) via (delta >= 0 ? +1 : -1); shifts toward the closer eigenvalue. const double sgn = (delta >= 0.0) ? 1.0 : -1.0; return d - sgn * (b * b) / denom; } @@ -284,24 +249,20 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) Matrix Ak = A; - // n_found: next write position (filled from index n-1 downward). std::size_t n_found = n; std::size_t active = n; // live subproblem is rows/cols 0..active-1 - // Store one eigenvalue (real) from the current trailing position. auto store_real = [&](double re) { --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = 0.0; }; - // Store a complex-conjugate pair. auto store_pair = [&](double re, double im) { --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = im; --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = -im; }; - // Extract eigenvalues from a 2×2 block and store them. auto close_2x2 = [&]() { const double a = Ak(active - 2, active - 2); const double b = Ak(active - 2, active - 1); @@ -321,7 +282,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) for (int k = 0; k < opts.max_iterations; ++k) { // --- Deflation sweep --- - // Shrink active as many times as the trailing subdiagonal allows. while (active >= 2) { const double sub = std::abs(Ak(active - 1, active - 2)); const double scale = std::abs(Ak(active - 2, active - 2)) @@ -340,7 +300,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) if (active == 2) { close_2x2(); break; } // --- Wilkinson-shifted QR step on the active × active subblock --- - // Extract submatrix (copy in). Matrix sub_mat(active, active); for (std::size_t i = 0; i < active; ++i) for (std::size_t j = 0; j < active; ++j) @@ -348,13 +307,11 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) const double sigma = wilkinson_shift(sub_mat); - // Shift, factor, unshift. for (std::size_t i = 0; i < active; ++i) sub_mat(i, i) -= sigma; const QRResult qr = qr_householder(sub_mat); sub_mat = qr.R * qr.Q; for (std::size_t i = 0; i < active; ++i) sub_mat(i, i) += sigma; - // Copy back. for (std::size_t i = 0; i < active; ++i) for (std::size_t j = 0; j < active; ++j) Ak(i, j) = sub_mat(i, j); @@ -374,9 +331,7 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) return result; } -// --------------------------------------------------------------------------- -// Stage 3a: Givens rotation -// --------------------------------------------------------------------------- +// --- Givens rotation --- GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) { const double r = std::hypot(x, y); @@ -409,10 +364,7 @@ void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const { } } -// --------------------------------------------------------------------------- -// Stage 3b: Hessenberg reduction -// --------------------------------------------------------------------------- -// +// --- Hessenberg reduction --- // For k = 0, 1, ..., n-3: // Build a Householder reflector H_k that zeros A[k+2:n, k]. // Apply from left: A[k+1:n, k:n] ← H_k * A[k+1:n, k:n] @@ -430,7 +382,6 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { Matrix Q = Matrix::identity(n); for (std::size_t k = 0; k + 2 <= n; ++k) { - // Length of the sub-vector to be zeroed: rows k+1..n-1, column k. const std::size_t p = n - k - 1; // p = n - (k+1) if (p == 0) break; @@ -438,7 +389,6 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { std::vector<double> u(p); for (std::size_t i = 0; i < p; ++i) u[i] = H(k + 1 + i, k); - // ||x|| and sigma = sign(u[0]) * ||x||. double x_norm = 0.0; for (double v : u) x_norm += v * v; x_norm = std::sqrt(x_norm); @@ -476,16 +426,13 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { for (std::size_t i = 0; i < p; ++i) Q(j, k + 1 + i) -= coeff * u[i]; } - // Zero out the numerical noise below the subdiagonal explicitly. for (std::size_t i = 1; i < p; ++i) H(k + 1 + i, k) = 0.0; } return HessenbergResult{std::move(H), std::move(Q)}; } -// --------------------------------------------------------------------------- -// Stage 3c: Hessenberg QR step via Givens rotations -// --------------------------------------------------------------------------- +// --- Hessenberg QR step via Givens rotations --- // // One shifted QR step on the upper Hessenberg matrix H: // 1. Shift: H ← H - σI. @@ -497,15 +444,13 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { // 4. Unshift: H ← H + σI. // // After the step H is again upper Hessenberg (GVL §7.4.2, Theorem 7.4.1). -// Total cost: O(n²). Ref: GVL §7.4.2; T&B Lecture 29. +// Total cost: O(n²). Ref: GVL §7.4.2. void hessenberg_qr_step(Matrix& H, double sigma) { const std::size_t n = H.rows(); - // Shift. for (std::size_t j = 0; j < n; ++j) H(j, j) -= sigma; - // Accumulate Givens rotations; apply from left as we go. std::vector<GivensRotation> gs; gs.reserve(n - 1); @@ -518,9 +463,6 @@ void hessenberg_qr_step(Matrix& H, double sigma) { gs.push_back(g); } - // Apply accumulated Givens from right (G_k^T on cols k, k+1). - // After all left applications H is upper triangular R; exploiting this, - // G_k^T only has nonzero effect on rows 0..k+1. for (std::size_t k = 0; k + 1 < n; ++k) { gs[k].apply_right(H, std::min(k + 2, n)); } @@ -529,9 +471,7 @@ void hessenberg_qr_step(Matrix& H, double sigma) { for (std::size_t j = 0; j < n; ++j) H(j, j) += sigma; } -// --------------------------------------------------------------------------- -// Stage 3d: Full practical QR algorithm -// --------------------------------------------------------------------------- +// --- Full QR algorithm --- // // Same outer deflation loop as eigenvalues_shifted, but each QR step uses // hessenberg_qr_step (O(n²) Givens rotations) instead of full Householder QR @@ -560,11 +500,9 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, return result; } - // One-time O(n³) Hessenberg reduction. HessenbergResult hr = hessenberg_reduction(A); Matrix& H = hr.H; - // Deflation bookkeeping — mirrors eigenvalues_shifted exactly. std::size_t n_found = n; std::size_t active = n; @@ -615,18 +553,15 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, if (active == 2) { close_2x2(); break; } // Wilkinson shift from trailing 2×2 of the active block. - // Inlined from wilkinson_shift() to avoid a temporary Matrix copy. - const double a_w = H(active - 2, active - 2); - const double b_w = H(active - 1, active - 2); - const double d_w = H(active - 1, active - 1); + const double a_w = H(active - 2, active - 2); + const double b_w = H(active - 1, active - 2); + const double d_w = H(active - 1, active - 1); const double delta = 0.5 * (a_w - d_w); const double denom = std::abs(delta) + std::hypot(delta, b_w); const double sigma = (denom == 0.0) ? d_w : d_w - ((delta >= 0.0) ? 1.0 : -1.0) * (b_w * b_w) / denom; // O(n²) Givens step on the active×active Hessenberg subblock. - // Copy in, step, copy out — preserves entries for already-deflated - // eigenvalues stored in the lower-right corner of H. Matrix sub_H(active, active); for (std::size_t ii = 0; ii < active; ++ii) for (std::size_t jj = 0; jj < active; ++jj) @@ -639,7 +574,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, H(ii, jj) = sub_H(ii, jj); if (opts.track_convergence) { - // Record the lower-triangle norm of the active subblock only. double s = 0.0; for (std::size_t ii = 1; ii < active; ++ii) for (std::size_t jj = 0; jj < ii; ++jj) |