From 4602b36e9d5ea08656e3222846a1f161bbb1cec1 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Mon, 27 Apr 2026 07:24:41 +0300 Subject: Module refactor --- src/linalgebra.cpp | 10 ++ src/linalgebra_error.cpp | 29 +++++ src/lu.cpp | 34 +++--- src/matrix.cpp | 162 ++++++++++---------------- src/norms.cpp | 14 ++- src/qr.cpp | 87 +++++++++----- src/qr_iteration.cpp | 287 ++++++++++++++++------------------------------- src/triangular_solve.cpp | 65 ++++++----- src/vector.cpp | 59 ++++++++-- 9 files changed, 365 insertions(+), 382 deletions(-) create mode 100644 src/linalgebra.cpp create mode 100644 src/linalgebra_error.cpp (limited to 'src') 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/src/linalgebra_error.cpp b/src/linalgebra_error.cpp new file mode 100644 index 0000000..367c79d --- /dev/null +++ b/src/linalgebra_error.cpp @@ -0,0 +1,29 @@ +export module linalgebra:error; +import std; + +export namespace linalgebra { + +class LinAlgError : public std::runtime_error { +public: + using std::runtime_error::runtime_error; +}; + +class DimensionMismatchError : public LinAlgError { +public: + explicit DimensionMismatchError(const std::string& message) + : LinAlgError(message) {} +}; + +class SingularMatrixError : public LinAlgError { +public: + explicit SingularMatrixError(const std::string& message) + : LinAlgError(message) {} +}; + +class NonConvergenceError : public LinAlgError { +public: + explicit NonConvergenceError(const std::string& message) + : LinAlgError(message) {} +}; + +} // namespace linalgebra diff --git a/src/lu.cpp b/src/lu.cpp index f6840ea..6d12eda 100644 --- a/src/lu.cpp +++ b/src/lu.cpp @@ -1,14 +1,26 @@ -#include "lu.hpp" +export module linalgebra:lu; +import std; +import :error; +import :vector; +import :matrix; +import :triangular_solve; -#include -#include -#include -#include +export namespace linalgebra { -#include "linalg_error.hpp" -#include "triangular_solve.hpp" +struct LUResult { + Matrix L; + Matrix U; + std::vector 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 -#include -#include -#include - -#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \ - (defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)) -#include -#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 #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> 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 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> values) : rows_(values.size()) { +Matrix::Matrix(std::initializer_list> values) + : rows_(values.size()) { if (rows_ == 0) { cols_ = 0; return; @@ -178,7 +131,6 @@ Matrix::Matrix(std::initializer_list> 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 +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 diff --git a/src/qr.cpp b/src/qr.cpp index 06770b2..7e32ce3 100644 --- a/src/qr.cpp +++ b/src/qr.cpp @@ -1,23 +1,55 @@ -#include "qr.hpp" - -#include -#include - -#include "linalg_error.hpp" - -namespace linalg { +export module linalgebra:qr; +import std; +import :error; +import :vector; +import :matrix; + +export namespace linalgebra { + +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 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); // + 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); // + 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 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 -#include -#include -#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 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 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= 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 -#include -#include +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 -#include -#include -#include +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 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 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 -- cgit v1.2.3