1 files changed, 279 insertions, 0 deletions
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<double>(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<double>(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