From 750f276a1403c5defd58b06f13c703b5e3d59245 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sat, 16 May 2026 12:10:05 +0300 Subject: Finish TODOs --- src/precond.cpp | 279 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100644 src/precond.cpp (limited to 'src/precond.cpp') diff --git a/src/precond.cpp b/src/precond.cpp new file mode 100644 index 0000000..a4ed503 --- /dev/null +++ b/src/precond.cpp @@ -0,0 +1,279 @@ +export module linalgebra:precond; +import std; +import :error; +import :vector; +import :matrix; +import :norms; +import :triangular_solve; +import :lu; +import :qr; + +// References used throughout this file: +// T&B — Trefethen & Bau, "Numerical Linear Algebra" +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. +// Hager — Hager, W.W. (1984), SIAM J. Sci. Stat. Comput. 5(2):311-316 + +export namespace linalgebra { + +// Condition number estimation (1-norm, Hager/LINPACK power iteration) +// Reference: Hager (1984); GVL §2.3.3 + +[[nodiscard]] double condition_number_1norm(const Matrix& A, + double singular_tolerance = 1e-12); + +// Jacobi (diagonal) preconditioner + +struct JacobiPrecond { + Vector inv_diag; +}; + +[[nodiscard]] JacobiPrecond precond_jacobi(const Matrix& A, + double zero_tolerance = 1e-14); + +[[nodiscard]] Vector apply(const JacobiPrecond& P, const Vector& x); + +// ILU(0) preconditioner (for dense matrices = LU without pivoting) +// Reference: Saad, "Iterative Methods for Sparse Linear Systems" §10.3 + +struct ILU0Precond { + Matrix LU; // combined: strict lower = L multipliers, upper = U +}; + +[[nodiscard]] ILU0Precond precond_ilu0(const Matrix& A, + double zero_tolerance = 1e-14); + +[[nodiscard]] Vector apply(const ILU0Precond& P, const Vector& b); + +// Least-squares solver via column-pivoting QR +// Reference: T&B Lecture 11; GVL §5.5 + +struct LstsqResult { + Vector x; + std::size_t rank; + double residual_norm; +}; + +struct LstsqOptions { + double rank_tolerance = 1e-12; +}; + +[[nodiscard]] LstsqResult lstsq(const Matrix& A, const Vector& b, + LstsqOptions opts = {}); + +} // namespace linalgebra + +namespace { + +// Solve A^T z = rhs given a pre-computed LU factorization of A. +// PA = LU => A^T = U^T L^T P +// Steps: (1) solve U^T q = rhs, (2) solve L^T w = q, (3) z[perm[i]] = w[i] +linalgebra::Vector solve_transpose(const linalgebra::LUResult& lu, + const linalgebra::Vector& rhs) { + const std::size_t n = lu.L.rows(); + + // Solve U^T q = rhs (U^T is lower triangular) + const linalgebra::Matrix Ut = linalgebra::transpose(lu.U); + const linalgebra::Vector q = linalgebra::forward_substitution(Ut, rhs); + + // Solve L^T w = q (L^T is upper triangular, unit diagonal) + const linalgebra::Matrix Lt = linalgebra::transpose(lu.L); + const linalgebra::Vector w = linalgebra::backward_substitution(Lt, q, 1e-14, true); + + // Apply inverse permutation: z[perm[i]] = w[i] + linalgebra::Vector z(n); + for (std::size_t i = 0; i < n; ++i) z[lu.perm[i]] = w[i]; + return z; +} + +} // namespace + +namespace linalgebra { + +double condition_number_1norm(const Matrix& A, double singular_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "condition_number_1norm requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + + // Exact 1-norm of A: max column sum of absolute values. + double norm_A = 0.0; + for (std::size_t j = 0; j < n; ++j) { + double col_sum = 0.0; + for (std::size_t i = 0; i < n; ++i) col_sum += std::abs(A(i, j)); + norm_A = std::max(norm_A, col_sum); + } + + const LUResult lu = lu_factor(A, singular_tolerance); + + // Estimate ||A^{-1}||_1 via the power-iteration method (Hager 1984). + // Start with x = [1/n, ..., 1/n]. + Vector x(n, 1.0 / static_cast(n)); + double est = 0.0; + + for (int iter = 0; iter < 5; ++iter) { + const Vector y = lu_solve(lu, x); // y = A^{-1} x + + // 1-norm of y. + double y1 = 0.0; + for (std::size_t i = 0; i < n; ++i) y1 += std::abs(y[i]); + + if (y1 <= est) break; + est = y1; + + Vector xi(n); + for (std::size_t i = 0; i < n; ++i) xi[i] = (y[i] >= 0.0) ? 1.0 : -1.0; + + // z = A^{-T} xi + const Vector z = solve_transpose(lu, xi); + + // Find the index maximizing |z[j]|. + std::size_t j_max = 0; + double max_z = std::abs(z[0]); + for (std::size_t i = 1; i < n; ++i) { + if (std::abs(z[i]) > max_z) { + max_z = std::abs(z[i]); + j_max = i; + } + } + + // Convergence check. + double xz = 0.0; + for (std::size_t i = 0; i < n; ++i) xz += std::abs(z[i]) / static_cast(n); + if (max_z <= xz) break; + + // New starting vector: e_{j_max}. + x.fill(0.0); + x[j_max] = 1.0; + } + + return norm_A * est; +} + +JacobiPrecond precond_jacobi(const Matrix& A, double zero_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "precond_jacobi requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + Vector inv_diag(n); + for (std::size_t i = 0; i < n; ++i) { + if (std::abs(A(i, i)) <= zero_tolerance) { + std::ostringstream oss; + oss << "precond_jacobi: zero diagonal entry at index " << i; + throw SingularMatrixError(oss.str()); + } + inv_diag[i] = 1.0 / A(i, i); + } + return JacobiPrecond{std::move(inv_diag)}; +} + +Vector apply(const JacobiPrecond& P, const Vector& x) { + const std::size_t n = x.size(); + if (n != P.inv_diag.size()) { + throw DimensionMismatchError("apply(JacobiPrecond): size mismatch"); + } + Vector result(n); + for (std::size_t i = 0; i < n; ++i) result[i] = P.inv_diag[i] * x[i]; + return result; +} + +ILU0Precond precond_ilu0(const Matrix& A, double zero_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "precond_ilu0 requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + Matrix LU = A; + + for (std::size_t k = 0; k < n; ++k) { + if (std::abs(LU(k, k)) <= zero_tolerance) { + std::ostringstream oss; + oss << "precond_ilu0: near-zero pivot at step " << k; + throw SingularMatrixError(oss.str()); + } + for (std::size_t i = k + 1; i < n; ++i) { + LU(i, k) /= LU(k, k); + for (std::size_t j = k + 1; j < n; ++j) { + LU(i, j) -= LU(i, k) * LU(k, j); + } + } + } + return ILU0Precond{std::move(LU)}; +} + +Vector apply(const ILU0Precond& P, const Vector& b) { + const std::size_t n = P.LU.rows(); + if (b.size() != n) { + throw DimensionMismatchError("apply(ILU0Precond): size mismatch"); + } + + // Extract L (unit lower) and U (upper) from combined storage. + Matrix L = Matrix::zeros(n, n); + Matrix U = Matrix::zeros(n, n); + for (std::size_t i = 0; i < n; ++i) { + L(i, i) = 1.0; + for (std::size_t j = 0; j < i; ++j) L(i, j) = P.LU(i, j); + for (std::size_t j = i; j < n; ++j) U(i, j) = P.LU(i, j); + } + + const Vector y = forward_substitution(L, b, 1e-14, true); + return backward_substitution(U, y); +} + +LstsqResult lstsq(const Matrix& A, const Vector& b, LstsqOptions opts) { + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + + if (m < n) { + std::ostringstream oss; + oss << "lstsq requires rows >= cols, got " << m << "x" << n; + throw DimensionMismatchError(oss.str()); + } + if (b.size() != m) { + std::ostringstream oss; + oss << "lstsq: rhs size " << b.size() << " does not match rows " << m; + throw DimensionMismatchError(oss.str()); + } + + const QRColPivResult qr = qr_colpiv(A, opts.rank_tolerance); + const std::size_t r = qr.rank; + + // c = Q^T b (Q is m×n with orthonormal columns) + Vector c(n, 0.0); + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < m; ++i) d += qr.Q(i, j) * b[i]; + c[j] = d; + } + + Vector x(n, 0.0); + + if (r > 0) { + Matrix Rr(r, r); + for (std::size_t i = 0; i < r; ++i) + for (std::size_t j = 0; j < r; ++j) + Rr(i, j) = qr.R(i, j); + + Vector cr(r); + for (std::size_t i = 0; i < r; ++i) cr[i] = c[i]; + + const Vector y = backward_substitution(Rr, cr); + + // Permuted solution: x_perm[0:r] = y, x_perm[r:n] = 0, then un-permute. + for (std::size_t j = 0; j < r; ++j) x[qr.perm[j]] = y[j]; + } + + const Vector res = A * x - b; + const double residual_norm = norm2(res); + + return LstsqResult{std::move(x), r, residual_norm}; +} + +} // namespace linalgebra -- cgit v1.2.3