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/qr.cpp | 87 ++++++++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 56 insertions(+), 31 deletions(-) (limited to 'src/qr.cpp') 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 -- cgit v1.2.3