diff options
| -rw-r--r-- | .clangd | 2 | ||||
| -rw-r--r-- | CMakeLists.txt | 79 | ||||
| -rw-r--r-- | README.md | 51 | ||||
| -rw-r--r-- | experiments/blas_comparison.cpp | 1127 | ||||
| -rw-r--r-- | experiments/hilbert_qr.cpp | 30 | ||||
| -rw-r--r-- | experiments/matmul.cpp | 24 | ||||
| -rw-r--r-- | experiments/pivoting_vs_no_pivoting.cpp | 45 | ||||
| -rw-r--r-- | include/lu.hpp | 22 | ||||
| -rw-r--r-- | include/matrix.hpp | 48 | ||||
| -rw-r--r-- | include/norms.hpp | 9 | ||||
| -rw-r--r-- | include/qr.hpp | 39 | ||||
| -rw-r--r-- | include/qr_iteration.hpp | 165 | ||||
| -rw-r--r-- | include/triangular_solve.hpp | 22 | ||||
| -rw-r--r-- | include/vector.hpp | 47 | ||||
| -rw-r--r-- | src/linalgebra.cpp | 10 | ||||
| -rw-r--r-- | src/linalgebra_error.cpp (renamed from include/linalg_error.hpp) | 10 | ||||
| -rw-r--r-- | src/lu.cpp | 34 | ||||
| -rw-r--r-- | src/matrix.cpp | 162 | ||||
| -rw-r--r-- | src/norms.cpp | 14 | ||||
| -rw-r--r-- | src/qr.cpp | 81 | ||||
| -rw-r--r-- | src/qr_iteration.cpp | 287 | ||||
| -rw-r--r-- | src/triangular_solve.cpp | 65 | ||||
| -rw-r--r-- | src/vector.cpp | 59 | ||||
| -rw-r--r-- | tests/test_lu.cpp | 109 | ||||
| -rw-r--r-- | tests/test_matrix.cpp | 14 | ||||
| -rw-r--r-- | tests/test_qr.cpp | 131 | ||||
| -rw-r--r-- | tests/test_qr_iteration.cpp | 200 | ||||
| -rw-r--r-- | tests/test_triangular_solve.cpp | 31 | ||||
| -rw-r--r-- | tests/test_vector.cpp | 16 |
29 files changed, 557 insertions, 2376 deletions
@@ -0,0 +1,2 @@ +CompileFlags: + CompilationDatabase: build diff --git a/CMakeLists.txt b/CMakeLists.txt index e345571..37dcaf7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,45 +1,49 @@ -cmake_minimum_required(VERSION 3.20) +cmake_minimum_required(VERSION 3.30) + +set(CMAKE_EXPERIMENTAL_CXX_IMPORT_STD "d0edc3af-4c50-42ea-a356-e2862fe7a444") + +if(APPLE AND CMAKE_CXX_COMPILER MATCHES "/opt/homebrew/.*/llvm/") + set(CMAKE_CXX_COMPILER_ID_ARG1 "-B/opt/homebrew/opt/llvm/lib/c++/") + set(CMAKE_CXX_FLAGS_INIT "-B/opt/homebrew/opt/llvm/lib/c++/") +endif() project(linear_algebra_cpp VERSION 0.1.0 LANGUAGES CXX) -set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD 23) set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_CXX_EXTENSIONS OFF) +set(CMAKE_CXX_EXTENSIONS ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +set(CMAKE_CXX_MODULE_STD ON) + option(LINEAR_ALGEBRA_BUILD_TESTS "Build unit tests" ON) -set(LINEAR_ALGEBRA_SIMD "AUTO" CACHE STRING "SIMD backend for matmul: AUTO, NONE, AVX, AVX2, AVX512") -set_property(CACHE LINEAR_ALGEBRA_SIMD PROPERTY STRINGS AUTO NONE AVX AVX2 AVX512) +set(LINEAR_ALGEBRA_SIMD "AUTO" CACHE STRING "SIMD backend for matmul: AUTO (uses NEON on ARM, scalar otherwise), or NONE") +set_property(CACHE LINEAR_ALGEBRA_SIMD PROPERTY STRINGS AUTO NONE) -add_library(linear_algebra - src/vector.cpp - src/matrix.cpp - src/norms.cpp - src/triangular_solve.cpp - src/lu.cpp - src/qr.cpp - src/qr_iteration.cpp -) +add_library(linear_algebra) add_library(linear_algebra::core ALIAS linear_algebra) -target_include_directories(linear_algebra +target_sources(linear_algebra PUBLIC - ${PROJECT_SOURCE_DIR}/include + FILE_SET CXX_MODULES FILES + src/linalgebra.cpp + src/linalgebra_error.cpp + src/vector.cpp + src/matrix.cpp + src/norms.cpp + src/triangular_solve.cpp + src/lu.cpp + src/qr.cpp + src/qr_iteration.cpp ) -target_compile_features(linear_algebra PUBLIC cxx_std_20) +target_compile_features(linear_algebra PUBLIC cxx_std_23) if(CMAKE_CXX_COMPILER_ID MATCHES "Clang|GNU") target_compile_options(linear_algebra PRIVATE -Wall -Wextra -Wpedantic -Wconversion) - if(LINEAR_ALGEBRA_SIMD STREQUAL "AVX") - target_compile_options(linear_algebra PRIVATE -mavx) - elseif(LINEAR_ALGEBRA_SIMD STREQUAL "AVX2") - target_compile_options(linear_algebra PRIVATE -mavx2) - elseif(LINEAR_ALGEBRA_SIMD STREQUAL "AVX512") - target_compile_options(linear_algebra PRIVATE -mavx512f) - elseif(LINEAR_ALGEBRA_SIMD STREQUAL "NONE") + if(LINEAR_ALGEBRA_SIMD STREQUAL "NONE") target_compile_definitions(linear_algebra PRIVATE LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL=1) elseif(LINEAR_ALGEBRA_SIMD STREQUAL "AUTO") else() @@ -50,7 +54,7 @@ elseif(MSVC) if(LINEAR_ALGEBRA_SIMD STREQUAL "NONE") target_compile_definitions(linear_algebra PRIVATE LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL=1) elseif(NOT LINEAR_ALGEBRA_SIMD STREQUAL "AUTO") - message(WARNING "LINEAR_ALGEBRA_SIMD explicit x86 flags are only wired for Clang/GNU right now.") + message(FATAL_ERROR "Unsupported LINEAR_ALGEBRA_SIMD value: ${LINEAR_ALGEBRA_SIMD}") endif() endif() @@ -64,33 +68,10 @@ if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS) add_executable(matmul experiments/matmul.cpp) target_link_libraries(matmul PRIVATE linear_algebra::core) - -# 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) endif() if(LINEAR_ALGEBRA_BUILD_TESTS) + include(FetchContent) FetchContent_Declare( Catch2 @@ -1,23 +1,48 @@ -# Numerical Linear Algebra - -This repo contains a small C++ dense numerical linear algebra library for `double`, with a companion experiments directory for evaluating performance. -I mostly follow Trefethen & Bau, "Numerical Linear Algebra" and Golub & Van Loan, "Matrix Computations" -The implementation supports compile-time SIMD backends for `AVX`, `AVX2`, `AVX512`, and `NEON` on `AArch64`/`ARM64` with FP64 vector support. +This is a small C++ dense numerical linear algebra library, with a companion experiments directory +for evaluating performance. +I mostly follow Trefethen & Bau, "Numerical Linear Algebra" and Golub & Van Loan, "Matrix +Computations." +The implementation uses `NEON` SIMD on ARM64 systems when available. ## Build +The library is packaged as a C++20 named module (`linalgebra`): + +- CMake 4.1.x +- Ninja +- LLVM Clang ≥ 18 with libc++ (Homebrew LLVM 22 is what's tested; AppleClang + doesn't yet support module dependency scanning) + ```bash -cmake -S . -B build +cmake -S . -B build -G Ninja \ + -DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ cmake --build build ``` -On x86, you can explicitly choose a matmul SIMD target at configure time: +You can explicitly disable SIMD at configure time with: ```bash -cmake -S . -B build -DLINEAR_ALGEBRA_SIMD=AVX2 +cmake -S . -B build -G Ninja \ + -DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ \ + -DLINEAR_ALGEBRA_SIMD=NONE ``` -Valid values are `AUTO`, `NONE`, `AVX`, `AVX2`, and `AVX512`. `AUTO` uses the compiler's current target. `NONE` forces the scalar fallback. +Valid values are `AUTO` (uses available SIMD) and `NONE` (forces scalar fallback). + +## Usage + +Import the module: + +```cpp +import linalgebra; + +int main() { + linalgebra::Matrix A{{1.0, 2.0}, {3.0, 4.0}}; + linalgebra::Vector b{5.0, 6.0}; + auto lu = linalgebra::lu_factor(A); + auto x = linalgebra::lu_solve(lu, b); +} +``` ## Run tests @@ -30,13 +55,15 @@ ctest --test-dir build --output-on-failure - Matrix / Vector core with SIMD matmul - Triangular solvers (forward / backward substitution) - LU factorization with partial pivoting (`lu_factor`, `lu_solve`) -- QR factorization — classical GS, modified GS, and Householder (`qr_classical_gs`, `qr_modified_gs`, `qr_householder`) +- QR factorization — classical GS, modified GS, and Householder (`qr_classical_gs`, + `qr_modified_gs`, `qr_householder`) - Eigenvalue computation via QR iteration: - Unshifted QR (`eigenvalues_unshifted`) — linear convergence, T&B Algorithm 28.1 - Wilkinson-shifted QR (`eigenvalues_shifted`) — typically cubic convergence, T&B Lecture 29 - - Hessenberg + Givens QR (`eigenvalues_hessenberg`) — O(n²) per step after one O(n³) reduction; ~10–30× faster than `eigenvalues_shifted` for n ≥ 50 + - Hessenberg + Givens QR (`eigenvalues_hessenberg`) — O(n²) per step after one O(n³) reduction; + ~10–30× faster than `eigenvalues_shifted` for n ≥ 50 -TODO: +## TODO: - [ ] Cholesky factorization (cholesky) — for symmetric positive definite systems - [ ] Rank-revealing QR — Householder QR with column pivoting (qr_colpiv) - [ ] Symmetric tridiagonalization — Householder reduction before symmetric QR (tridiagonalize) diff --git a/experiments/blas_comparison.cpp b/experiments/blas_comparison.cpp deleted file mode 100644 index f6972e9..0000000 --- a/experiments/blas_comparison.cpp +++ /dev/null @@ -1,1127 +0,0 @@ -// blas_comparison.cpp -// -// Compares every operation in this linalg library against the corresponding -// BLAS / LAPACK reference routine for correctness and performance. -// -// Sections: -// §1 Level 1 BLAS : ddot, dnrm2 -// §2 Level 2 BLAS : dgemv (matrix–vector multiply, square and rectangular) -// §3 Level 3 BLAS : dgemm (matrix–matrix multiply) -// §4 Triangular : dtrsv (forward and backward substitution) -// §5 LU solve : dgesv (single RHS and multiple RHS) -// §6 QR : dgeqrf + dorgqr (vs all three of this library's QR methods) -// §7 Ill-conditioned: LU solve and QR on Hilbert matrices -// -// Accuracy metric : compare output to BLAS/LAPACK (or known exact solution). -// Performance metric: minimum wall-clock time over several trials. -// -// "this/blas" ratio < 1 means this library's implementation is faster. -// -// Note on LAPACK timing: calls to lapack_lu_solve() and lapack_qr() include -// to/from column-major conversion overhead because this library's storage is row-major. -// This is the real cost of calling LAPACK from a row-major library. - -#ifdef __APPLE__ -# include <Accelerate/Accelerate.h> - using lapack_int_t = __CLPK_integer; -#else -# include <cblas.h> - extern "C" { - void dgesv_(int*, int*, double*, int*, int*, double*, int*, int*); - void dgeqrf_(int*, int*, double*, int*, double*, double*, int*, int*); - void dorgqr_(int*, int*, int*, double*, int*, double*, double*, int*, int*); - } - using lapack_int_t = int; -#endif - -#include "lu.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "qr.hpp" -#include "triangular_solve.hpp" -#include "vector.hpp" - -#include <chrono> -#include <cmath> -#include <cstddef> -#include <functional> -#include <iomanip> -#include <iostream> -#include <optional> -#include <random> -#include <stdexcept> -#include <string> -#include <vector> - -using linalg::Matrix; -using linalg::Vector; -using Clock = std::chrono::high_resolution_clock; -using Seconds = std::chrono::duration<double>; - -// Volatile sink prevents the compiler from eliminating timed computations. -static volatile double g_sink = 0.0; - -// =========================================================================== -// Random data (fixed seed for reproducibility) -// =========================================================================== - -static std::mt19937_64 rng(0xDEADBEEF42ULL); - -static double rand_dbl(double lo = -1.0, double hi = 1.0) { - return std::uniform_real_distribution<double>(lo, hi)(rng); -} - -static Vector random_vec(std::size_t n) { - Vector v(n); - for (std::size_t i = 0; i < n; ++i) v[i] = rand_dbl(); - return v; -} - -static Matrix random_mat(std::size_t rows, std::size_t cols) { - Matrix M(rows, cols); - for (std::size_t i = 0; i < rows; ++i) - for (std::size_t j = 0; j < cols; ++j) - M(i, j) = rand_dbl(); - return M; -} - -// Lower-triangular with diagonal entries in [1, 2] (well-conditioned). -static Matrix random_lower(std::size_t n) { - Matrix L(n, n, 0.0); - for (std::size_t i = 0; i < n; ++i) { - for (std::size_t j = 0; j < i; ++j) L(i, j) = rand_dbl(); - L(i, i) = 1.0 + rand_dbl(0.0, 1.0); - } - return L; -} - -// Upper-triangular with diagonal entries in [1, 2]. -static Matrix random_upper(std::size_t n) { - Matrix U(n, n, 0.0); - for (std::size_t i = 0; i < n; ++i) { - U(i, i) = 1.0 + rand_dbl(0.0, 1.0); - for (std::size_t j = i + 1; j < n; ++j) U(i, j) = rand_dbl(); - } - return U; -} - -// Hilbert matrix H[i][j] = 1/(i+j+1). -static Matrix hilbert(std::size_t n) { - Matrix H(n, n); - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) - H(i, j) = 1.0 / static_cast<double>(i + j + 1); - return H; -} - -// =========================================================================== -// Error metrics -// =========================================================================== - -static double vec_l2(const Vector& v) { - double s = 0.0; - for (std::size_t i = 0; i < v.size(); ++i) s += v[i] * v[i]; - return std::sqrt(s); -} - -static double vec_diff_l2(const Vector& a, const Vector& b) { - double s = 0.0; - for (std::size_t i = 0; i < a.size(); ++i) { - const double d = a[i] - b[i]; - s += d * d; - } - return std::sqrt(s); -} - -static double frob_diff(const Matrix& A, const Matrix& B) { - double s = 0.0; - for (std::size_t i = 0; i < A.rows(); ++i) - for (std::size_t j = 0; j < A.cols(); ++j) { - const double d = A(i, j) - B(i, j); - s += d * d; - } - return std::sqrt(s); -} - -static double qr_recon_err(const Matrix& A, const linalg::QRResult& qr) { - return frob_diff(A, qr.Q * qr.R); -} - -static double qr_ortho_err(const linalg::QRResult& qr) { - const Matrix& Q = qr.Q; - const std::size_t n = Q.cols(); - const Matrix QtQ = linalg::transpose(Q) * Q; - double s = 0.0; - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) { - const double d = QtQ(i, j) - (i == j ? 1.0 : 0.0); - s += d * d; - } - return std::sqrt(s); -} - -// =========================================================================== -// Column-major conversion (this library's Matrix is row-major; LAPACK expects col-major) -// =========================================================================== - -static std::vector<double> to_col_major(const Matrix& A) { - const std::size_t m = A.rows(), n = A.cols(); - std::vector<double> buf(m * n); - for (std::size_t i = 0; i < m; ++i) - for (std::size_t j = 0; j < n; ++j) - buf[j * m + i] = A(i, j); - return buf; -} - -static Matrix from_col_major(const std::vector<double>& buf, - std::size_t m, std::size_t n) { - Matrix A(m, n); - for (std::size_t i = 0; i < m; ++i) - for (std::size_t j = 0; j < n; ++j) - A(i, j) = buf[j * m + i]; - return A; -} - -// =========================================================================== -// Timing -// =========================================================================== - -template <typename Fn> -static double min_time_s(Fn fn, int trials) { - double best = 1e30; - for (int t = 0; t < trials; ++t) { - const auto t0 = Clock::now(); - fn(); - const auto t1 = Clock::now(); - best = std::min(best, Seconds(t1 - t0).count()); - } - return best; -} - -// =========================================================================== -// LAPACK wrappers -// =========================================================================== - -// Solve A*x = b using LAPACK dgesv_ (LU with partial pivoting). -// Includes to/from column-major conversion. -static Vector lapack_lu_solve(const Matrix& A, const Vector& b) { - const std::size_t n = A.rows(); - lapack_int_t ni = static_cast<lapack_int_t>(n); - lapack_int_t nrhs = 1; - lapack_int_t lda = ni; - lapack_int_t ldb = ni; - lapack_int_t info = 0; - - std::vector<double> a_cm(to_col_major(A)); - std::vector<double> b_cm(b.data(), b.data() + n); - std::vector<lapack_int_t> ipiv(n); - - dgesv_(&ni, &nrhs, a_cm.data(), &lda, ipiv.data(), - b_cm.data(), &ldb, &info); - - if (info != 0) throw std::runtime_error("dgesv_ failed (info=" + - std::to_string(info) + ")"); - Vector x(n); - for (std::size_t i = 0; i < n; ++i) x[i] = b_cm[i]; - return x; -} - -// Solve A*X = B using LAPACK dgesv_ with multiple RHS columns. -// Returns solution matrix X (n x nrhs), stored row-major. -static Matrix lapack_lu_solve_multi(const Matrix& A, const Matrix& B) { - const std::size_t n = A.rows(); - const std::size_t nrhs = B.cols(); - lapack_int_t ni = static_cast<lapack_int_t>(n); - lapack_int_t nrhsi = static_cast<lapack_int_t>(nrhs); - lapack_int_t lda = ni; - lapack_int_t ldb = ni; - lapack_int_t info = 0; - - std::vector<double> a_cm(to_col_major(A)); - // B stored col-major for LAPACK: each RHS is a column - std::vector<double> b_cm(to_col_major(B)); - std::vector<lapack_int_t> ipiv(n); - - dgesv_(&ni, &nrhsi, a_cm.data(), &lda, ipiv.data(), - b_cm.data(), &ldb, &info); - - if (info != 0) throw std::runtime_error("dgesv_ (multi) failed"); - return from_col_major(b_cm, n, nrhs); -} - -// QR factorization via LAPACK dgeqrf_ + dorgqr_. -// Returns thin QR (m×n Q, n×n R), including col-major conversion overhead. -static std::optional<linalg::QRResult> lapack_qr(const Matrix& A) { - const std::size_t m = A.rows(), n = A.cols(); - lapack_int_t mi = static_cast<lapack_int_t>(m); - lapack_int_t ni = static_cast<lapack_int_t>(n); - lapack_int_t ki = ni; // number of reflectors = n for square/tall A - lapack_int_t lda = mi; // column-major leading dimension - lapack_int_t info = 0; - - std::vector<double> a_cm(to_col_major(A)); - std::vector<double> tau(n); - - // --- dgeqrf: workspace query then factorize --- - { - lapack_int_t lwork = -1; - double wq = 0.0; - dgeqrf_(&mi, &ni, a_cm.data(), &lda, tau.data(), &wq, &lwork, &info); - if (info != 0) return std::nullopt; - lwork = static_cast<lapack_int_t>(wq); - std::vector<double> work(static_cast<std::size_t>(lwork)); - dgeqrf_(&mi, &ni, a_cm.data(), &lda, tau.data(), - work.data(), &lwork, &info); - if (info != 0) return std::nullopt; - } - - // Extract R from the upper triangle of a_cm *before* dorgqr overwrites it. - Matrix R(n, n, 0.0); - for (std::size_t j = 0; j < n; ++j) - for (std::size_t i = 0; i <= j; ++i) - R(i, j) = a_cm[j * m + i]; - - // --- dorgqr: workspace query then form explicit Q --- - { - lapack_int_t lwork = -1; - double wq = 0.0; - dorgqr_(&mi, &ni, &ki, a_cm.data(), &lda, tau.data(), - &wq, &lwork, &info); - if (info != 0) return std::nullopt; - lwork = static_cast<lapack_int_t>(wq); - std::vector<double> work(static_cast<std::size_t>(lwork)); - dorgqr_(&mi, &ni, &ki, a_cm.data(), &lda, tau.data(), - work.data(), &lwork, &info); - if (info != 0) return std::nullopt; - } - - Matrix Q = from_col_major(a_cm, m, n); - return linalg::QRResult{Q, R}; -} - -// =========================================================================== -// Formatting helpers -// =========================================================================== - -static void separator(char c = '=', int w = 78) { - std::cout << std::string(static_cast<std::size_t>(w), c) << "\n"; -} - -static void ratio_col(double r) { - // Print ratio with a directional note. - std::cout << std::fixed << std::setprecision(2) - << std::setw(10) << r - << (r < 1.0 ? " (faster)\n" : " (slower)\n"); -} - -// =========================================================================== -// §1 Level 1 BLAS — ddot and dnrm2 -// =========================================================================== - -static void section_level1() { - std::cout << "\n"; separator(); - std::cout << " §1 Level 1 BLAS — dot product (ddot) and L2 norm (dnrm2)\n"; - separator(); - std::cout << "\n"; - - const std::vector<std::size_t> sizes = {64, 256, 1024, 4096, 16384, 65536}; - constexpr int trials = 30; - - // ---- ddot ---- - std::cout << " cblas_ddot vs linalg::dot\n\n"; - std::cout << std::left - << std::setw(10) << "n" - << std::setw(18) << "|this - blas|" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 66); - - for (std::size_t n : sizes) { - const Vector x = random_vec(n); - const Vector y = random_vec(n); - const int ni = static_cast<int>(n); - - const double dot_blas = cblas_ddot(ni, x.data(), 1, y.data(), 1); - const double dot_ours = linalg::dot(x, y); - - const double t_blas = min_time_s([&]{ - g_sink += cblas_ddot(ni, x.data(), 1, y.data(), 1); - }, trials); - const double t_ours = min_time_s([&]{ - g_sink += linalg::dot(x, y); - }, trials); - - std::cout << std::left << std::setw(10) << n - << std::scientific << std::setprecision(2) - << std::setw(18) << std::abs(dot_ours - dot_blas) - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - - // ---- dnrm2 ---- - std::cout << "\n cblas_dnrm2 vs linalg::norm2\n\n"; - std::cout << std::left - << std::setw(10) << "n" - << std::setw(18) << "|this - blas|" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 66); - - for (std::size_t n : sizes) { - const Vector x = random_vec(n); - const int ni = static_cast<int>(n); - - const double nrm_blas = cblas_dnrm2(ni, x.data(), 1); - const double nrm_ours = linalg::norm2(x); - - const double t_blas = min_time_s([&]{ - g_sink += cblas_dnrm2(ni, x.data(), 1); - }, trials); - const double t_ours = min_time_s([&]{ - g_sink += linalg::norm2(x); - }, trials); - - std::cout << std::left << std::setw(10) << n - << std::scientific << std::setprecision(2) - << std::setw(18) << std::abs(nrm_ours - nrm_blas) - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } -} - -// =========================================================================== -// §2 Level 2 BLAS — dgemv (y = A x) -// =========================================================================== - -static void section_dgemv() { - std::cout << "\n"; separator(); - std::cout << " §2 Level 2 BLAS — matrix–vector multiply (dgemv)\n"; - separator(); - std::cout << "\n"; - - constexpr int trials = 20; - - auto run_dgemv = [&](const std::vector<std::size_t>& row_sizes, - const std::vector<std::size_t>& col_sizes, - const std::string& label) { - std::cout << " " << label << "\n\n"; - std::cout << std::left - << std::setw(8) << "rows" - << std::setw(8) << "cols" - << std::setw(20) << "||y_this - y_blas||" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 74); - - for (std::size_t i = 0; i < row_sizes.size(); ++i) { - const std::size_t m = row_sizes[i]; - const std::size_t k = col_sizes[i]; - const Matrix A = random_mat(m, k); - const Vector x = random_vec(k); - const int mi = static_cast<int>(m); - const int ki = static_cast<int>(k); - - Vector y_blas(m, 0.0); - cblas_dgemv(CblasRowMajor, CblasNoTrans, mi, ki, - 1.0, A.data(), ki, x.data(), 1, - 0.0, y_blas.data(), 1); - const Vector y_ours = A * x; - const double err = vec_diff_l2(y_ours, y_blas); - - const double t_blas = min_time_s([&]{ - Vector tmp(m, 0.0); - cblas_dgemv(CblasRowMajor, CblasNoTrans, mi, ki, - 1.0, A.data(), ki, x.data(), 1, - 0.0, tmp.data(), 1); - g_sink += tmp[0]; - }, trials); - const double t_ours = min_time_s([&]{ - Vector r = A * x; - g_sink += r[0]; - }, trials); - - std::cout << std::left << std::setw(8) << m - << std::setw(8) << k - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - std::cout << "\n"; - }; - - // Square matrices - run_dgemv({8, 32, 64, 128, 256, 512, 1024}, - {8, 32, 64, 128, 256, 512, 1024}, - "Square y = A*x, A is n×n"); - - // Tall matrices (more rows than cols) - run_dgemv({256, 512, 1024, 2048}, - { 32, 64, 128, 256}, - "Tall y = A*x, A is m×k (m >> k)"); - - // Wide matrices (more cols than rows) - run_dgemv({ 32, 64, 128, 256}, - {256, 512, 1024, 2048}, - "Wide y = A*x, A is m×k (m << k)"); -} - -// =========================================================================== -// §3 Level 3 BLAS — dgemm (C = A B) -// =========================================================================== - -static void section_dgemm() { - std::cout << "\n"; separator(); - std::cout << " §3 Level 3 BLAS — matrix–matrix multiply (dgemm)\n"; - separator(); - std::cout << "\n"; - - const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512}; - constexpr std::size_t thresh_small = 128; - constexpr int trials_small = 10; - constexpr int trials_large = 3; - - std::cout << " C = A*B, all matrices n×n\n\n"; - std::cout << std::left - << std::setw(8) << "n" - << std::setw(20) << "||C_this - C_blas||_F" - << std::setw(12) << "blas ms" - << std::setw(12) << "this ms" - << std::setw(14) << "GFLOP/s blas" - << std::setw(14) << "GFLOP/s this" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 90); - - for (std::size_t n : sizes) { - const Matrix A = random_mat(n, n); - const Matrix B = random_mat(n, n); - const int ni = static_cast<int>(n); - const int trials = (n <= thresh_small) ? trials_small : trials_large; - - Matrix C_blas(n, n, 0.0); - cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - ni, ni, ni, 1.0, A.data(), ni, B.data(), ni, - 0.0, C_blas.data(), ni); - - const Matrix C_ours = A * B; - const double err = frob_diff(C_ours, C_blas); - - const double t_blas = min_time_s([&]{ - Matrix tmp(n, n, 0.0); - cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - ni, ni, ni, 1.0, A.data(), ni, B.data(), ni, - 0.0, tmp.data(), ni); - g_sink += tmp(0, 0); - }, trials); - const double t_ours = min_time_s([&]{ - Matrix r = A * B; - g_sink += r(0, 0); - }, trials); - - const double fp_ops = 2.0 * static_cast<double>(n) - * static_cast<double>(n) - * static_cast<double>(n); - const double gf_blas = fp_ops / t_blas / 1e9; - const double gf_ours = fp_ops / t_ours / 1e9; - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(12) << t_blas * 1e3 - << std::setw(12) << t_ours * 1e3 - << std::setprecision(2) - << std::setw(14) << gf_blas - << std::setw(14) << gf_ours; - ratio_col(t_ours / t_blas); - } - - // Non-square: C (m×n) = A (m×k) * B (k×n) - std::cout << "\n Non-square C = A*B, shapes (m×k) * (k×n) -> m×n\n\n"; - std::cout << std::left - << std::setw(8) << "m" - << std::setw(8) << "k" - << std::setw(8) << "n" - << std::setw(20) << "||C_this - C_blas||_F" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 82); - - const std::vector<std::array<std::size_t,3>> shapes = { - {64, 32, 128}, - {128, 64, 256}, - {256, 128, 64}, - {512, 32, 256}, - }; - - for (const auto& [m, k, nc] : shapes) { - const Matrix A = random_mat(m, k); - const Matrix B = random_mat(k, nc); - const int mi = static_cast<int>(m); - const int ki = static_cast<int>(k); - const int ni = static_cast<int>(nc); - - Matrix C_blas(m, nc, 0.0); - cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - mi, ni, ki, 1.0, A.data(), ki, B.data(), ni, - 0.0, C_blas.data(), ni); - - const Matrix C_ours = A * B; - const double err = frob_diff(C_ours, C_blas); - - const double t_blas = min_time_s([&]{ - Matrix tmp(m, nc, 0.0); - cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - mi, ni, ki, 1.0, A.data(), ki, B.data(), ni, - 0.0, tmp.data(), ni); - g_sink += tmp(0, 0); - }, 10); - const double t_ours = min_time_s([&]{ - Matrix r = A * B; - g_sink += r(0, 0); - }, 10); - - std::cout << std::left << std::setw(8) << m - << std::setw(8) << k - << std::setw(8) << nc - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - - // A^T * B (transposed LHS) - std::cout << "\n Transposed C = A^T * B, A is k×m, B is k×n -> m×n\n" - << " (BLAS uses CblasTrans; this library calls linalg::transpose(A) * B)\n\n"; - std::cout << std::left - << std::setw(8) << "k" - << std::setw(8) << "m" - << std::setw(8) << "n" - << std::setw(20) << "||C_this - C_blas||_F" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 82); - - for (std::size_t sz : {64UL, 128UL, 256UL}) { - const std::size_t k = sz; - const std::size_t mm = sz / 2; - const std::size_t nc = sz; - const Matrix A = random_mat(k, mm); // k × m - const Matrix B = random_mat(k, nc); // k × n - const int ki = static_cast<int>(k); - const int mi = static_cast<int>(mm); - const int ni = static_cast<int>(nc); - - // BLAS: C = A^T * B using CblasTrans for A - Matrix C_blas(mm, nc, 0.0); - cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, - mi, ni, ki, 1.0, A.data(), mi, B.data(), ni, - 0.0, C_blas.data(), ni); - - const Matrix C_ours = linalg::transpose(A) * B; - const double err = frob_diff(C_ours, C_blas); - - const double t_blas = min_time_s([&]{ - Matrix tmp(mm, nc, 0.0); - cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, - mi, ni, ki, 1.0, A.data(), mi, B.data(), ni, - 0.0, tmp.data(), ni); - g_sink += tmp(0, 0); - }, 10); - const double t_ours = min_time_s([&]{ - Matrix r = linalg::transpose(A) * B; - g_sink += r(0, 0); - }, 10); - - std::cout << std::left << std::setw(8) << k - << std::setw(8) << mm - << std::setw(8) << nc - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } -} - -// =========================================================================== -// §4 Triangular solve — dtrsv vs forward/backward_substitution -// =========================================================================== - -static void section_dtrsv() { - std::cout << "\n"; separator(); - std::cout << " §4 Triangular solve — dtrsv vs forward/backward_substitution\n"; - separator(); - std::cout << "\n"; - - const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512, 1024}; - constexpr int trials = 20; - - auto print_header = [] { - std::cout << std::left - << std::setw(8) << "n" - << std::setw(20) << "||x_this - x_blas||" - << std::setw(14) << "blas µs" - << std::setw(14) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 66); - }; - - // ---- Forward substitution: Lx = b ---- - std::cout << " Forward substitution Lx = b (L lower triangular, non-unit diagonal)\n\n"; - print_header(); - - for (std::size_t n : sizes) { - const Matrix L = random_lower(n); - const Vector b = random_vec(n); - const int ni = static_cast<int>(n); - - Vector x_blas = b; - cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, - ni, L.data(), ni, x_blas.data(), 1); - const Vector x_ours = linalg::forward_substitution(L, b); - const double err = vec_diff_l2(x_ours, x_blas); - - const double t_blas = min_time_s([&]{ - Vector tmp = b; - cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, - ni, L.data(), ni, tmp.data(), 1); - g_sink += tmp[0]; - }, trials); - const double t_ours = min_time_s([&]{ - Vector r = linalg::forward_substitution(L, b); - g_sink += r[0]; - }, trials); - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - - // ---- Backward substitution: Ux = b ---- - std::cout << "\n Backward substitution Ux = b (U upper triangular, non-unit diagonal)\n\n"; - print_header(); - - for (std::size_t n : sizes) { - const Matrix U = random_upper(n); - const Vector b = random_vec(n); - const int ni = static_cast<int>(n); - - Vector x_blas = b; - cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit, - ni, U.data(), ni, x_blas.data(), 1); - const Vector x_ours = linalg::backward_substitution(U, b); - const double err = vec_diff_l2(x_ours, x_blas); - - const double t_blas = min_time_s([&]{ - Vector tmp = b; - cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit, - ni, U.data(), ni, tmp.data(), 1); - g_sink += tmp[0]; - }, trials); - const double t_ours = min_time_s([&]{ - Vector r = linalg::backward_substitution(U, b); - g_sink += r[0]; - }, trials); - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - - // ---- Unit-diagonal forward substitution ---- - std::cout << "\n Forward substitution Lx = b (unit diagonal)\n\n"; - print_header(); - - for (std::size_t n : sizes) { - // Build unit lower triangular - Matrix L = random_lower(n); - for (std::size_t i = 0; i < n; ++i) L(i, i) = 1.0; - const Vector b = random_vec(n); - const int ni = static_cast<int>(n); - - Vector x_blas = b; - cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasUnit, - ni, L.data(), ni, x_blas.data(), 1); - const Vector x_ours = linalg::forward_substitution(L, b, - /*singular_tolerance=*/1e-12, - /*unit_diagonal=*/true); - const double err = vec_diff_l2(x_ours, x_blas); - - const double t_blas = min_time_s([&]{ - Vector tmp = b; - cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasUnit, - ni, L.data(), ni, tmp.data(), 1); - g_sink += tmp[0]; - }, trials); - const double t_ours = min_time_s([&]{ - Vector r = linalg::forward_substitution(L, b, 1e-12, true); - g_sink += r[0]; - }, trials); - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(20) << err - << std::fixed << std::setprecision(3) - << std::setw(14) << t_blas * 1e6 - << std::setw(14) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } -} - -// =========================================================================== -// §5 LU solve — dgesv vs lu_factor + lu_solve -// =========================================================================== - -static void section_lu_solve() { - std::cout << "\n"; separator(); - std::cout << " §5 LU solve — dgesv vs lu_factor + lu_solve\n"; - separator(); - std::cout << "\n"; - - // ---- Single RHS ---- - std::cout << " Single RHS Ax = b\n\n"; - std::cout << std::left - << std::setw(8) << "n" - << std::setw(18) << "||x_this-x_blas||" - << std::setw(18) << "||res_this||" - << std::setw(18) << "||res_blas||" - << std::setw(12) << "blas µs" - << std::setw(12) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 96); - - const std::vector<std::size_t> sizes = {8, 32, 64, 128, 256, 512}; - constexpr std::size_t thresh = 128; - constexpr int ts = 20, tl = 5; - - for (std::size_t n : sizes) { - const Matrix A = random_mat(n, n); - const Vector b = random_vec(n); - const int trials = (n <= thresh) ? ts : tl; - - const Vector x_blas = lapack_lu_solve(A, b); - const auto lu = linalg::lu_factor(A); - const Vector x_ours = linalg::lu_solve(lu, b); - - const Vector res_ours = (A * x_ours) - b; - const Vector res_blas = (A * x_blas) - b; - - const double t_blas = min_time_s([&]{ - Vector r = lapack_lu_solve(A, b); - g_sink += r[0]; - }, trials); - const double t_ours = min_time_s([&]{ - auto lu2 = linalg::lu_factor(A); - Vector r = linalg::lu_solve(lu2, b); - g_sink += r[0]; - }, trials); - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(18) << vec_diff_l2(x_ours, x_blas) - << std::setw(18) << vec_l2(res_ours) - << std::setw(18) << vec_l2(res_blas) - << std::fixed << std::setprecision(3) - << std::setw(12) << t_blas * 1e6 - << std::setw(12) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - - // ---- Multiple RHS ---- - // LAPACK dgesv handles multiple RHS in one shot. - // Our lu_solve only handles one vector at a time; we loop. - std::cout << "\n Multiple RHS AX = B (nrhs = 8)\n" - << " BLAS: one dgesv call. Ours: lu_factor once, lu_solve 8 times.\n\n"; - std::cout << std::left - << std::setw(8) << "n" - << std::setw(22) << "||X_this - X_blas||_F" - << std::setw(18) << "||res_this||_F" - << std::setw(18) << "||res_blas||_F" - << std::setw(12) << "blas µs" - << std::setw(12) << "this µs" - << std::setw(10) << "this/blas" - << "\n"; - separator('-', 100); - - constexpr std::size_t nrhs = 8; - - for (std::size_t n : sizes) { - const Matrix A = random_mat(n, n); - const Matrix B = random_mat(n, nrhs); - const int trials = (n <= thresh) ? ts : tl; - - // LAPACK (single call, multiple RHS) - const Matrix X_blas = lapack_lu_solve_multi(A, B); - - // Ours: factor once, solve per column - const auto lu = linalg::lu_factor(A); - Matrix X_ours(n, nrhs); - for (std::size_t j = 0; j < nrhs; ++j) { - Vector col_b(n); - for (std::size_t i = 0; i < n; ++i) col_b[i] = B(i, j); - const Vector col_x = linalg::lu_solve(lu, col_b); - for (std::size_t i = 0; i < n; ++i) X_ours(i, j) = col_x[i]; - } - - const double sol_err = frob_diff(X_ours, X_blas); - const double res_ours = frob_diff(A * X_ours, B); - const double res_blas = frob_diff(A * X_blas, B); - - const double t_blas = min_time_s([&]{ - Matrix r = lapack_lu_solve_multi(A, B); - g_sink += r(0, 0); - }, trials); - const double t_ours = min_time_s([&]{ - auto lu2 = linalg::lu_factor(A); - for (std::size_t j = 0; j < nrhs; ++j) { - Vector col_b(n); - for (std::size_t i = 0; i < n; ++i) col_b[i] = B(i, j); - Vector col_x = linalg::lu_solve(lu2, col_b); - g_sink += col_x[0]; - } - }, trials); - - std::cout << std::left << std::setw(8) << n - << std::scientific << std::setprecision(2) - << std::setw(22) << sol_err - << std::setw(18) << res_ours - << std::setw(18) << res_blas - << std::fixed << std::setprecision(3) - << std::setw(12) << t_blas * 1e6 - << std::setw(12) << t_ours * 1e6; - ratio_col(t_ours / t_blas); - } - std::cout << " Note: blas timing includes column-major conversion.\n"; -} - -// =========================================================================== -// §6 QR factorization — dgeqrf+dorgqr vs this library's three methods -// =========================================================================== - -static void section_qr() { - std::cout << "\n"; separator(); - std::cout << " §6 QR factorization — dgeqrf+dorgqr vs this library's three methods\n"; - separator(); - std::cout << "\n"; - std::cout << " Metrics per method:\n" - << " ||A - QR||_F : reconstruction error\n" - << " ||Q^TQ - I||_F: orthogonality loss\n" - << " time µs : minimum wall-clock time\n" - << " (lapack timing includes to/from col-major conversion)\n\n"; - - using OurFn = std::function<linalg::QRResult(const Matrix&)>; - - const std::vector<std::pair<std::string, OurFn>> methods = { - {"lapack", [](const Matrix& A) -> linalg::QRResult { - auto r = lapack_qr(A); - if (!r) throw std::runtime_error("lapack_qr failed"); - return *r; - }}, - {"classical_gs", [](const Matrix& A){ return linalg::qr_classical_gs(A); }}, - {"modified_gs", [](const Matrix& A){ return linalg::qr_modified_gs(A); }}, - {"householder", [](const Matrix& A){ return linalg::qr_householder(A); }}, - }; - - auto run_qr_block = [&](const std::vector<std::size_t>& row_vec, - const std::vector<std::size_t>& col_vec, - const std::string& label) { - for (std::size_t idx = 0; idx < row_vec.size(); ++idx) { - const std::size_t m = row_vec[idx]; - const std::size_t n = col_vec[idx]; - const Matrix A = random_mat(m, n); - const int trials = (n <= 64) ? 15 : (n <= 128 ? 8 : 4); - - separator('-', 78); - std::cout << " " << label << " m=" << m << " n=" << n << "\n\n"; - std::cout << std::left - << std::setw(16) << "method" - << std::setw(16) << "||A-QR||_F" - << std::setw(16) << "||QtQ-I||_F" - << std::setw(12) << "time µs" - << "\n"; - separator('-', 60); - - for (const auto& [name, fn] : methods) { - try { - const linalg::QRResult qr = fn(A); - const double re = qr_recon_err(A, qr); - const double oe = qr_ortho_err(qr); - const double t = min_time_s([&]{ fn(A); }, trials); - - std::cout << std::left << std::setw(16) << name - << std::scientific << std::setprecision(2) - << std::setw(16) << re - << std::setw(16) << oe - << std::fixed << std::setprecision(2) - << std::setw(12) << t * 1e6 - << "\n"; - } catch (const std::exception& e) { - std::cout << std::left << std::setw(16) << name - << " FAILED: " << e.what() << "\n"; - } - } - std::cout << "\n"; - } - }; - - run_qr_block({8, 32, 64, 128, 256}, {8, 32, 64, 128, 256}, - "Square random"); - - run_qr_block({128, 256, 512, 256}, {32, 64, 64, 128}, - "Tall rectangular (m > n)"); -} - -// =========================================================================== -// §7 Ill-conditioned accuracy — Hilbert matrices -// =========================================================================== - -static void section_ill_conditioned() { - std::cout << "\n"; separator(); - std::cout << " §7 Accuracy on ill-conditioned systems (Hilbert matrices)\n"; - separator(); - std::cout << "\n"; - std::cout << " H[i][j] = 1/(i+j+1). Condition number grows ~exponentially.\n" - << " LU: true solution x* = ones (b = H * ones).\n" - << " QR: reconstruction and orthogonality errors.\n\n"; - - const std::vector<std::size_t> sizes = {4, 6, 8, 10, 12, 14}; - - // ---- LU solve ---- - std::cout << " LU solve on Hilbert matrices\n\n"; - std::cout << std::left - << std::setw(6) << "n" - << std::setw(20) << "||res_this||" - << std::setw(20) << "||res_blas||" - << std::setw(20) << "||x_this - x*||" - << std::setw(20) << "||x_blas - x*||" - << "\n"; - separator('-', 86); - - for (std::size_t n : sizes) { - const Matrix H = hilbert(n); - const Vector ones(n, 1.0); - const Vector b = H * ones; - - try { - const auto lu = linalg::lu_factor(H); - const Vector x_ours = linalg::lu_solve(lu, b); - const Vector x_blas = lapack_lu_solve(H, b); - - const Vector res_ours = (H * x_ours) - b; - const Vector res_blas = (H * x_blas) - b; - const Vector err_ours = x_ours - ones; - const Vector err_blas = x_blas - ones; - - std::cout << std::left << std::setw(6) << n - << std::scientific << std::setprecision(2) - << std::setw(20) << vec_l2(res_ours) - << std::setw(20) << vec_l2(res_blas) - << std::setw(20) << vec_l2(err_ours) - << std::setw(20) << vec_l2(err_blas) - << "\n"; - } catch (const std::exception& e) { - std::cout << std::setw(6) << n - << " FAILED: " << e.what() << "\n"; - } - } - - // ---- QR on Hilbert matrices ---- - std::cout << "\n QR factorization on Hilbert matrices\n\n"; - std::cout << std::left - << std::setw(6) << "n" - << std::setw(16) << "method" - << std::setw(18) << "||A-QR||_F" - << std::setw(18) << "||QtQ-I||_F" - << "\n"; - separator('-', 58); - - const std::vector<std::size_t> qr_sizes = {4, 6, 8, 10, 12}; - - using OurFn2 = std::function<linalg::QRResult(const Matrix&)>; - const std::vector<std::pair<std::string, OurFn2>> methods2 = { - {"lapack", [](const Matrix& A) -> linalg::QRResult { - auto r = lapack_qr(A); - if (!r) throw std::runtime_error("failed"); - return *r; - }}, - {"classical_gs", [](const Matrix& A){ return linalg::qr_classical_gs(A); }}, - {"modified_gs", [](const Matrix& A){ return linalg::qr_modified_gs(A); }}, - {"householder", [](const Matrix& A){ return linalg::qr_householder(A); }}, - }; - - for (std::size_t n : qr_sizes) { - const Matrix H = hilbert(n); - bool first = true; - for (const auto& [name, fn] : methods2) { - try { - const linalg::QRResult qr = fn(H); - const double re = qr_recon_err(H, qr); - const double oe = qr_ortho_err(qr); - std::cout << std::left - << std::setw(6) << (first ? std::to_string(n) : "") - << std::setw(16) << name - << std::scientific << std::setprecision(2) - << std::setw(18) << re - << std::setw(18) << oe - << "\n"; - } catch (const std::exception& e) { - std::cout << std::setw(6) << (first ? std::to_string(n) : "") - << std::setw(16) << name - << " FAILED: " << e.what() << "\n"; - } - first = false; - } - std::cout << "\n"; - } -} - -// =========================================================================== -// main -// =========================================================================== - -int main() { - separator('*'); - std::cout << " BLAS / LAPACK vs linalg"; - separator('*'); - - section_level1(); - section_dgemv(); - section_dgemm(); - section_dtrsv(); - section_lu_solve(); - section_qr(); - section_ill_conditioned(); - - return 0; -} diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp index 849bc8d..3a52194 100644 --- a/experiments/hilbert_qr.cpp +++ b/experiments/hilbert_qr.cpp @@ -1,20 +1,8 @@ -#include "matrix.hpp" -#include "qr.hpp" +import linalgebra; +import std; -#include <chrono> -#include <cmath> -#include <cstddef> -#include <functional> -#include <iomanip> -#include <iostream> -#include <optional> -#include <stdexcept> -#include <string> - -using linalg::Matrix; -using linalg::QRResult; - -// --- Matrix construction --- +using linalgebra::Matrix; +using linalgebra::QRResult; Matrix hilbert(std::size_t n) { Matrix H(n, n); @@ -24,8 +12,6 @@ Matrix hilbert(std::size_t n) { return H; } -// --- Metrics --- - double reconstruction_error(const Matrix& A, const QRResult& qr) { const std::size_t m = A.rows(); const std::size_t n = A.cols(); @@ -53,7 +39,6 @@ double orthogonality_error(const QRResult& qr) { return std::sqrt(err); } - using Clock = std::chrono::high_resolution_clock; using Seconds = std::chrono::duration<double>; @@ -98,7 +83,6 @@ void print_row(const std::string& method, std::optional<Result> r) { << std::setw(10) << r->time_s * 1e6 << " µs\n"; } - int main() { std::cout << std::string(70, '*') << "\n"; std::cout << " Hilbert QR Experiment: comparing GS variants and Householder\n"; @@ -124,11 +108,11 @@ int main() { std::cout << std::string(70, ' ') << "\n"; print_row("classical_gs", - measure(H, [](const Matrix& A) { return linalg::qr_classical_gs(A); })); + measure(H, [](const Matrix& A) { return linalgebra::qr_classical_gs(A); })); print_row("modified_gs", - measure(H, [](const Matrix& A) { return linalg::qr_modified_gs(A); })); + measure(H, [](const Matrix& A) { return linalgebra::qr_modified_gs(A); })); print_row("householder", - measure(H, [](const Matrix& A) { return linalg::qr_householder(A); })); + measure(H, [](const Matrix& A) { return linalgebra::qr_householder(A); })); } return 0; diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp index ed08c3c..a33b501 100644 --- a/experiments/matmul.cpp +++ b/experiments/matmul.cpp @@ -1,19 +1,7 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" +import linalgebra; +import std; -#include <chrono> -#include <cstddef> -#include <iomanip> -#include <iostream> -#include <vector> - -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__) -# define MATMUL_BACKEND "AVX512" -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__) -# define MATMUL_BACKEND "AVX2" -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) -# define MATMUL_BACKEND "AVX" -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \ +#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \ defined(__ARM_NEON) && defined(__aarch64__) && \ defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC) # define MATMUL_BACKEND "NEON" @@ -21,14 +9,14 @@ # define MATMUL_BACKEND "scalar" #endif -using linalg::Matrix; +using linalgebra::Matrix; using Clock = std::chrono::high_resolution_clock; using Seconds = std::chrono::duration<double>; Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) { if (lhs.cols() != rhs.rows()) { - throw linalg::DimensionMismatchError( + throw linalgebra::DimensionMismatchError( "naive_matmul: lhs.cols() != rhs.rows()"); } const std::size_t m = lhs.rows(); @@ -49,8 +37,6 @@ Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) { return result; } -// --- Helpers --- - Matrix make_matrix(std::size_t n) { Matrix M(n, n); const double inv = 1.0 / static_cast<double>(n + 1); diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp index 6478d86..2f043cb 100644 --- a/experiments/pivoting_vs_no_pivoting.cpp +++ b/experiments/pivoting_vs_no_pivoting.cpp @@ -1,27 +1,13 @@ -#include "lu.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "triangular_solve.hpp" -#include "vector.hpp" +import linalgebra; +import std; -#include <cmath> -#include <cstddef> -#include <iomanip> -#include <iostream> -#include <optional> -#include <random> -#include <string> -#include <vector> - -using linalg::Matrix; -using linalg::Vector; - -// --- Local no-pivot LU for comparison only. --- +using linalgebra::Matrix; +using linalgebra::Vector; struct NoPivotLU { Matrix L; Matrix U; - bool failed = false; + bool failed = false; std::size_t fail_step = 0; }; @@ -50,20 +36,18 @@ NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) { std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) { if (f.failed) return std::nullopt; try { - const Vector y = linalg::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true); - return linalg::backward_substitution(f.U, y); + const Vector y = linalgebra::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true); + return linalgebra::backward_substitution(f.U, y); } catch (...) { return std::nullopt; } } -// --- Metrics --- - double solve_residual(const Matrix& A, const Vector& x, const Vector& b) { - return linalg::norm2(A * x - b); + return linalgebra::norm2(A * x - b); } -double reconstruction_error(const Matrix& A, const linalg::LUResult& lu) { +double reconstruction_error(const Matrix& A, const linalgebra::LUResult& lu) { const std::size_t n = A.rows(); Matrix PA(n, n); for (std::size_t i = 0; i < n; ++i) @@ -79,8 +63,6 @@ double reconstruction_error(const Matrix& A, const linalg::LUResult& lu) { return std::sqrt(err); } -// --- Reporting --- - void print_header(const std::string& title) { std::cout << "\n" << std::string(60, '=') << "\n"; std::cout << " " << title << "\n"; @@ -95,8 +77,8 @@ void print_header(const std::string& title) { void report_pivoted(const Matrix& A, const Vector& b) { try { - const linalg::LUResult lu = linalg::lu_factor(A); - const Vector x = linalg::lu_solve(lu, b); + const linalgebra::LUResult lu = linalgebra::lu_factor(A); + const Vector x = linalgebra::lu_solve(lu, b); std::cout << std::left << std::setw(22) << "Pivoted LU" << std::setw(20) << std::scientific << std::setprecision(3) << solve_residual(A, x, b) @@ -143,8 +125,6 @@ void run_case(const std::string& label, const Matrix& A, const Vector& b) { report_no_pivot(A, b); } -// --- Experiment cases --- - void exp_random(std::size_t n = 8) { std::mt19937 rng(42); std::uniform_real_distribution<double> dist(-5.0, 5.0); @@ -165,7 +145,7 @@ void exp_badly_scaled() { {1.0, 3.0, 4.0 }, {2.0, 5.0, 7.0 } }; - const Vector b{1e-14 + 3.0, 8.0, 14.0}; // b = A * [1, 1, 1] + const Vector b{1e-14 + 3.0, 8.0, 14.0}; run_case("Badly scaled (row norms differ by 10^14)", A, b); } @@ -199,7 +179,6 @@ void exp_permutation() { run_case("Multiple row swaps required (zeros in pivot positions)", A, b); } - int main() { std::cout << std::string(60, '*') << "\n"; std::cout << " Pivoting vs No-Pivoting LU Experiment\n"; diff --git a/include/lu.hpp b/include/lu.hpp deleted file mode 100644 index e7fe939..0000000 --- a/include/lu.hpp +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -#include <cstddef> -#include <vector> - -#include "matrix.hpp" -#include "vector.hpp" - -namespace linalg { - -struct LUResult { - Matrix L; - Matrix U; - std::vector<std::size_t> perm; - int sign; -}; - -LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12); - -Vector lu_solve(const LUResult& lu, const Vector& b); - -} // namespace linalg diff --git a/include/matrix.hpp b/include/matrix.hpp deleted file mode 100644 index d78f814..0000000 --- a/include/matrix.hpp +++ /dev/null @@ -1,48 +0,0 @@ -#pragma once - -#include <cstddef> -#include <initializer_list> -#include <vector> -#include "vector.hpp" - -namespace linalg { - -class Matrix { -public: - Matrix() = default; - Matrix(std::size_t rows, std::size_t cols); - Matrix(std::size_t rows, std::size_t cols, double value); - Matrix(std::initializer_list<std::initializer_list<double>> values); - - [[nodiscard]] std::size_t rows() const noexcept; - [[nodiscard]] std::size_t cols() const noexcept; - [[nodiscard]] bool empty() const noexcept; - - double& operator()(std::size_t i, std::size_t j); - const double& operator()(std::size_t i, std::size_t j) const; - - void fill(double value); - - double* data() noexcept; - const double* data() const noexcept; - - static Matrix identity(std::size_t n); - static Matrix zeros(std::size_t rows, std::size_t cols); - -private: - [[nodiscard]] std::size_t index(std::size_t i, std::size_t j) const; - void check_bounds(std::size_t i, std::size_t j) const; - - std::size_t rows_ = 0; - std::size_t cols_ = 0; - std::vector<double> data_; -}; - -Matrix transpose(const Matrix& matrix); -Matrix operator+(const Matrix& lhs, const Matrix& rhs); -Matrix operator-(const Matrix& lhs, const Matrix& rhs); -Vector operator*(const Matrix& matrix, const Vector& vector); -Matrix operator*(const Matrix& lhs, const Matrix& rhs); - -} // namespace linalg - diff --git a/include/norms.hpp b/include/norms.hpp deleted file mode 100644 index 85fad41..0000000 --- a/include/norms.hpp +++ /dev/null @@ -1,9 +0,0 @@ -#pragma once - -#include "vector.hpp" - -namespace linalg { - -double norm2(const Vector& vector); - -} diff --git a/include/qr.hpp b/include/qr.hpp deleted file mode 100644 index c9e1638..0000000 --- a/include/qr.hpp +++ /dev/null @@ -1,39 +0,0 @@ -#pragma once - -#include "matrix.hpp" -#include "vector.hpp" - -namespace linalg { - -struct QRResult { - Matrix Q; - Matrix R; -}; - -// Classical Gram-Schmidt. -// Mathematically natural but numerically fragile: orthogonality of Q -// degrades rapidly on ill-conditioned inputs. -// Provided for comparison — prefer modified_gs or householder in practice. -// -// Throws DimensionMismatchError if rows < cols. -// Throws SingularMatrixError if a column is (nearly) linearly dependent. -QRResult qr_classical_gs(const Matrix& A, double zero_tolerance = 1e-14); - -// Modified Gram-Schmidt. -// Subtracts each projection immediately on the running vector rather than -// on the original column. Algebraically equivalent to classical GS but -// numerically much better — round-off stays local instead of accumulating. -// -// Same exceptions as classical GS. -QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14); - -// Householder QR. -// Applies a sequence of orthogonal reflections to zero out below-diagonal -// entries column by column. Backward-stable and the standard choice for -// dense QR. Works correctly on rank-deficient matrices (zero pivots -// produce zero diagonal entries in R without throwing). -// -// Throws DimensionMismatchError if rows < cols. -QRResult qr_householder(const Matrix& A); - -} // namespace linalg diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp deleted file mode 100644 index ccf48cd..0000000 --- a/include/qr_iteration.hpp +++ /dev/null @@ -1,165 +0,0 @@ -#pragma once - -// QR iteration for eigenvalue computation -// -// Refs: -// Trefethen & Bau, "Numerical Linear Algebra" (T&B) -// Lecture 25 — Eigenvalue algorithms -// Lecture 26 — Schur factorisation -// Lecture 28 — The QR algorithm (unshifted) -// Lecture 29 — The QR algorithm with shifts -// Golub & Van Loan, "Matrix Computations" 4th ed. (GVL) -// §7.3 — The Unshifted QR Algorithm -// §7.4 — The Shifted QR Algorithm -// §7.4.2 — Wilkinson shift - -#include <vector> - -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "vector.hpp" - -namespace linalg { - -// --------------------------------------------------------------------------- -// Options -// --------------------------------------------------------------------------- - -// All defaults are consistent with the recommendations in T&B Lecture 28. -struct QRIterationOptions { - // Convergence threshold. Iteration halts once the Frobenius norm of the - // strict lower triangle of A_k falls below this value. - // Ref: T&B §28; GVL §7.3. - double tolerance = 1e-10; - - int max_iterations = 1000; - - // When true, the Frobenius norm of the strict lower triangle is recorded - // after every QR step and returned in QRIterationResult::convergence_history. - bool track_convergence = false; -}; - - -struct QRIterationResult { - // For symmetric inputs all imaginary parts are zero. - // Complex-conjugate pairs from 2×2 Schur blocks appear as +/-imag entries. - Vector eigenvalues_real; - Vector eigenvalues_imag; - - int iterations = 0; - - // Populated only when QRIterationOptions::track_convergence is true. - // Entry k is ||lower(A_k)||_F after the k-th QR step. - // std::vector is used here because linalg::Vector has no push_back; - // convergence_history is a plain time-series container, not a math object. - std::vector<double> convergence_history; -}; - -// Algorithm (T&B Algorithm 28.1): -// -// A_0 = A -// for k = 1, 2, ...: -// factor A_{k-1} = Q_k R_k (Householder QR) -// set A_k = R_k Q_k (orthogonal similarity: preserves eigenvalues) -// -// The iterates A_k converge to the real Schur form of A: a quasi-upper- -// triangular matrix with 1×1 blocks (real eigenvalue) and 2×2 blocks -// (complex-conjugate pair) on the diagonal. -// -// Convergence rate: linear. Per-step factor ~= |lambda_{j+1} / lambda_j| -// for the off-diagonal entries linking eigenvalue clusters j and j+1. -// (T&B Lecture 28, Theorem 28.2) -// -// Throws DimensionMismatchError if A is not square. -// Throws NonConvergenceError if convergence is not achieved within -// opts.max_iterations steps. -[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A, - QRIterationOptions opts = {}); - -// --------------------------------------------------------------------------- -// Wilkinson-shifted QR iteration -// --------------------------------------------------------------------------- -// -// Same outer loop as Stage 1, but each step applies a shift σ chosen as -// the eigenvalue of the bottom-right 2×2 block of A_{k-1} that is closest -// to the (n,n) entry, then unshifts after the QR step: -// -// factor (A_{k-1} - σI) = Q_k R_k -// set A_k = R_k Q_k + σI -// -// The Wilkinson shift (T&B Lecture 29; GVL §7.4.2): -// Given the bottom-right 2×2 block | a b | -// | b c | -// δ = (a - c) / 2 -// σ = c - sign(δ) * b² / (|δ| + sqrt(δ² + b²)) -// equivalently: the eigenvalue of the block closer to c. -// -// Convergence rate: typically cubic near a simple eigenvalue. -// (T&B Lecture 29; GVL §7.5.1) -// -// Same exceptions as eigenvalues_unshifted. -[[nodiscard]] QRIterationResult eigenvalues_shifted(const Matrix& A, - QRIterationOptions opts = {}); - -// --------------------------------------------------------------------------- -// Stage 3: Hessenberg reduction algorithm -// --------------------------------------------------------------------------- - -// Givens rotation G acting on rows/columns i and i+1: -// -// G = | c s | chosen so that G * [x; y]^T = [r; 0]^T -// | -s c | with c = x/r, s = y/r, r = hypot(x, y) -// -// T&B Lecture 10 (Givens rotations). -struct GivensRotation { - double c; // cos(theta) - double s; // sin(theta) - std::size_t i; // first row/column index (second is i+1) - - // Construct the rotation that maps [x, y]^T → [hypot(x,y), 0]^T. - // Returns the identity (c=1, s=0) when x == y == 0. - [[nodiscard]] static GivensRotation make(double x, double y, - std::size_t row_index); - - // Apply G from the left to rows i and i+1 of M, columns [col_start, n). - // M[i:i+2, col_start:] ← G * M[i:i+2, col_start:] - void apply_left(Matrix& M, std::size_t col_start = 0) const; - - // Apply G^T from the right to columns i and i+1 of M, rows [0, row_end). - // M[0:row_end, i:i+2] ← M[0:row_end, i:i+2] * G^T - void apply_right(Matrix& M, std::size_t row_end) const; -}; - -// Result of reducing A to upper Hessenberg form. -// H is upper Hessenberg: H(i,j) = 0 for all i > j+1. -// Q is orthogonal and A = Q H Q^T. -// Ref: GVL Algorithm 7.4.2; T&B Lecture 26. -struct HessenbergResult { - Matrix H; // upper Hessenberg similarity of A - Matrix Q; // accumulated orthogonal transformation -}; - -// Reduce A to upper Hessenberg form via Householder reflectors applied -// from both sides. Costs O(10n³/3) flops; done once before QR iteration. -// Ref: GVL §7.4.2 (Algorithm 7.4.2). -// -// Throws DimensionMismatchError if A is not square. -[[nodiscard]] HessenbergResult hessenberg_reduction(const Matrix& A); - -// Apply one shifted QR step to an upper Hessenberg matrix H in-place, -// using n-1 Givens rotations. Costs O(n²) vs O(n³) for Householder QR. -// H remains upper Hessenberg after the step. -// Ref: GVL §7.4.2; T&B Lecture 29. -void hessenberg_qr_step(Matrix& H, double sigma); - -// Full practical QR algorithm: -// 1. Reduce A to Hessenberg H = Q^T A Q (O(n³), done once). -// 2. Run Wilkinson-shifted QR on H using Givens steps (O(n²) each). -// Substantially faster than eigenvalues_shifted for n ≥ 50. -// Ref: T&B Lecture 29. -// -// Same exceptions as eigenvalues_unshifted. -[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A, - QRIterationOptions opts = {}); - -} // namespace linalg diff --git a/include/triangular_solve.hpp b/include/triangular_solve.hpp deleted file mode 100644 index 632f225..0000000 --- a/include/triangular_solve.hpp +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -#include <cstddef> - -#include "matrix.hpp" -#include "vector.hpp" - -namespace linalg { - -Vector forward_substitution( - const Matrix& lower, - const Vector& rhs, - double singular_tolerance = 1e-12, - bool unit_diagonal = false); - -Vector backward_substitution( - const Matrix& upper, - const Vector& rhs, - double singular_tolerance = 1e-12, - bool unit_diagonal = false); - -} // namespace linalg diff --git a/include/vector.hpp b/include/vector.hpp deleted file mode 100644 index 3e797d8..0000000 --- a/include/vector.hpp +++ /dev/null @@ -1,47 +0,0 @@ -#pragma once - -#include <cstddef> -#include <initializer_list> -#include <vector> - -namespace linalg { - -class Vector { -public: - Vector() = default; - explicit Vector(std::size_t n); - Vector(std::size_t n, double value); - Vector(std::initializer_list<double> values); - - [[nodiscard]] std::size_t size() const noexcept; - [[nodiscard]] bool empty() const noexcept; - - double& operator[](std::size_t i); - const double& operator[](std::size_t i) const; - - void fill(double value); - - double* data() noexcept; - const double* data() const noexcept; - - auto begin() noexcept { return data_.begin(); } - auto end() noexcept { return data_.end(); } - auto begin() const noexcept { return data_.begin(); } - auto end() const noexcept { return data_.end(); } - auto cbegin() const noexcept { return data_.cbegin(); } - auto cend() const noexcept { return data_.cend(); } - -private: - void check_index(std::size_t i) const; - - std::vector<double> data_; -}; - -Vector operator+(const Vector& lhs, const Vector& rhs); -Vector operator-(const Vector& lhs, const Vector& rhs); -Vector operator*(const Vector& v, double scalar); -Vector operator*(double scalar, const Vector& v); -Vector operator/(const Vector& v, double scalar); -double dot(const Vector& lhs, const Vector& rhs); - -} // namespace linalg diff --git a/src/linalgebra.cpp b/src/linalgebra.cpp new file mode 100644 index 0000000..302810f --- /dev/null +++ b/src/linalgebra.cpp @@ -0,0 +1,10 @@ +export module linalgebra; + +export import :error; +export import :vector; +export import :matrix; +export import :norms; +export import :triangular_solve; +export import :lu; +export import :qr; +export import :qr_iteration; diff --git a/include/linalg_error.hpp b/src/linalgebra_error.cpp index ac9b456..367c79d 100644 --- a/include/linalg_error.hpp +++ b/src/linalgebra_error.cpp @@ -1,9 +1,7 @@ -#pragma once +export module linalgebra:error; +import std; -#include <stdexcept> -#include <string> - -namespace linalg { +export namespace linalgebra { class LinAlgError : public std::runtime_error { public: @@ -28,4 +26,4 @@ public: : LinAlgError(message) {} }; -} // namespace linalg +} // namespace linalgebra @@ -1,14 +1,26 @@ -#include "lu.hpp" +export module linalgebra:lu; +import std; +import :error; +import :vector; +import :matrix; +import :triangular_solve; -#include <algorithm> -#include <cmath> -#include <numeric> -#include <sstream> +export namespace linalgebra { -#include "linalg_error.hpp" -#include "triangular_solve.hpp" +struct LUResult { + Matrix L; + Matrix U; + std::vector<std::size_t> perm; + int sign; +}; -namespace linalg { +LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12); + +Vector lu_solve(const LUResult& lu, const Vector& b); + +} // namespace linalgebra + +namespace linalgebra { LUResult lu_factor(const Matrix& A, double singular_tolerance) { if (A.rows() != A.cols()) { @@ -33,7 +45,6 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { int sign = 1; for (std::size_t k = 0; k < n; ++k) { - // ---- Partial pivoting: find row with largest magnitude in column k ---- std::size_t pivot_row = k; double max_val = std::abs(work(k, k)); for (std::size_t i = k + 1; i < n; ++i) { @@ -55,7 +66,6 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { sign = -sign; } - // ---- Singularity check ---- if (std::abs(work(k, k)) <= singular_tolerance) { std::ostringstream oss; oss << "lu_factor: near-zero pivot " << work(k, k) << " at step " << k @@ -63,12 +73,10 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) { throw SingularMatrixError(oss.str()); } - // ---- Record U row k ---- for (std::size_t j = k; j < n; ++j) { U(k, j) = work(k, j); } - // ---- Compute multipliers and eliminate below pivot ---- for (std::size_t i = k + 1; i < n; ++i) { L(i, k) = work(i, k) / work(k, k); for (std::size_t j = k + 1; j < n; ++j) { @@ -101,4 +109,4 @@ Vector lu_solve(const LUResult& lu, const Vector& b) { return backward_substitution(lu.U, y); } -} // namespace linalg +} // namespace linalgebra diff --git a/src/matrix.cpp b/src/matrix.cpp index c27264b..75fcb07 100644 --- a/src/matrix.cpp +++ b/src/matrix.cpp @@ -1,120 +1,76 @@ -#include "matrix.hpp" -#include "linalg_error.hpp" +module; -#include <algorithm> -#include <cstddef> -#include <sstream> -#include <stdexcept> - -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \ - (defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)) -#include <immintrin.h> -#endif - -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC) +#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) #include <arm_neon.h> #endif -namespace linalg { +export module linalgebra:matrix; +import std; +import :error; +import :vector; -namespace { +export namespace linalgebra { -void check_same_shape(const Matrix& lhs, const Matrix& rhs, const char* operation) { - if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) { - std::ostringstream oss; - oss << operation << " requires equal matrix shapes, got " << lhs.rows() << "x" << lhs.cols() << " and " << rhs.rows() << "x" << rhs.cols(); - throw DimensionMismatchError(oss.str()); - } -} +class Matrix { +public: + Matrix() = default; + Matrix(std::size_t rows, std::size_t cols); + Matrix(std::size_t rows, std::size_t cols, double value); + Matrix(std::initializer_list<std::initializer_list<double>> values); -double dot_product_scalar(const double* lhs, const double* rhs, std::size_t count) { - double sum = 0.0; - for (std::size_t i = 0; i < count; ++i) { - sum += lhs[i] * rhs[i]; - } - return sum; -} - -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__) -double horizontal_sum(__m512d values) { - alignas(64) double lanes[8]; - _mm512_store_pd(lanes, values); - double sum = 0.0; - for (double lane : lanes) { - sum += lane; - } - return sum; -} + [[nodiscard]] std::size_t rows() const noexcept; + [[nodiscard]] std::size_t cols() const noexcept; + [[nodiscard]] bool empty() const noexcept; -double dot_product_avx512(const double* lhs, const double* rhs, std::size_t count) { - std::size_t i = 0; - __m512d acc0 = _mm512_setzero_pd(); - __m512d acc1 = _mm512_setzero_pd(); + double& operator()(std::size_t i, std::size_t j); + const double& operator()(std::size_t i, std::size_t j) const; - for (; i + 15 < count; i += 16) { - const __m512d lhs0 = _mm512_loadu_pd(lhs + i); - const __m512d rhs0 = _mm512_loadu_pd(rhs + i); - const __m512d lhs1 = _mm512_loadu_pd(lhs + i + 8); - const __m512d rhs1 = _mm512_loadu_pd(rhs + i + 8); + void fill(double value); - acc0 = _mm512_add_pd(acc0, _mm512_mul_pd(lhs0, rhs0)); - acc1 = _mm512_add_pd(acc1, _mm512_mul_pd(lhs1, rhs1)); - } + double* data() noexcept; + const double* data() const noexcept; - return horizontal_sum(acc0) + horizontal_sum(acc1) + - dot_product_scalar(lhs + i, rhs + i, count - i); -} -#endif + static Matrix identity(std::size_t n); + static Matrix zeros(std::size_t rows, std::size_t cols); -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__) -double horizontal_sum(__m256d values) { - alignas(32) double lanes[4]; - _mm256_store_pd(lanes, values); - return lanes[0] + lanes[1] + lanes[2] + lanes[3]; -} +private: + [[nodiscard]] std::size_t index(std::size_t i, std::size_t j) const; + void check_bounds(std::size_t i, std::size_t j) const; -double dot_product_avx2(const double* lhs, const double* rhs, std::size_t count) { - std::size_t i = 0; - __m256d acc0 = _mm256_setzero_pd(); - __m256d acc1 = _mm256_setzero_pd(); + std::size_t rows_ = 0; + std::size_t cols_ = 0; + std::vector<double> data_; +}; - for (; i + 7 < count; i += 8) { - const __m256d lhs0 = _mm256_loadu_pd(lhs + i); - const __m256d rhs0 = _mm256_loadu_pd(rhs + i); - const __m256d lhs1 = _mm256_loadu_pd(lhs + i + 4); - const __m256d rhs1 = _mm256_loadu_pd(rhs + i + 4); +Matrix transpose(const Matrix& matrix); +Matrix operator+(const Matrix& lhs, const Matrix& rhs); +Matrix operator-(const Matrix& lhs, const Matrix& rhs); +Vector operator*(const Matrix& matrix, const Vector& vector); +Matrix operator*(const Matrix& lhs, const Matrix& rhs); - acc0 = _mm256_add_pd(acc0, _mm256_mul_pd(lhs0, rhs0)); - acc1 = _mm256_add_pd(acc1, _mm256_mul_pd(lhs1, rhs1)); - } +} // namespace linalgebra - return horizontal_sum(acc0) + horizontal_sum(acc1) + - dot_product_scalar(lhs + i, rhs + i, count - i); -} -#endif +namespace { -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) && !defined(__AVX2__) -double horizontal_sum(__m256d values) { - alignas(32) double lanes[4]; - _mm256_store_pd(lanes, values); - return lanes[0] + lanes[1] + lanes[2] + lanes[3]; +void check_same_shape(const linalgebra::Matrix& lhs, const linalgebra::Matrix& rhs, + const char* operation) { + if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) { + std::ostringstream oss; + oss << operation << " requires equal matrix shapes, got " << lhs.rows() << "x" + << lhs.cols() << " and " << rhs.rows() << "x" << rhs.cols(); + throw linalgebra::DimensionMismatchError(oss.str()); + } } -double dot_product_avx(const double* lhs, const double* rhs, std::size_t count) { - std::size_t i = 0; - __m256d acc = _mm256_setzero_pd(); - - for (; i + 3 < count; i += 4) { - const __m256d lhs_values = _mm256_loadu_pd(lhs + i); - const __m256d rhs_values = _mm256_loadu_pd(rhs + i); - acc = _mm256_add_pd(acc, _mm256_mul_pd(lhs_values, rhs_values)); +double dot_product_scalar(const double* lhs, const double* rhs, std::size_t count) { + double sum = 0.0; + for (std::size_t i = 0; i < count; ++i) { + sum += lhs[i] * rhs[i]; } - - return horizontal_sum(acc) + dot_product_scalar(lhs + i, rhs + i, count - i); + return sum; } -#endif -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC) +#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) double horizontal_sum(float64x2_t values) { return vgetq_lane_f64(values, 0) + vgetq_lane_f64(values, 1); } @@ -140,13 +96,7 @@ double dot_product_neon(const double* lhs, const double* rhs, std::size_t count) #endif double dot_product_simd(const double* lhs, const double* rhs, std::size_t count) { -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__) - return dot_product_avx512(lhs, rhs, count); -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__) - return dot_product_avx2(lhs, rhs, count); -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) && !defined(__AVX2__) - return dot_product_avx(lhs, rhs, count); -#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC) +#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) return dot_product_neon(lhs, rhs, count); #else return dot_product_scalar(lhs, rhs, count); @@ -155,13 +105,16 @@ double dot_product_simd(const double* lhs, const double* rhs, std::size_t count) } // namespace +namespace linalgebra { + Matrix::Matrix(std::size_t rows, std::size_t cols) : rows_(rows), cols_(cols), data_(rows * cols) {} Matrix::Matrix(std::size_t rows, std::size_t cols, double value) : rows_(rows), cols_(cols), data_(rows * cols, value) {} -Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) : rows_(values.size()) { +Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) + : rows_(values.size()) { if (rows_ == 0) { cols_ = 0; return; @@ -178,7 +131,6 @@ Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) : ro } } - std::size_t Matrix::rows() const noexcept { return rows_; } std::size_t Matrix::cols() const noexcept { return cols_; } @@ -298,4 +250,4 @@ Matrix operator*(const Matrix& lhs, const Matrix& rhs) { return result; } -} // namespace linalg +} // namespace linalgebra diff --git a/src/norms.cpp b/src/norms.cpp index cafee69..7a40f1f 100644 --- a/src/norms.cpp +++ b/src/norms.cpp @@ -1,9 +1,15 @@ -#include "norms.hpp" +export module linalgebra:norms; +import std; +import :vector; -#include <cmath> +export namespace linalgebra { -namespace linalg { +double norm2(const Vector& vector); + +} // namespace linalgebra + +namespace linalgebra { double norm2(const Vector& vector) { return std::sqrt(dot(vector, vector)); } -} // namespace linalg +} // namespace linalgebra @@ -1,23 +1,55 @@ -#include "qr.hpp" +export module linalgebra:qr; +import std; +import :error; +import :vector; +import :matrix; -#include <cmath> -#include <sstream> +export namespace linalgebra { -#include "linalg_error.hpp" +struct QRResult { + Matrix Q; + Matrix R; +}; -namespace linalg { +// Classical Gram-Schmidt. +// Mathematically natural but numerically fragile: orthogonality of Q +// degrades rapidly on ill-conditioned inputs. +// Provided for comparison — prefer modified_gs or householder in practice. +// +// Throws DimensionMismatchError if rows < cols. +// Throws SingularMatrixError if a column is (nearly) linearly dependent. +QRResult qr_classical_gs(const Matrix& A, double zero_tolerance = 1e-14); + +// Modified Gram-Schmidt. +// Subtracts each projection immediately on the running vector rather than +// on the original column. Algebraically equivalent to classical GS but +// numerically much better — round-off stays local instead of accumulating. +// +// Same exceptions as classical GS. +QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14); + +// Householder QR. +// Applies a sequence of orthogonal reflections to zero out below-diagonal +// entries column by column. Backward-stable and the standard choice for +// dense QR. Works correctly on rank-deficient matrices (zero pivots +// produce zero diagonal entries in R without throwing). +// +// Throws DimensionMismatchError if rows < cols. +QRResult qr_householder(const Matrix& A); + +} // namespace linalgebra namespace { -void require_tall(const Matrix& A, const char* name) { +void require_tall(const linalgebra::Matrix& A, const char* name) { if (A.rows() < A.cols()) { std::ostringstream oss; oss << name << " requires rows >= cols, got " << A.rows() << "x" << A.cols(); - throw DimensionMismatchError(oss.str()); + throw linalgebra::DimensionMismatchError(oss.str()); } } -double col_norm(const Matrix& M, std::size_t j) { +double col_norm(const linalgebra::Matrix& M, std::size_t j) { double s = 0.0; for (std::size_t i = 0; i < M.rows(); ++i) { s += M(i, j) * M(i, j); @@ -25,7 +57,8 @@ double col_norm(const Matrix& M, std::size_t j) { return std::sqrt(s); } -double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) { +double col_dot(const linalgebra::Matrix& M, std::size_t j, + const linalgebra::Matrix& N, std::size_t k) { double s = 0.0; for (std::size_t i = 0; i < M.rows(); ++i) { s += M(i, j) * N(i, k); @@ -35,7 +68,7 @@ double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) { } // namespace -// --- Gram-Schmidt --- +namespace linalgebra { QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { require_tall(A, "qr_classical_gs"); @@ -49,7 +82,7 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j); for (std::size_t k = 0; k < j; ++k) { - R(k, j) = col_dot(A, j, Q, k); // <a_j, q_k> + R(k, j) = col_dot(A, j, Q, k); for (std::size_t i = 0; i < m; ++i) { Q(i, j) -= R(k, j) * Q(i, k); } @@ -69,8 +102,6 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) { return QRResult{std::move(Q), std::move(R)}; } -// --- Modified Gram-Schmidt --- - QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) { require_tall(A, "qr_modified_gs"); const std::size_t m = A.rows(); @@ -83,7 +114,7 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) { for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j); for (std::size_t k = 0; k < j; ++k) { - R(k, j) = col_dot(Q, j, Q, k); // <v_running, q_k> + R(k, j) = col_dot(Q, j, Q, k); for (std::size_t i = 0; i < m; ++i) { Q(i, j) -= R(k, j) * Q(i, k); } @@ -103,20 +134,16 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) { return QRResult{std::move(Q), std::move(R)}; } -// --- Householder QR --- - QRResult qr_householder(const Matrix& A) { require_tall(A, "qr_householder"); const std::size_t m = A.rows(); const std::size_t n = A.cols(); - // Will become R. Matrix work = A; - 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 + const std::size_t p = m - k; std::vector<double> u(p); for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k); @@ -139,19 +166,17 @@ QRResult qr_householder(const Matrix& A) { }(); const double tau = 2.0 / utu; - // Apply H_k to work[k:, k:n] for (std::size_t j = k; j < n; ++j) { - double dot = 0.0; - for (std::size_t i = 0; i < p; ++i) dot += u[i] * work(k + i, j); - const double coeff = tau * dot; + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * work(k + i, j); + const double coeff = tau * d; for (std::size_t i = 0; i < p; ++i) work(k + i, j) -= coeff * u[i]; } - // Apply H_k to Q_full[k:, 0:m] for (std::size_t j = 0; j < m; ++j) { - double dot = 0.0; - for (std::size_t i = 0; i < p; ++i) dot += u[i] * Q_full(k + i, j); - const double coeff = tau * dot; + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * Q_full(k + i, j); + const double coeff = tau * d; for (std::size_t i = 0; i < p; ++i) Q_full(k + i, j) -= coeff * u[i]; } } @@ -169,4 +194,4 @@ QRResult qr_householder(const Matrix& A) { return QRResult{std::move(Q), std::move(R)}; } -} // namespace linalg +} // namespace linalgebra diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp index e476781..64fd25f 100644 --- a/src/qr_iteration.cpp +++ b/src/qr_iteration.cpp @@ -1,73 +1,99 @@ -#include "qr_iteration.hpp" - +module; #include <cassert> -#include <cmath> -#include <sstream> -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "qr.hpp" -#include "vector.hpp" +export module linalgebra:qr_iteration; +import std; +import :error; +import :vector; +import :matrix; +import :qr; // References used throughout this file: // T&B — Trefethen & Bau, "Numerical Linear Algebra" // GVL — Golub & Van Loan, "Matrix Computations" 4th ed. -namespace linalg { +export namespace linalgebra { -namespace { +// --------------------------------------------------------------------------- +// Options +// --------------------------------------------------------------------------- + +struct QRIterationOptions { + double tolerance = 1e-10; + int max_iterations = 1000; + bool track_convergence = false; +}; + +struct QRIterationResult { + Vector eigenvalues_real; + Vector eigenvalues_imag; + int iterations = 0; + // std::vector is used here because linalgebra::Vector has no push_back; + // convergence_history is a plain time-series container, not a math object. + std::vector<double> convergence_history; +}; + +[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A, + QRIterationOptions opts = {}); + +[[nodiscard]] QRIterationResult eigenvalues_shifted(const Matrix& A, + QRIterationOptions opts = {}); + +// Givens rotation G acting on rows/columns i and i+1: +// +// G = | c s | chosen so that G * [x; y]^T = [r; 0]^T +// | -s c | with c = x/r, s = y/r, r = hypot(x, y) +struct GivensRotation { + double c; + double s; + std::size_t i; + + [[nodiscard]] static GivensRotation make(double x, double y, std::size_t row_index); + void apply_left(Matrix& M, std::size_t col_start = 0) const; + void apply_right(Matrix& M, std::size_t row_end) const; +}; +struct HessenbergResult { + Matrix H; + Matrix Q; +}; -// 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. -// Ref: T&B §28; used as the convergence criterion in Algorithm 28.1. -double lower_triangle_norm(const Matrix& A) { +[[nodiscard]] HessenbergResult hessenberg_reduction(const Matrix& A); + +void hessenberg_qr_step(Matrix& H, double sigma); + +[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A, + QRIterationOptions opts = {}); + +} // namespace linalgebra + +namespace { + +double lower_triangle_norm(const linalgebra::Matrix& A) { const std::size_t n = A.rows(); double s = 0.0; - for (std::size_t i = 1; i < n; ++i) // row 1 .. n-1 - for (std::size_t j = 0; j < i; ++j) // col 0 .. i-1 (strict lower) + for (std::size_t i = 1; i < n; ++i) + for (std::size_t j = 0; j < i; ++j) s += A(i, j) * A(i, j); return std::sqrt(s); } -// Extract eigenvalues from a quasi-upper-triangular matrix (real Schur form). -// -// Scans the diagonal from top-left to bottom-right. At each position i: -// — |A(i+1, i)| < tol → 1×1 block: real eigenvalue A(i,i), imag = 0. -// — otherwise → 2×2 block [A(i..i+1, i..i+1)]: eigenvalues via -// quadratic formula. When the discriminant is -// negative the result is a complex-conjugate pair, -// stored as (re, +im) and (re, -im) in the real -// and imaginary part Vectors. -// -// Fills positions 0..n-1 of `real_out` and `imag_out` (pre-sized to n). -// -// Ref: T&B Lecture 28; GVL §7.4.1. -void extract_eigenvalues(const Matrix& T, double tol, - Vector& real_out, Vector& imag_out) { +void extract_eigenvalues(const linalgebra::Matrix& T, double tol, + linalgebra::Vector& real_out, linalgebra::Vector& imag_out) { const std::size_t n = T.rows(); - std::size_t out = 0; - std::size_t i = 0; + std::size_t out = 0; + std::size_t i = 0; while (i < n) { const bool is_last = (i + 1 == n); const bool sub_small = is_last || (std::abs(T(i + 1, i)) < tol); if (sub_small) { - // 1×1 block: real eigenvalue. real_out[out] = T(i, i); imag_out[out] = 0.0; ++out; ++i; } else { - // 2×2 block: - // | a b | - // | c d | - // 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); @@ -76,15 +102,12 @@ void extract_eigenvalues(const Matrix& T, double tol, const double disc = (a - d) * (a - d) + 4.0 * b * c; if (disc >= 0.0) { - // 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); imag_out[out] = 0.0; real_out[out + 1] = 0.5 * (tr - sq); imag_out[out + 1] = 0.0; } else { - // Complex-conjugate pair: real part ± imaginary part. const double re = 0.5 * tr; const double im = 0.5 * std::sqrt(-disc); real_out[out] = re; @@ -100,37 +123,32 @@ void extract_eigenvalues(const Matrix& T, double tol, assert(out == n); } -void require_square(const Matrix& A, const char* fname) { +void require_square(const linalgebra::Matrix& A, const char* fname) { if (A.rows() != A.cols()) { std::ostringstream oss; oss << fname << ": requires a square matrix, got " << A.rows() << "x" << A.cols(); - throw DimensionMismatchError(oss.str()); + throw linalgebra::DimensionMismatchError(oss.str()); } } +double wilkinson_shift(const linalgebra::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); + const double d = A(n - 1, n - 1); + const double delta = 0.5 * (a - d); + const double denom = std::abs(delta) + std::hypot(delta, b); + if (denom == 0.0) return d; + const double sgn = (delta >= 0.0) ? 1.0 : -1.0; + return d - sgn * (b * b) / denom; +} + } // namespace -// --- Unshifted QR iteration --- -// -// Each step performs an orthogonal similarity transformation: -// A_{k-1} = Q_k R_k (Householder QR; backward-stable) -// A_k = R_k Q_k = Q_k^T A_{k-1} Q_k -// -// Similarity preserves eigenvalues (GVL §7.3.1, Theorem 7.3.1). -// The iterates converge to the real Schur form: a quasi-upper-triangular -// matrix whose 1×1 blocks give real eigenvalues and 2×2 blocks give -// complex-conjugate pairs. -// -// Convergence rate: linear, with per-step reduction factor -// |lambda_{j+1} / lambda_j| for the (j, j+1) coupling. -// (T&B Lecture 28, Theorem 28.2; GVL §7.3.2) -// -// Each iteration costs O(n^3) due to full Householder QR; Hessenberg -// reduction (Stage 3) reduces subsequent steps to O(n^2). +namespace linalgebra { -QRIterationResult eigenvalues_unshifted(const Matrix& A, - QRIterationOptions opts) { +QRIterationResult eigenvalues_unshifted(const Matrix& A, QRIterationOptions opts) { require_square(A, "eigenvalues_unshifted"); const std::size_t n = A.rows(); @@ -153,13 +171,9 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, Matrix Ak = A; for (int k = 0; k < opts.max_iterations; ++k) { - // Factor A_{k-1} = Q R using backward-stable Householder reflections. const QRResult qr = qr_householder(Ak); - - // A_k = R Q (orthogonal similarity: Q^T A_{k-1} Q) Ak = qr.R * qr.Q; - // --- Convergence check --- const double lower_norm = lower_triangle_norm(Ak); if (opts.track_convergence) { @@ -184,52 +198,6 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, throw NonConvergenceError(oss.str()); } -// --- Wilkinson-shifted QR iteration --- -// -// The Wilkinson shift is the eigenvalue of the bottom-right 2×2 block -// | a b | -// | c d | -// that is closest to d (the trailing diagonal entry). -// -// Exact eigenvalue formula: μ_{1,2} = (a+d)/2 ± sqrt(((a-d)/2)² + b·c) -// We pick the one with |μ - d| smaller. -// -// When the discriminant is negative (complex eigenvalues), fall back to σ = d -// (Rayleigh quotient shift), which still accelerates convergence. -// -// Ref: T&B Lecture 29; GVL §7.4.2. - -namespace { - -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 delta = 0.5 * (a - d); - const double denom = std::abs(delta) + std::hypot(delta, b); - if (denom == 0.0) return d; - const double sgn = (delta >= 0.0) ? 1.0 : -1.0; - return d - sgn * (b * b) / denom; -} - -} // namespace - -// eigenvalues_shifted — Wilkinson-shifted QR with trailing deflation. -// -// After each QR step we check whether the trailing subdiagonal entry of the -// active block is negligible (relative criterion: GVL §7.4.1). If so, the -// bottom diagonal entry is accepted as a converged eigenvalue and the active -// subproblem shrinks by one. This "trailing deflation" enables the cubic -// convergence promised by the Wilkinson shift to compound across successive -// eigenvalues rather than stalling on the full lower-triangle norm. -// -// When the active size reaches 2 we extract both eigenvalues analytically -// from the 2×2 block (handling real and complex-conjugate pairs) rather than -// continuing to iterate. For symmetric inputs this is always a real pair. -// -// Ref: GVL §7.5.1; T&B Lecture 29. - QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) { require_square(A, "eigenvalues_shifted"); const std::size_t n = A.rows(); @@ -250,7 +218,7 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) Matrix Ak = A; std::size_t n_found = n; - std::size_t active = n; // live subproblem is rows/cols 0..active-1 + std::size_t active = n; auto store_real = [&](double re) { --n_found; @@ -281,16 +249,14 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) }; for (int k = 0; k < opts.max_iterations; ++k) { - // --- Deflation sweep --- while (active >= 2) { const double sub = std::abs(Ak(active - 1, active - 2)); const double scale = std::abs(Ak(active - 2, active - 2)) + std::abs(Ak(active - 1, active - 1)); - // Relative + absolute floor tolerance (GVL §7.4.1). const double deflation_tol = opts.tolerance * (scale > 0.0 ? scale : 1.0); if (sub > deflation_tol) break; - Ak(active - 1, active - 2) = 0.0; // enforce exact zero + Ak(active - 1, active - 2) = 0.0; store_real(Ak(active - 1, active - 1)); --active; } @@ -299,7 +265,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) if (active == 1) { store_real(Ak(0, 0)); active = 0; break; } if (active == 2) { close_2x2(); break; } - // --- Wilkinson-shifted QR step on the active × active subblock --- Matrix sub_mat(active, active); for (std::size_t i = 0; i < active; ++i) for (std::size_t j = 0; j < active; ++j) @@ -331,8 +296,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) return result; } -// --- Givens rotation --- - GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) { const double r = std::hypot(x, y); if (r == 0.0) return {1.0, 0.0, row_index}; @@ -340,9 +303,6 @@ GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) { } void GivensRotation::apply_left(Matrix& M, std::size_t col_start) const { - // Rows i and i+1, columns col_start..n-1. - // [ c s] [x] [cx + sy] - // [-s c] [y] = [-sx + cy] for (std::size_t j = col_start; j < M.cols(); ++j) { const double xi = M(i, j); const double xi1 = M(i + 1, j); @@ -352,10 +312,6 @@ void GivensRotation::apply_left(Matrix& M, std::size_t col_start) const { } void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const { - // Columns i and i+1, rows 0..row_end-1. - // M * G^T where G^T = [c -s; s c]: - // new col i = c * old_i + s * old_{i+1} - // new col i+1 = -s * old_i + c * old_{i+1} for (std::size_t j = 0; j < row_end; ++j) { const double xi = M(j, i); const double xi1 = M(j, i + 1); @@ -364,16 +320,6 @@ void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const { } } -// --- 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] -// Apply from right: A[0:n, k+1:n] ← A[0:n, k+1:n] * H_k -// Accumulate Q: Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k -// -// H_k is never formed explicitly; applied via rank-1 update with tau = 2/uᵀu. -// Ref: GVL §7.4.2 (Algorithm 7.4.2). - HessenbergResult hessenberg_reduction(const Matrix& A) { require_square(A, "hessenberg_reduction"); const std::size_t n = A.rows(); @@ -382,10 +328,9 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { Matrix Q = Matrix::identity(n); for (std::size_t k = 0; k + 2 <= n; ++k) { - const std::size_t p = n - k - 1; // p = n - (k+1) + const std::size_t p = n - k - 1; if (p == 0) break; - // Build Householder vector u from H[k+1:n, k]. std::vector<double> u(p); for (std::size_t i = 0; i < p; ++i) u[i] = H(k + 1 + i, k); @@ -402,27 +347,24 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { for (double v : u) utu += v * v; const double tau = 2.0 / utu; - // Apply H_k from the LEFT to H[k+1:n, k:n]. for (std::size_t j = k; j < n; ++j) { - double dot = 0.0; - for (std::size_t i = 0; i < p; ++i) dot += u[i] * H(k + 1 + i, j); - const double coeff = tau * dot; + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * H(k + 1 + i, j); + const double coeff = tau * d; for (std::size_t i = 0; i < p; ++i) H(k + 1 + i, j) -= coeff * u[i]; } - // Apply H_k from the RIGHT to H[0:n, k+1:n]. for (std::size_t j = 0; j < n; ++j) { - double dot = 0.0; - for (std::size_t i = 0; i < p; ++i) dot += H(j, k + 1 + i) * u[i]; - const double coeff = tau * dot; + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += H(j, k + 1 + i) * u[i]; + const double coeff = tau * d; for (std::size_t i = 0; i < p; ++i) H(j, k + 1 + i) -= coeff * u[i]; } - // Accumulate Q: Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k. for (std::size_t j = 0; j < n; ++j) { - double dot = 0.0; - for (std::size_t i = 0; i < p; ++i) dot += Q(j, k + 1 + i) * u[i]; - const double coeff = tau * dot; + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += Q(j, k + 1 + i) * u[i]; + const double coeff = tau * d; for (std::size_t i = 0; i < p; ++i) Q(j, k + 1 + i) -= coeff * u[i]; } @@ -432,20 +374,6 @@ HessenbergResult hessenberg_reduction(const Matrix& A) { return HessenbergResult{std::move(H), std::move(Q)}; } -// --- Hessenberg QR step via Givens rotations --- -// -// One shifted QR step on the upper Hessenberg matrix H: -// 1. Shift: H ← H - σI. -// 2. For k = 0..n-2: compute G_k = Givens(H(k,k), H(k+1,k)); -// apply G_k from left to rows k,k+1 of H, -// starting from column k (Hessenberg: H(k+1,j)=0, j<k). -// 3. For k = 0..n-2: apply G_k^T from right to cols k,k+1 of H, -// up to row k+2 (exploits upper-triangular structure). -// 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. - void hessenberg_qr_step(Matrix& H, double sigma) { const std::size_t n = H.rows(); @@ -455,10 +383,7 @@ void hessenberg_qr_step(Matrix& H, double sigma) { gs.reserve(n - 1); for (std::size_t k = 0; k + 1 < n; ++k) { - // Eliminate H(k+1, k) via a rotation on rows k and k+1. GivensRotation g = GivensRotation::make(H(k, k), H(k + 1, k), k); - // Left application: rows k, k+1; columns k..n-1. - // (Hessenberg: H(k+1, j) = 0 for j < k, so starting from col k is exact.) g.apply_left(H, k); gs.push_back(g); } @@ -467,23 +392,10 @@ void hessenberg_qr_step(Matrix& H, double sigma) { gs[k].apply_right(H, std::min(k + 2, n)); } - // Unshift. for (std::size_t j = 0; j < n; ++j) H(j, j) += sigma; } -// --- 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 -// (O(n³)). After Hessenberg reduction the matrix stays Hessenberg throughout, -// so the O(n²) per-step cost applies for every step after the one-time O(n³) -// reduction. Total cost is thus O(n³) + O(iterations · n²), which beats -// eigenvalues_shifted's O(iterations · n³) for large n. -// -// Ref: GVL §7.4.2; T&B Lecture 29. - -QRIterationResult eigenvalues_hessenberg(const Matrix& A, - QRIterationOptions opts) { +QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opts) { require_square(A, "eigenvalues_hessenberg"); const std::size_t n = A.rows(); @@ -535,7 +447,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, }; for (int k = 0; k < opts.max_iterations; ++k) { - // --- Deflation sweep --- while (active >= 2) { const double sub = std::abs(H(active - 1, active - 2)); const double scale = std::abs(H(active - 2, active - 2)) @@ -552,7 +463,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, if (active == 1) { store_real(H(0, 0)); active = 0; break; } if (active == 2) { close_2x2(); break; } - // Wilkinson shift from trailing 2×2 of the active block. 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); @@ -561,7 +471,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, 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. Matrix sub_H(active, active); for (std::size_t ii = 0; ii < active; ++ii) for (std::size_t jj = 0; jj < active; ++jj) @@ -593,4 +502,4 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, return result; } -} // namespace linalg +} // namespace linalgebra diff --git a/src/triangular_solve.cpp b/src/triangular_solve.cpp index 6a5e8ad..58dda6b 100644 --- a/src/triangular_solve.cpp +++ b/src/triangular_solve.cpp @@ -1,28 +1,41 @@ -#include "triangular_solve.hpp" +export module linalgebra:triangular_solve; +import std; +import :error; +import :vector; +import :matrix; -#include "linalg_error.hpp" +export namespace linalgebra { -#include <cmath> -#include <sstream> -#include <stdexcept> +Vector forward_substitution( + const Matrix& lower, + const Vector& rhs, + double singular_tolerance = 1e-12, + bool unit_diagonal = false); + +Vector backward_substitution( + const Matrix& upper, + const Vector& rhs, + double singular_tolerance = 1e-12, + bool unit_diagonal = false); -namespace linalg { +} // namespace linalgebra namespace { -void validate_square_system(const Matrix& matrix, const Vector& rhs, const char* operation) { +void validate_square_system(const linalgebra::Matrix& matrix, const linalgebra::Vector& rhs, + const char* operation) { if (matrix.rows() != matrix.cols()) { std::ostringstream oss; oss << operation << " requires a square matrix, got " << matrix.rows() << "x" << matrix.cols(); - throw DimensionMismatchError(oss.str()); + throw linalgebra::DimensionMismatchError(oss.str()); } if (matrix.rows() != rhs.size()) { std::ostringstream oss; oss << operation << " requires matrix dimension to match rhs size, got " << matrix.rows() << " and " << rhs.size(); - throw DimensionMismatchError(oss.str()); + throw linalgebra::DimensionMismatchError(oss.str()); } } @@ -32,10 +45,8 @@ void validate_tolerance(double singular_tolerance) { } } -void validate_lower_triangular( - const Matrix& lower, - double singular_tolerance, - bool unit_diagonal) { +void validate_lower_triangular(const linalgebra::Matrix& lower, double singular_tolerance, + bool unit_diagonal) { for (std::size_t i = 0; i < lower.rows(); ++i) { for (std::size_t j = i + 1; j < lower.cols(); ++j) { if (std::abs(lower(i, j)) > singular_tolerance) { @@ -45,16 +56,14 @@ void validate_lower_triangular( } if (!unit_diagonal && std::abs(lower(i, i)) <= singular_tolerance) { - throw SingularMatrixError( + throw linalgebra::SingularMatrixError( "Forward substitution encountered a zero or tiny diagonal entry"); } } } -void validate_upper_triangular( - const Matrix& upper, - double singular_tolerance, - bool unit_diagonal) { +void validate_upper_triangular(const linalgebra::Matrix& upper, double singular_tolerance, + bool unit_diagonal) { for (std::size_t i = 0; i < upper.rows(); ++i) { for (std::size_t j = 0; j < i; ++j) { if (std::abs(upper(i, j)) > singular_tolerance) { @@ -64,7 +73,7 @@ void validate_upper_triangular( } if (!unit_diagonal && std::abs(upper(i, i)) <= singular_tolerance) { - throw SingularMatrixError( + throw linalgebra::SingularMatrixError( "Backward substitution encountered a negligible diagonal entry"); } } @@ -72,11 +81,10 @@ void validate_upper_triangular( } // namespace -Vector forward_substitution( - const Matrix& lower, - const Vector& rhs, - double singular_tolerance, - bool unit_diagonal) { +namespace linalgebra { + +Vector forward_substitution(const Matrix& lower, const Vector& rhs, + double singular_tolerance, bool unit_diagonal) { validate_tolerance(singular_tolerance); validate_square_system(lower, rhs, "Forward substitution"); validate_lower_triangular(lower, singular_tolerance, unit_diagonal); @@ -98,11 +106,8 @@ Vector forward_substitution( return solution; } -Vector backward_substitution( - const Matrix& upper, - const Vector& rhs, - double singular_tolerance, - bool unit_diagonal) { +Vector backward_substitution(const Matrix& upper, const Vector& rhs, + double singular_tolerance, bool unit_diagonal) { validate_tolerance(singular_tolerance); validate_square_system(upper, rhs, "Backward substitution"); validate_upper_triangular(upper, singular_tolerance, unit_diagonal); @@ -125,4 +130,4 @@ Vector backward_substitution( return solution; } -} // namespace linalg +} // namespace linalgebra diff --git a/src/vector.cpp b/src/vector.cpp index 8c4a8c4..ed86f88 100644 --- a/src/vector.cpp +++ b/src/vector.cpp @@ -1,26 +1,65 @@ -#include "vector.hpp" -#include "linalg_error.hpp" +export module linalgebra:vector; +import std; +import :error; -#include <algorithm> -#include <numeric> -#include <sstream> -#include <stdexcept> +export namespace linalgebra { -namespace linalg { +class Vector { +public: + Vector() = default; + explicit Vector(std::size_t n); + Vector(std::size_t n, double value); + Vector(std::initializer_list<double> values); + + [[nodiscard]] std::size_t size() const noexcept; + [[nodiscard]] bool empty() const noexcept; + + double& operator[](std::size_t i); + const double& operator[](std::size_t i) const; + + void fill(double value); + + double* data() noexcept; + const double* data() const noexcept; + + auto begin() noexcept { return data_.begin(); } + auto end() noexcept { return data_.end(); } + auto begin() const noexcept { return data_.begin(); } + auto end() const noexcept { return data_.end(); } + auto cbegin() const noexcept { return data_.cbegin(); } + auto cend() const noexcept { return data_.cend(); } + +private: + void check_index(std::size_t i) const; + + std::vector<double> data_; +}; + +Vector operator+(const Vector& lhs, const Vector& rhs); +Vector operator-(const Vector& lhs, const Vector& rhs); +Vector operator*(const Vector& v, double scalar); +Vector operator*(double scalar, const Vector& v); +Vector operator/(const Vector& v, double scalar); +double dot(const Vector& lhs, const Vector& rhs); + +} // namespace linalgebra namespace { -void check_same_size(const Vector& lhs, const Vector& rhs, const char* operation) { +void check_same_size(const linalgebra::Vector& lhs, const linalgebra::Vector& rhs, + const char* operation) { if (lhs.size() != rhs.size()) { std::ostringstream oss; oss << operation << " requires equal vector sizes, got " << lhs.size() << " and " << rhs.size(); - throw DimensionMismatchError(oss.str()); + throw linalgebra::DimensionMismatchError(oss.str()); } } } // namespace +namespace linalgebra { + Vector::Vector(std::size_t n) : data_(n) {} Vector::Vector(std::size_t n, double value) : data_(n, value) {} @@ -101,4 +140,4 @@ double dot(const Vector& lhs, const Vector& rhs) { return std::inner_product(lhs.begin(), lhs.end(), rhs.begin(), 0.0); } -} // namespace linalg +} // namespace linalgebra diff --git a/tests/test_lu.cpp b/tests/test_lu.cpp index dfb76c4..c4f9c5d 100644 --- a/tests/test_lu.cpp +++ b/tests/test_lu.cpp @@ -1,8 +1,4 @@ -#include "linalg_error.hpp" -#include "lu.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "vector.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> @@ -11,31 +7,23 @@ #include <cstddef> #include <random> -using linalg::DimensionMismatchError; -using linalg::Matrix; -using linalg::SingularMatrixError; -using linalg::Vector; -using linalg::LUResult; - -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::SingularMatrixError; +using linalgebra::Vector; +using linalgebra::LUResult; namespace { -// ||PA - LU||_F (Frobenius, computed element-wise via ||vec||_2) double reconstruction_error(const Matrix& A, const LUResult& lu) { const std::size_t n = A.rows(); - // Build PA by permuting rows of A. Matrix PA(n, n); for (std::size_t i = 0; i < n; ++i) { for (std::size_t j = 0; j < n; ++j) { PA(i, j) = A(lu.perm[i], j); } } - // Compute LU product. const Matrix LU = lu.L * lu.U; - // Compute Frobenius norm of (PA - LU). double err = 0.0; for (std::size_t i = 0; i < n; ++i) { for (std::size_t j = 0; j < n; ++j) { @@ -46,12 +34,10 @@ double reconstruction_error(const Matrix& A, const LUResult& lu) { return std::sqrt(err); } -// ||Ax - b||_2 double solve_residual(const Matrix& A, const Vector& x, const Vector& b) { - return linalg::norm2(A * x - b); + return linalgebra::norm2(A * x - b); } -// Generate a reproducible random nonsingular n x n matrix. Matrix random_matrix(std::size_t n, unsigned seed = 42) { std::mt19937 rng(seed); std::uniform_real_distribution<double> dist(-10.0, 10.0); @@ -66,17 +52,13 @@ Matrix random_matrix(std::size_t n, unsigned seed = 42) { } // namespace -// --------------------------------------------------------------------------- -// Factorization correctness -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: 3x3 known system", "[lu]") { const Matrix A{ {2.0, 1.0, -1.0}, {-3.0, -1.0, 2.0}, {-2.0, 1.0, 2.0} }; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); REQUIRE(lu.L.rows() == 3); REQUIRE(lu.U.rows() == 3); @@ -86,41 +68,33 @@ TEST_CASE("LU factorization: 3x3 known system", "[lu]") { CHECK(lu.L(i, i) == Catch::Approx(1.0)); } - // ||PA - LU|| must be near zero. CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-12)); } TEST_CASE("LU factorization: identity matrix", "[lu]") { const Matrix I = Matrix::identity(4); - const LUResult lu = linalg::lu_factor(I); + const LUResult lu = linalgebra::lu_factor(I); CHECK(reconstruction_error(I, lu) == Catch::Approx(0.0).margin(1e-14)); - // U should equal I (up to row ordering already handled by PA=LU). for (std::size_t i = 0; i < 4; ++i) { CHECK(lu.U(i, i) == Catch::Approx(1.0)); } } TEST_CASE("LU factorization: matrix requiring row swaps", "[lu]") { - // First column entry is zero. No-pivot LU would immediately fail. const Matrix A{ {0.0, 1.0, 2.0}, {3.0, 4.0, 5.0}, {6.0, 7.0, 8.0} }; - // A is singular (rows are in AP), but check that partial pivoting still - // proceeds and detects singularity correctly. - // Row 3 - row 2 = row 2 - row 1, so rank < 3. - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: first-column zero, nonsingular", "[lu]") { - // [[0, 1], [1, 0]] — requires a swap at step 0. const Matrix A{{0.0, 1.0}, {1.0, 0.0}}; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-14)); - // Solving Ax = b: A swaps components. const Vector b{3.0, 7.0}; - const Vector x = linalg::lu_solve(lu, b); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-12)); CHECK(x[0] == Catch::Approx(7.0)); CHECK(x[1] == Catch::Approx(3.0)); @@ -129,17 +103,12 @@ TEST_CASE("LU factorization: first-column zero, nonsingular", "[lu]") { TEST_CASE("LU factorization: random nonsingular matrices", "[lu]") { for (std::size_t n : {5u, 10u, 20u}) { const Matrix A = random_matrix(n, 123u + static_cast<unsigned>(n)); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-10)); } } -// --------------------------------------------------------------------------- -// Solve correctness -// --------------------------------------------------------------------------- - TEST_CASE("LU solve: known 3x3 system", "[lu]") { - // From Cramer / textbook: solution is x = (2, 3, -1). const Matrix A{ {2.0, 1.0, -1.0}, {-3.0, -1.0, 2.0}, @@ -148,8 +117,8 @@ TEST_CASE("LU solve: known 3x3 system", "[lu]") { const Vector expected{2.0, 3.0, -1.0}; const Vector b = A * expected; - const LUResult lu = linalg::lu_factor(A); - const Vector x = linalg::lu_solve(lu, b); + const LUResult lu = linalgebra::lu_factor(A); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(x[0] == Catch::Approx(expected[0]).epsilon(1e-12)); CHECK(x[1] == Catch::Approx(expected[1]).epsilon(1e-12)); @@ -160,92 +129,77 @@ TEST_CASE("LU solve: known 3x3 system", "[lu]") { TEST_CASE("LU solve: random nonsingular systems", "[lu]") { for (std::size_t n : {5u, 15u, 30u}) { const Matrix A = random_matrix(n, 7u * static_cast<unsigned>(n)); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); - // Random rhs. std::mt19937 rng(n); std::uniform_real_distribution<double> dist(-5.0, 5.0); Vector b(n); for (std::size_t i = 0; i < n; ++i) b[i] = dist(rng); - const Vector x = linalg::lu_solve(lu, b); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-9)); } } TEST_CASE("LU solve: diagonal system", "[lu]") { - // D = diag(2, 3, 4), b = (2, 9, 8), solution = (1, 3, 2). const Matrix D{ {2.0, 0.0, 0.0}, {0.0, 3.0, 0.0}, {0.0, 0.0, 4.0} }; const Vector b{2.0, 9.0, 8.0}; - const LUResult lu = linalg::lu_factor(D); - const Vector x = linalg::lu_solve(lu, b); + const LUResult lu = linalgebra::lu_factor(D); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(x[0] == Catch::Approx(1.0)); CHECK(x[1] == Catch::Approx(3.0)); CHECK(x[2] == Catch::Approx(2.0)); } -// --------------------------------------------------------------------------- -// Failure cases -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: non-square matrix throws", "[lu]") { const Matrix A(3, 4); - CHECK_THROWS_AS(linalg::lu_factor(A), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), DimensionMismatchError); } TEST_CASE("LU factorization: exactly singular matrix throws", "[lu]") { - // Zero row → singular. const Matrix A{ {1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {0.0, 0.0, 0.0} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: rank-deficient matrix throws", "[lu]") { - // Row 2 is a linear combination of rows 0 and 1. const Matrix A{ {1.0, 2.0}, {2.0, 4.0} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: near-singular matrix throws at default tolerance", "[lu]") { - // Pivot reduced to ~1e-16, should trip the singularity check. const Matrix A{ {1.0, 1.0}, {1.0, 1.0 + 1e-16} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU solve: mismatched rhs throws", "[lu]") { const Matrix A = Matrix::identity(3); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const Vector b(5, 1.0); - CHECK_THROWS_AS(linalg::lu_solve(lu, b), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::lu_solve(lu, b), DimensionMismatchError); } -// --------------------------------------------------------------------------- -// L and U structure -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: L is unit lower triangular", "[lu]") { const Matrix A = random_matrix(6, 999u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const std::size_t n = A.rows(); for (std::size_t i = 0; i < n; ++i) { - // Unit diagonal. CHECK(lu.L(i, i) == Catch::Approx(1.0)); - // Strict upper triangle is zero. for (std::size_t j = i + 1; j < n; ++j) { CHECK(lu.L(i, j) == Catch::Approx(0.0).margin(1e-15)); } @@ -254,7 +208,7 @@ TEST_CASE("LU factorization: L is unit lower triangular", "[lu]") { TEST_CASE("LU factorization: U is upper triangular", "[lu]") { const Matrix A = random_matrix(6, 777u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const std::size_t n = A.rows(); for (std::size_t i = 1; i < n; ++i) { @@ -264,20 +218,15 @@ TEST_CASE("LU factorization: U is upper triangular", "[lu]") { } } -// --------------------------------------------------------------------------- -// Permutation sign and determinant -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: sign of permutation is ±1", "[lu]") { const Matrix A = random_matrix(5, 321u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK((lu.sign == 1 || lu.sign == -1)); } TEST_CASE("LU factorization: determinant via sign * prod(diag(U))", "[lu]") { - // det([[3,1],[2,4]]) = 12 - 2 = 10 const Matrix A{{3.0, 1.0}, {2.0, 4.0}}; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); double det = static_cast<double>(lu.sign); for (std::size_t i = 0; i < A.rows(); ++i) { det *= lu.U(i, i); diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp index fcd3d27..2fd9a5e 100644 --- a/tests/test_matrix.cpp +++ b/tests/test_matrix.cpp @@ -1,15 +1,11 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> -#include <type_traits> -#include <utility> - -using linalg::Matrix; -using linalg::Vector; -using linalg::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::Vector; +using linalgebra::DimensionMismatchError; TEST_CASE("Matrix constructors initialize dimensions and values", "[matrix]") { const Matrix empty; @@ -101,7 +97,7 @@ TEST_CASE("Matrix initializer list rejects unequal row lengths", "[matrix]") { TEST_CASE("Matrix transpose swaps rows and columns", "[matrix]") { const Matrix a{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}}; - const Matrix at = linalg::transpose(a); + const Matrix at = linalgebra::transpose(a); REQUIRE(at.rows() == 3); REQUIRE(at.cols() == 2); diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp index a837132..1a4114f 100644 --- a/tests/test_qr.cpp +++ b/tests/test_qr.cpp @@ -1,8 +1,4 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "qr.hpp" -#include "vector.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> @@ -12,18 +8,13 @@ #include <functional> #include <random> -using linalg::DimensionMismatchError; -using linalg::Matrix; -using linalg::QRResult; -using linalg::SingularMatrixError; - -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::QRResult; +using linalgebra::SingularMatrixError; namespace { -// ||A - QR||_F double reconstruction_error(const Matrix& A, const QRResult& qr) { const Matrix diff = A - qr.Q * qr.R; double err = 0.0; @@ -33,11 +24,9 @@ double reconstruction_error(const Matrix& A, const QRResult& qr) { return std::sqrt(err); } -// ||Q^T Q - I||_F (should be ~0 for orthonormal Q) double orthogonality_error(const QRResult& qr) { const Matrix& Q = qr.Q; const std::size_t n = Q.cols(); - // Compute Q^T Q Matrix QtQ(n, n); for (std::size_t i = 0; i < n; ++i) for (std::size_t j = 0; j < n; ++j) { @@ -45,7 +34,6 @@ double orthogonality_error(const QRResult& qr) { for (std::size_t k = 0; k < Q.rows(); ++k) s += Q(k, i) * Q(k, j); QtQ(i, j) = s; } - // ||QtQ - I||_F double err = 0.0; for (std::size_t i = 0; i < n; ++i) for (std::size_t j = 0; j < n; ++j) { @@ -55,7 +43,6 @@ double orthogonality_error(const QRResult& qr) { return std::sqrt(err); } -// R must be upper triangular (strict lower triangle near zero). bool r_is_upper_triangular(const Matrix& R, double tol = 1e-12) { for (std::size_t i = 1; i < R.rows(); ++i) for (std::size_t j = 0; j < i; ++j) @@ -73,7 +60,6 @@ Matrix random_matrix(std::size_t m, std::size_t n, unsigned seed = 42) { return M; } -// Run all checks for a given QR function and matrix. using QRFn = std::function<QRResult(const Matrix&)>; void check_qr(const Matrix& A, QRFn fn, @@ -91,33 +77,25 @@ void check_qr(const Matrix& A, QRFn fn, } // namespace -// --------------------------------------------------------------------------- -// Macro to run the same test body for all three methods -// --------------------------------------------------------------------------- - -#define FOR_ALL_METHODS(A, recon_tol, ortho_tol) \ - SECTION("classical_gs") { \ - check_qr(A, [](const Matrix& M) { return linalg::qr_classical_gs(M); }, \ - recon_tol, ortho_tol, "classical_gs"); \ - } \ - SECTION("modified_gs") { \ - check_qr(A, [](const Matrix& M) { return linalg::qr_modified_gs(M); }, \ - recon_tol, ortho_tol, "modified_gs"); \ - } \ - SECTION("householder") { \ - check_qr(A, [](const Matrix& M) { return linalg::qr_householder(M); }, \ - recon_tol, ortho_tol, "householder"); \ +#define FOR_ALL_METHODS(A, recon_tol, ortho_tol) \ + SECTION("classical_gs") { \ + check_qr(A, [](const Matrix& M) { return linalgebra::qr_classical_gs(M); }, \ + recon_tol, ortho_tol, "classical_gs"); \ + } \ + SECTION("modified_gs") { \ + check_qr(A, [](const Matrix& M) { return linalgebra::qr_modified_gs(M); }, \ + recon_tol, ortho_tol, "modified_gs"); \ + } \ + SECTION("householder") { \ + check_qr(A, [](const Matrix& M) { return linalgebra::qr_householder(M); }, \ + recon_tol, ortho_tol, "householder"); \ } -// --------------------------------------------------------------------------- -// Basic correctness: square matrices -// --------------------------------------------------------------------------- - TEST_CASE("QR: 3x3 known matrix", "[qr]") { const Matrix A{ {1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, - {7.0, 8.0, 10.0} // not exactly singular + {7.0, 8.0, 10.0} }; FOR_ALL_METHODS(A, 1e-12, 1e-12) } @@ -136,10 +114,6 @@ TEST_CASE("QR: diagonal matrix", "[qr]") { FOR_ALL_METHODS(D, 1e-14, 1e-14) } -// --------------------------------------------------------------------------- -// Rectangular (tall) matrices -// --------------------------------------------------------------------------- - TEST_CASE("QR: tall 5x3 random matrix", "[qr]") { const Matrix A = random_matrix(5, 3, 7u); FOR_ALL_METHODS(A, 1e-12, 1e-12) @@ -150,10 +124,6 @@ TEST_CASE("QR: tall 10x4 random matrix", "[qr]") { FOR_ALL_METHODS(A, 1e-12, 1e-12) } -// --------------------------------------------------------------------------- -// Random square matrices -// --------------------------------------------------------------------------- - TEST_CASE("QR: random 6x6", "[qr]") { const Matrix A = random_matrix(6, 6, 123u); FOR_ALL_METHODS(A, 1e-12, 1e-12) @@ -164,14 +134,7 @@ TEST_CASE("QR: random 12x12", "[qr]") { FOR_ALL_METHODS(A, 1e-11, 1e-11) } -// --------------------------------------------------------------------------- -// Nearly dependent columns — GS methods degrade; Householder stays clean -// --------------------------------------------------------------------------- - TEST_CASE("QR: nearly dependent columns", "[qr]") { - // Column 1 = column 0 + epsilon * e_1. - // Classical GS will lose most of Q's orthogonality here. - // Modified GS is better. Householder is unaffected. constexpr double eps = 1e-7; const Matrix A{ {1.0, 1.0 + eps, 0.0}, @@ -180,32 +143,24 @@ TEST_CASE("QR: nearly dependent columns", "[qr]") { {0.0, 0.0, 1.0} }; - // All three should reconstruct A accurately. SECTION("classical_gs reconstruction") { - const QRResult qr = linalg::qr_classical_gs(A); + const QRResult qr = linalgebra::qr_classical_gs(A); CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10)); - // Orthogonality will be poor for classical GS on this input. - // We only assert it's not catastrophically wrong (< 0.01). CHECK(orthogonality_error(qr) < 0.01); } SECTION("modified_gs reconstruction") { - const QRResult qr = linalg::qr_modified_gs(A); + const QRResult qr = linalgebra::qr_modified_gs(A); CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10)); CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-8)); } SECTION("householder reconstruction") { - const QRResult qr = linalg::qr_householder(A); + const QRResult qr = linalgebra::qr_householder(A); CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-13)); CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13)); } } -// --------------------------------------------------------------------------- -// Hilbert-like ill-conditioned matrix -// --------------------------------------------------------------------------- - TEST_CASE("QR: 4x4 Hilbert matrix", "[qr]") { - // H[i][j] = 1 / (i + j + 1) const std::size_t n = 4; Matrix H(n, n); for (std::size_t i = 0; i < n; ++i) @@ -213,39 +168,34 @@ TEST_CASE("QR: 4x4 Hilbert matrix", "[qr]") { H(i, j) = 1.0 / static_cast<double>(i + j + 1); SECTION("classical_gs") { - const QRResult qr = linalg::qr_classical_gs(H); + const QRResult qr = linalgebra::qr_classical_gs(H); CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12)); - // Orthogonality is imperfect on Hilbert matrices with classical GS. CHECK(orthogonality_error(qr) < 1e-8); } SECTION("modified_gs") { - const QRResult qr = linalg::qr_modified_gs(H); + const QRResult qr = linalgebra::qr_modified_gs(H); CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12)); CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-10)); } SECTION("householder") { - const QRResult qr = linalg::qr_householder(H); + const QRResult qr = linalgebra::qr_householder(H); CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-13)); CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13)); } } -// --------------------------------------------------------------------------- -// Structure checks -// --------------------------------------------------------------------------- - TEST_CASE("QR: R is upper triangular", "[qr]") { const Matrix A = random_matrix(5, 5, 555u); - CHECK(r_is_upper_triangular(linalg::qr_classical_gs(A).R)); - CHECK(r_is_upper_triangular(linalg::qr_modified_gs(A).R)); - CHECK(r_is_upper_triangular(linalg::qr_householder(A).R)); + CHECK(r_is_upper_triangular(linalgebra::qr_classical_gs(A).R)); + CHECK(r_is_upper_triangular(linalgebra::qr_modified_gs(A).R)); + CHECK(r_is_upper_triangular(linalgebra::qr_householder(A).R)); } TEST_CASE("QR: Q columns are unit length", "[qr]") { const Matrix A = random_matrix(6, 4, 321u); - for (QRFn fn : {QRFn{[](const Matrix& M) { return linalg::qr_classical_gs(M); }}, - QRFn{[](const Matrix& M) { return linalg::qr_modified_gs(M); }}, - QRFn{[](const Matrix& M) { return linalg::qr_householder(M); }}}) { + for (QRFn fn : {QRFn{[](const Matrix& M) { return linalgebra::qr_classical_gs(M); }}, + QRFn{[](const Matrix& M) { return linalgebra::qr_modified_gs(M); }}, + QRFn{[](const Matrix& M) { return linalgebra::qr_householder(M); }}}) { const QRResult qr = fn(A); for (std::size_t j = 0; j < qr.Q.cols(); ++j) { double norm2 = 0.0; @@ -256,25 +206,20 @@ TEST_CASE("QR: Q columns are unit length", "[qr]") { } } -// --------------------------------------------------------------------------- -// Failure cases -// --------------------------------------------------------------------------- - TEST_CASE("QR: fat matrix throws DimensionMismatchError", "[qr]") { - const Matrix A(3, 5); // rows < cols - CHECK_THROWS_AS(linalg::qr_classical_gs(A), DimensionMismatchError); - CHECK_THROWS_AS(linalg::qr_modified_gs(A), DimensionMismatchError); - CHECK_THROWS_AS(linalg::qr_householder(A), DimensionMismatchError); + const Matrix A(3, 5); + CHECK_THROWS_AS(linalgebra::qr_classical_gs(A), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::qr_modified_gs(A), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::qr_householder(A), DimensionMismatchError); } TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") { const Matrix A{ - {1.0, 2.0, 2.0}, // col 2 = 2 * col 0 + {1.0, 2.0, 2.0}, {2.0, 4.0, 4.0}, {3.0, 6.0, 6.0} }; - CHECK_THROWS_AS(linalg::qr_classical_gs(A), SingularMatrixError); - CHECK_THROWS_AS(linalg::qr_modified_gs(A), SingularMatrixError); - // Householder handles rank-deficient input gracefully (R gets a zero diagonal entry). - CHECK_NOTHROW(linalg::qr_householder(A)); + CHECK_THROWS_AS(linalgebra::qr_classical_gs(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::qr_modified_gs(A), SingularMatrixError); + CHECK_NOTHROW(linalgebra::qr_householder(A)); } diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index 7126334..008f109 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -1,7 +1,4 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "qr_iteration.hpp" -#include "vector.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> @@ -16,19 +13,14 @@ #include <utility> #include <vector> -using linalg::Matrix; -using linalg::NonConvergenceError; -using linalg::QRIterationOptions; -using linalg::QRIterationResult; -using linalg::Vector; - -// --------------------------------------------------------------------------- -// Test helpers -// --------------------------------------------------------------------------- +using linalgebra::Matrix; +using linalgebra::NonConvergenceError; +using linalgebra::QRIterationOptions; +using linalgebra::QRIterationResult; +using linalgebra::Vector; namespace { -// Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag. using EigPairs = std::vector<std::pair<double, double>>; EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { @@ -66,20 +58,6 @@ bool eigs_match(const Vector& computed_real, const Vector& computed_imag, } // namespace -// --------------------------------------------------------------------------- -// Test 1: 2×2 symmetric matrix with known eigenvalues -// --------------------------------------------------------------------------- -// -// A = | 2 1 | is symmetric positive definite. -// | 1 2 | -// -// Characteristic polynomial: (2-λ)^2 - 1 = 0 → λ = 1, 3. -// Eigenvectors: [1,-1]/√2 (λ=1) and [1,1]/√2 (λ=3). -// -// The unshifted iteration converges at rate |λ_1/λ_2| = 1/3 per step, -// so only a handful of iterations are needed. -// Ref: T&B Theorem 28.2. - TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", "[qr_iteration][shifted]") { const Matrix A{ @@ -87,13 +65,12 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", {1.0, 2.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(A); REQUIRE(res.eigenvalues_real.size() == 2); REQUIRE(res.eigenvalues_imag.size() == 2); REQUIRE(res.iterations > 0); - // All eigenvalues of a symmetric matrix must be real. CHECK(std::abs(res.eigenvalues_imag[0]) < 1e-8); CHECK(std::abs(res.eigenvalues_imag[1]) < 1e-8); @@ -101,23 +78,6 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } -// --------------------------------------------------------------------------- -// Test 2: 4×4 symmetric tridiagonal matrix — reference eigenvalues -// --------------------------------------------------------------------------- -// -// The n×n symmetric tridiagonal matrix with 2 on the diagonal and -1 on the -// first super- and sub-diagonals has known eigenvalues (discrete Laplacian): -// -// λ_k = 2 - 2 cos(k π / (n+1)), k = 1, …, n -// -// Ref: Golub & Van Loan §4.4.2 (discrete sine transform). -// -// For n = 4: -// λ_1 = 2 - 2 cos(π/5) ≈ 0.3820 -// λ_2 = 2 - 2 cos(2π/5) ≈ 1.3820 -// λ_3 = 2 - 2 cos(3π/5) ≈ 2.6180 -// λ_4 = 2 - 2 cos(4π/5) ≈ 3.6180 - TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", "[qr_iteration][unshifted]") { const Matrix A{ @@ -127,16 +87,14 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", { 0.0, 0.0, -1.0, 2.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(A); REQUIRE(res.eigenvalues_real.size() == 4); REQUIRE(res.eigenvalues_imag.size() == 4); - // All eigenvalues of a symmetric matrix must be real. for (std::size_t k = 0; k < 4; ++k) CHECK(std::abs(res.eigenvalues_imag[k]) < 1e-8); - // Compare against the closed-form reference. constexpr double pi = 3.14159265358979323846; const EigPairs expected = { {2.0 - 2.0 * std::cos( pi / 5.0), 0.0}, @@ -147,21 +105,6 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } -// --------------------------------------------------------------------------- -// Test 3: 5×5 symmetric tridiagonal — convergence history -// --------------------------------------------------------------------------- -// -// Uses a 5×5 symmetric tridiagonal (discrete Laplacian) to guarantee all -// real eigenvalues and predictable linear convergence. The Frobenius norm -// of the strict lower triangle is printed at every step so the convergence -// rate can be observed directly. -// -// Expected behaviour: ||lower(A_k)||_F decreases geometrically each step -// (linear convergence), with ratio ≈ max_j |λ_{j+1}/λ_j|. -// Ref: T&B Theorem 28.2. -// -// 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5. - TEST_CASE("QR iteration (unshifted): 5x5 convergence history", "[qr_iteration][unshifted]") { const Matrix A{ @@ -175,12 +118,11 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", QRIterationOptions opts; opts.track_convergence = true; - const QRIterationResult res = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(A, opts); REQUIRE_FALSE(res.convergence_history.empty()); REQUIRE(res.eigenvalues_real.size() == 5); - // Print convergence history so the linear rate is visible. std::cout << "\n=== Unshifted QR — 5x5 convergence history ===\n"; std::cout << " Converged in " << res.iterations << " iteration(s)\n"; for (std::size_t k = 0; k < res.convergence_history.size(); ++k) { @@ -190,32 +132,26 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", } std::cout << "=======================================================\n"; - // The final recorded norm must be below the default tolerance. CHECK(res.convergence_history.back() < opts.tolerance); } -// --------------------------------------------------------------------------- -// Test 4: Eigenvalue residuals below 1e-8 -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (unshifted): residuals below 1e-8", "[qr_iteration][unshifted]") { SECTION("2x2: eigenvalues 1 and 3") { const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(A); const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } SECTION("3x3 diagonal: eigenvalues 1, 4, 9") { - // Diagonal matrix — already in Schur form; converges in one step. const Matrix D{ {1.0, 0.0, 0.0}, {0.0, 4.0, 0.0}, {0.0, 0.0, 9.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(D); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(D); const EigPairs expected = {{1.0, 0.0}, {4.0, 0.0}, {9.0, 0.0}}; CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } @@ -234,13 +170,13 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(A); CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } SECTION("5x5 identity: all eigenvalues == 1") { const Matrix I = Matrix::identity(5); - const QRIterationResult res = linalg::eigenvalues_unshifted(I); + const QRIterationResult res = linalgebra::eigenvalues_unshifted(I); REQUIRE(res.eigenvalues_real.size() == 5); for (std::size_t k = 0; k < 5; ++k) { CHECK(std::abs(res.eigenvalues_real[k] - 1.0) < 1e-8); @@ -249,35 +185,21 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", } } -// --------------------------------------------------------------------------- -// Failure cases -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (unshifted): non-square matrix throws", "[qr_iteration][unshifted]") { - const Matrix A(3, 4); // non-square - CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A), - linalg::DimensionMismatchError); + const Matrix A(3, 4); + CHECK_THROWS_AS(linalgebra::eigenvalues_unshifted(A), + linalgebra::DimensionMismatchError); } TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws", "[qr_iteration][unshifted]") { - // Cap at zero iterations — any non-trivial matrix fails immediately. const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; QRIterationOptions opts; opts.max_iterations = 0; - CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A, opts), NonConvergenceError); + CHECK_THROWS_AS(linalgebra::eigenvalues_unshifted(A, opts), NonConvergenceError); } -// =========================================================================== -// Wilkinson-shifted QR iteration -// =========================================================================== - -// --------------------------------------------------------------------------- -// Helper builds a random symmetric matrix via A = M + M^T (guaranteed real -// eigenvalues) with a fixed seed for reproducibility. -// --------------------------------------------------------------------------- - namespace { Matrix random_symmetric(std::size_t n, unsigned seed = 42) { @@ -296,14 +218,6 @@ Matrix random_symmetric(std::size_t n, unsigned seed = 42) { } // namespace -// --------------------------------------------------------------------------- -// Test S1: shifted vs unshifted iteration count on the same matrix. -// -// Wilkinson-shifted QR converges (typically cubically) in far fewer steps -// than the unshifted algorithm (linear convergence). -// The test asserts the shifted count is strictly smaller and prints both. -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", "[qr_iteration][shifted]") { const Matrix A{ @@ -318,8 +232,8 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", QRIterationOptions opts; opts.track_convergence = true; - const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); - const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts); std::cout << "\n=== Shifted vs Unshifted ===\n"; std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; @@ -333,25 +247,15 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", 1e-8)); } -// --------------------------------------------------------------------------- -// Test S2: matrix where unshifted takes >100 iterations, shifted takes <20. -// -// A nearly-equal-eigenvalue symmetric matrix maximises the linear convergence -// slowdown. Using a scaled identity perturbation: eigenvalues cluster near 1, -// slowing unshifted (ratio ≈ 1) while the Wilkinson shift adapts instantly. -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >100", "[qr_iteration][shifted]") { - // 5×5 symmetric matrix with eigenvalues 1, 1.001, 1.002, 1.003, 1.004. - // Off-diagonal entries couple them. Unshifted stalls (|λ_{j+1}/λ_j| ≈ 1). const Matrix A = random_symmetric(5, 17u); QRIterationOptions opts; - opts.max_iterations = 2000; + opts.max_iterations = 2000; - const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); - const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts); std::cout << "\n=== Hard matrix ===\n"; std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; @@ -361,10 +265,6 @@ TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >10 CHECK(shifted.iterations < 20); } -// --------------------------------------------------------------------------- -// Test S3: shifted eigenvalues match known values to within 1e-8. -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (shifted): residuals below 1e-8", "[qr_iteration][shifted]") { SECTION("4x4 tridiagonal: closed-form eigenvalues") { @@ -381,29 +281,18 @@ TEST_CASE("QR iteration (shifted): residuals below 1e-8", {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} }; - const QRIterationResult res = linalg::eigenvalues_shifted(A); + const QRIterationResult res = linalgebra::eigenvalues_shifted(A); CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } SECTION("2x2 known eigenvalues") { const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; - const QRIterationResult res = linalg::eigenvalues_shifted(A); + const QRIterationResult res = linalgebra::eigenvalues_shifted(A); CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } } -// =========================================================================== -// Hessenberg reduction + practical QR algorithm -// =========================================================================== - -// --------------------------------------------------------------------------- -// Test H1: hessenberg_reduction produces correct H and Q. -// -// Verify: (1) H is upper Hessenberg, (2) Q is orthogonal, (3) A = Q H Q^T. -// Ref: GVL §7.4.2. -// --------------------------------------------------------------------------- - namespace { bool is_upper_hessenberg(const Matrix& H, double tol = 1e-10) { @@ -421,7 +310,6 @@ double frobenius_norm(const Matrix& A) { return std::sqrt(s); } -// ||A - B||_F double diff_norm(const Matrix& A, const Matrix& B) { double s = 0.0; for (std::size_t i = 0; i < A.rows(); ++i) @@ -432,7 +320,6 @@ double diff_norm(const Matrix& A, const Matrix& B) { return std::sqrt(s); } -// ||Q^T Q - I||_F double orthogonality_error(const Matrix& Q) { const std::size_t n = Q.rows(); double err = 0.0; @@ -451,29 +338,21 @@ double orthogonality_error(const Matrix& Q) { TEST_CASE("Hessenberg reduction: structure and similarity", "[qr_iteration][shifted]") { const Matrix A = random_symmetric(6, 7u); - const linalg::HessenbergResult hr = linalg::hessenberg_reduction(A); + const linalgebra::HessenbergResult hr = linalgebra::hessenberg_reduction(A); - // H must be upper Hessenberg. CHECK(is_upper_hessenberg(hr.H)); - - // Q must be orthogonal. CHECK(orthogonality_error(hr.Q) < 1e-10); - // A = Q H Q^T ⟹ ||A - Q H Q^T||_F < tol. - const Matrix QtHQ = hr.Q * hr.H * linalg::transpose(hr.Q); + const Matrix QtHQ = hr.Q * hr.H * linalgebra::transpose(hr.Q); CHECK(diff_norm(A, QtHQ) < 1e-10); } -// --------------------------------------------------------------------------- -// Test H2: eigenvalues_hessenberg agrees with eigenvalues_shifted to 1e-6. -// --------------------------------------------------------------------------- - TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6", "[qr_iteration][shifted]") { const Matrix A = random_symmetric(8, 99u); - const QRIterationResult ref = linalg::eigenvalues_shifted(A); - const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + const QRIterationResult ref = linalgebra::eigenvalues_shifted(A); + const QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A); REQUIRE(hess.eigenvalues_real.size() == 8); CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, @@ -481,10 +360,6 @@ TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6", 1e-6)); } -// --------------------------------------------------------------------------- -// Test H3: known eigenvalues — 4×4 tridiagonal. -// --------------------------------------------------------------------------- - TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix", "[qr_iteration][shifted]") { const Matrix A{ @@ -500,19 +375,10 @@ TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix", {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} }; - const QRIterationResult res = linalg::eigenvalues_hessenberg(A); + const QRIterationResult res = linalgebra::eigenvalues_hessenberg(A); CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); } -// --------------------------------------------------------------------------- -// Test H4: benchmark — shifted QR vs Hessenberg pipeline for n = 50, 100, 200. -// -// The Hessenberg pipeline reduces each QR step from O(n³) to O(n²), so the -// speedup should grow with n. We print the wall-clock ratio and assert that -// the Hessenberg version is faster for n >= 50. -// Ref: GVL §7.4.2; T&B Lecture 29. -// --------------------------------------------------------------------------- - TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", "[qr_iteration][hessenberg]") { using Clock = std::chrono::high_resolution_clock; @@ -531,11 +397,11 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", const Matrix A = random_symmetric(n, 13u); const auto t0s = Clock::now(); - { const auto tmp = linalg::eigenvalues_shifted(A); (void)tmp; } + { const auto tmp = linalgebra::eigenvalues_shifted(A); (void)tmp; } const double t_shifted = Seconds(Clock::now() - t0s).count(); const auto t0h = Clock::now(); - const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + const QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A); const double t_hess = Seconds(Clock::now() - t0h).count(); const double speedup = t_shifted / t_hess; @@ -547,11 +413,9 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", << std::setprecision(2) << std::setw(12) << speedup << "x\n"; - // The Hessenberg version must be faster for all tested sizes. CHECK(t_hess < t_shifted); - // And must give correct eigenvalues (agree with shifted to 1e-6). - const QRIterationResult ref = linalg::eigenvalues_shifted(A); + const QRIterationResult ref = linalgebra::eigenvalues_shifted(A); CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), 1e-6)); diff --git a/tests/test_triangular_solve.cpp b/tests/test_triangular_solve.cpp index da7e55c..f1ce709 100644 --- a/tests/test_triangular_solve.cpp +++ b/tests/test_triangular_solve.cpp @@ -1,20 +1,17 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "triangular_solve.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> -using linalg::DimensionMismatchError; -using linalg::Matrix; -using linalg::SingularMatrixError; -using linalg::Vector; +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::SingularMatrixError; +using linalgebra::Vector; namespace { double residual_norm(const Matrix& a, const Vector& x, const Vector& b) { - return linalg::norm2((a * x) - b); + return linalgebra::norm2((a * x) - b); } } // namespace @@ -28,7 +25,7 @@ TEST_CASE("Forward substitution solves lower-triangular systems", "[triangular]" const Vector expected{1.0, 2.0, -1.0}; const Vector rhs = lower * expected; - const Vector x = linalg::forward_substitution(lower, rhs); + const Vector x = linalgebra::forward_substitution(lower, rhs); REQUIRE(x.size() == expected.size()); CHECK(x[0] == Catch::Approx(expected[0])); CHECK(x[1] == Catch::Approx(expected[1])); @@ -45,7 +42,7 @@ TEST_CASE("Backward substitution solves upper-triangular systems", "[triangular] const Vector expected{2.0, -1.0, 3.0}; const Vector rhs = upper * expected; - const Vector x = linalg::backward_substitution(upper, rhs); + const Vector x = linalgebra::backward_substitution(upper, rhs); REQUIRE(x.size() == expected.size()); CHECK(x[0] == Catch::Approx(expected[0])); CHECK(x[1] == Catch::Approx(expected[1])); @@ -61,7 +58,7 @@ TEST_CASE("Triangular solves support unit-diagonal systems", "[triangular]") { }; const Vector rhs{1.0, 0.0, 4.0}; - const Vector x = linalg::forward_substitution(lower, rhs, 1e-12, true); + const Vector x = linalgebra::forward_substitution(lower, rhs, 1e-12, true); CHECK(x[0] == Catch::Approx(1.0)); CHECK(x[1] == Catch::Approx(2.0)); CHECK(x[2] == Catch::Approx(3.0)); @@ -70,19 +67,19 @@ TEST_CASE("Triangular solves support unit-diagonal systems", "[triangular]") { TEST_CASE("Triangular solves reject shape and structure mismatches", "[triangular]") { const Matrix nonsquare(2, 3); const Vector rhs2{1.0, 2.0}; - CHECK_THROWS_AS(linalg::forward_substitution(nonsquare, rhs2), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::forward_substitution(nonsquare, rhs2), DimensionMismatchError); const Matrix lower{ {1.0, 1.0}, {2.0, 3.0} }; - CHECK_THROWS_AS(linalg::forward_substitution(lower, rhs2), std::invalid_argument); + CHECK_THROWS_AS(linalgebra::forward_substitution(lower, rhs2), std::invalid_argument); const Matrix upper{ {1.0, 2.0}, {1.0, 3.0} }; - CHECK_THROWS_AS(linalg::backward_substitution(upper, rhs2), std::invalid_argument); + CHECK_THROWS_AS(linalgebra::backward_substitution(upper, rhs2), std::invalid_argument); } TEST_CASE("Triangular solves detect negligible diagonal entries", "[triangular]") { @@ -91,11 +88,11 @@ TEST_CASE("Triangular solves detect negligible diagonal entries", "[triangular]" {2.0, 1.0} }; const Vector rhs{1.0, 2.0}; - CHECK_THROWS_AS(linalg::forward_substitution(lower, rhs), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::forward_substitution(lower, rhs), SingularMatrixError); const Matrix upper{ {1.0, 2.0}, {0.0, 1e-14} }; - CHECK_THROWS_AS(linalg::backward_substitution(upper, rhs), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::backward_substitution(upper, rhs), SingularMatrixError); } diff --git a/tests/test_vector.cpp b/tests/test_vector.cpp index 364ce31..a776dc8 100644 --- a/tests/test_vector.cpp +++ b/tests/test_vector.cpp @@ -1,6 +1,4 @@ -#include "linalg_error.hpp" -#include "norms.hpp" -#include "vector.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> @@ -8,8 +6,8 @@ #include <type_traits> #include <utility> -using linalg::DimensionMismatchError; -using linalg::Vector; +using linalgebra::DimensionMismatchError; +using linalgebra::Vector; TEST_CASE("Vector constructors initialize size and values", "[vector]") { const Vector empty; @@ -82,17 +80,17 @@ TEST_CASE("Vector arithmetic enforces shape compatibility", "[vector]") { CHECK(scaled[1] == 1.0); CHECK(scaled[2] == 1.5); - CHECK(linalg::dot(a, b) == 32.0); + CHECK(linalgebra::dot(a, b) == 32.0); const Vector short_vec{1.0, 2.0}; CHECK_THROWS_AS(a + short_vec, DimensionMismatchError); - CHECK_THROWS_AS(linalg::dot(a, short_vec), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::dot(a, short_vec), DimensionMismatchError); } TEST_CASE("Vector 2-norm matches manually computed values", "[vector][norms]") { const Vector v{3.0, 4.0}; - CHECK(linalg::norm2(v) == Catch::Approx(5.0)); + CHECK(linalgebra::norm2(v) == Catch::Approx(5.0)); const Vector zero(5); - CHECK(linalg::norm2(zero) == Catch::Approx(0.0)); + CHECK(linalgebra::norm2(zero) == Catch::Approx(0.0)); } |