Finish TODOs

7 files changed, 1390 insertions, 4 deletions

diff --git a/src/expm.cpp b/src/expm.cpp
new file mode 100644
index 0000000..71bdbd9
--- /dev/null
+++ b/src/expm.cpp
@@ -0,0 +1,171 @@
+export module linalgebra:expm;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :lu;
+
+// References used throughout this file:
+//   Higham — "Functions of Matrices: Theory and Computation" (SIAM 2008)
+//   Al-Mohy & Higham (2009), SIAM J. Matrix Anal. Appl. 31(3)
+
+export namespace linalgebra {
+
+// Matrix exponential via Padé [13/13] approximant + scaling and squaring.
+[[nodiscard]] Matrix expm(const Matrix& A);
+
+}  // namespace linalgebra
+
+namespace {
+
+// Padé [13/13] coefficients, Higham (2008) Table 10.4.
+constexpr double pade_b[14] = {
+    64764752532480000.0,
+    32382376266240000.0,
+     7771770303897600.0,
+     1187353796428800.0,
+      129060195264000.0,
+       10559470521600.0,
+         670442572800.0,
+          33522128640.0,
+           1323241920.0,
+             40840800.0,
+               960960.0,
+                16380.0,
+                  182.0,
+                    1.0,
+};
+
+// theta_13: ||A||_1 threshold below which no scaling is needed.
+constexpr double theta_13 = 5.371920351148152;
+
+// 1-norm of a matrix.
+double one_norm(const linalgebra::Matrix& M) {
+    const std::size_t n = M.cols();
+    const std::size_t m = M.rows();
+    double result = 0.0;
+    for (std::size_t j = 0; j < n; ++j) {
+        double col_sum = 0.0;
+        for (std::size_t i = 0; i < m; ++i) col_sum += std::abs(M(i, j));
+        result = std::max(result, col_sum);
+    }
+    return result;
+}
+
+// Evaluate the Padé [13/13] numerator U and denominator V for matrix B.
+// Uses the factored evaluation from Higham Algorithm 10.20.
+//
+//   W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
+//   W2 = b[7]*A6  + b[5]*A4  + b[3]*A2 + b[1]*I
+//   Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
+//   Z2 = b[6]*A6  + b[4]*A4  + b[2]*A2 + b[0]*I
+//   W  = A6*W1 + W2
+//   U  = B * W
+//   V  = A6*Z1 + Z2
+//
+// expm(B) ≈ (V - U)^{-1} * (V + U)
+std::pair<linalgebra::Matrix, linalgebra::Matrix>
+pade13(const linalgebra::Matrix& B) {
+    const std::size_t n = B.rows();
+    const linalgebra::Matrix I = linalgebra::Matrix::identity(n);
+
+    const linalgebra::Matrix A2 = B * B;
+    const linalgebra::Matrix A4 = A2 * A2;
+    const linalgebra::Matrix A6 = A2 * A4;
+
+    // W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
+    linalgebra::Matrix W1(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            W1(i, j) = pade_b[13] * A6(i, j) + pade_b[11] * A4(i, j) + pade_b[9] * A2(i, j);
+
+    // W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*I
+    linalgebra::Matrix W2(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            W2(i, j) = pade_b[7] * A6(i, j) + pade_b[5] * A4(i, j)
+                     + pade_b[3] * A2(i, j) + pade_b[1] * I(i, j);
+
+    // Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
+    linalgebra::Matrix Z1(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            Z1(i, j) = pade_b[12] * A6(i, j) + pade_b[10] * A4(i, j) + pade_b[8] * A2(i, j);
+
+    // Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*I
+    linalgebra::Matrix Z2(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            Z2(i, j) = pade_b[6] * A6(i, j) + pade_b[4] * A4(i, j)
+                     + pade_b[2] * A2(i, j) + pade_b[0] * I(i, j);
+
+    // W = A6*W1 + W2
+    const linalgebra::Matrix W = A6 * W1 + W2;
+
+    // U = B * W  (numerator)
+    const linalgebra::Matrix U = B * W;
+
+    // V = A6*Z1 + Z2  (denominator)
+    const linalgebra::Matrix V = A6 * Z1 + Z2;
+
+    return {U, V};
+}
+
+}  // namespace
+
+namespace linalgebra {
+
+Matrix expm(const Matrix& A) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "expm requires a square matrix, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    const std::size_t n = A.rows();
+    if (n == 0) return Matrix::identity(0);
+
+    // Determine scaling factor s such that ||A / 2^s||_1 <= theta_13.
+    const double norm_A = one_norm(A);
+    int s = 0;
+    if (norm_A > theta_13) {
+        s = static_cast<int>(std::ceil(std::log2(norm_A / theta_13)));
+        if (s < 0) s = 0;
+    }
+
+    // Scale B = A / 2^s.
+    const double scale = 1.0 / std::ldexp(1.0, s);
+    Matrix B(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            B(i, j) = A(i, j) * scale;
+
+    // Compute Padé [13/13] approximant: R = (V - U)^{-1} * (V + U).
+    auto [U, V] = pade13(B);
+
+    // Numerator = V + U,  Denominator = V - U.
+    Matrix Numerator(n, n, 0.0);
+    Matrix Denominator(n, n, 0.0);
+    for (std::size_t i = 0; i < n; ++i) {
+        for (std::size_t j = 0; j < n; ++j) {
+            Numerator(i, j)   = V(i, j) + U(i, j);
+            Denominator(i, j) = V(i, j) - U(i, j);
+        }
+    }
+
+    const LUResult lu_denom = lu_factor(Denominator);
+
+    Matrix R(n, n, 0.0);
+    for (std::size_t j = 0; j < n; ++j) {
+        Vector col(n);
+        for (std::size_t i = 0; i < n; ++i) col[i] = Numerator(i, j);
+        const Vector sol = lu_solve(lu_denom, col);
+        for (std::size_t i = 0; i < n; ++i) R(i, j) = sol[i];
+    }
+
+    for (int i = 0; i < s; ++i) R = R * R;
+
+    return R;
+}
+
+}  // namespace linalgebra
diff --git a/src/iterative.cpp b/src/iterative.cpp
new file mode 100644
index 0000000..a0ff41d
--- /dev/null
+++ b/src/iterative.cpp
@@ -0,0 +1,404 @@
+export module linalgebra:iterative;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :norms;
+import :triangular_solve;
+
+// References used throughout this file:
+//   T&B  — Trefethen & Bau, "Numerical Linear Algebra"
+//   GVL  — Golub & Van Loan, "Matrix Computations" 4th ed.
+
+export namespace linalgebra {
+
+// Arnoldi iteration builds an orthonormal Krylov basis and the
+// corresponding upper Hessenberg matrix.
+// Reference: T&B Algorithm 33.1; GVL §6.3
+
+struct ArnoldiResult {
+    Matrix Q;          // n × (steps_taken + 1) orthonormal columns
+    Matrix H;          // (steps_taken + 1) × steps_taken upper Hessenberg
+    int    steps_taken;
+    bool   breakdown;  // true if invariant subspace found early
+};
+
+struct ArnoldiOptions {
+    double breakdown_tolerance = 1e-14;
+};
+
+[[nodiscard]] ArnoldiResult arnoldi(const Matrix& A, const Vector& b, int k,
+                                    ArnoldiOptions opts = {});
+
+// Conjugate Gradient — for symmetric positive definite systems
+// Reference: T&B Algorithm 38.1; GVL §11.3
+
+struct CGOptions {
+    double tolerance     = 1e-10;
+    int    max_iterations = 1000;
+};
+
+struct CGResult {
+    Vector x;
+    int    iterations;
+    double final_residual;
+};
+
+[[nodiscard]] CGResult solve_cg(const Matrix& A, const Vector& b, CGOptions opts = {});
+
+// Restarted GMRES — for general square systems
+// Reference: T&B Algorithm 35.1; GVL §11.4.2
+
+struct GMRESOptions {
+    double tolerance      = 1e-10;
+    int    max_iterations = 200;
+    int    restart        = 50;
+};
+
+struct GMRESResult {
+    Vector x;
+    int    iterations;
+    double final_residual;
+};
+
+[[nodiscard]] GMRESResult solve_gmres(const Matrix& A, const Vector& b,
+                                      GMRESOptions opts = {});
+
+// BiCGSTAB — for general square systems (van der Vorst 1992)
+// Reference: GVL §11.5.3
+
+struct BiCGSTABOptions {
+    double tolerance     = 1e-10;
+    int    max_iterations = 1000;
+};
+
+struct BiCGSTABResult {
+    Vector x;
+    int    iterations;
+    double final_residual;
+};
+
+[[nodiscard]] BiCGSTABResult solve_bicgstab(const Matrix& A, const Vector& b,
+                                             BiCGSTABOptions opts = {});
+
+}  // namespace linalgebra
+
+namespace linalgebra {
+
+ArnoldiResult arnoldi(const Matrix& A, const Vector& b, int k, ArnoldiOptions opts) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "arnoldi requires a square matrix, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+    const std::size_t n = A.rows();
+    if (b.size() != n) {
+        std::ostringstream oss;
+        oss << "arnoldi: b size " << b.size() << " must match matrix dimension " << n;
+        throw DimensionMismatchError(oss.str());
+    }
+    if (k <= 0) throw std::invalid_argument("arnoldi: k must be >= 1");
+
+    const auto kk = static_cast<std::size_t>(k);
+
+    Matrix Q(n, kk + 1, 0.0);
+    Matrix H(kk + 1, kk, 0.0);
+
+    const double b_norm = norm2(b);
+    if (b_norm == 0.0) {
+        return ArnoldiResult{std::move(Q), std::move(H), 0, true};
+    }
+
+    // q_0 = b / ||b||
+    for (std::size_t i = 0; i < n; ++i) Q(i, 0) = b[i] / b_norm;
+
+    int steps = 0;
+    for (std::size_t j = 0; j < kk; ++j) {
+        // z = A * Q[:, j]
+        Vector qj(n);
+        for (std::size_t i = 0; i < n; ++i) qj[i] = Q(i, j);
+        Vector z = A * qj;
+
+        // Modified Gram-Schmidt orthogonalization.
+        for (std::size_t i = 0; i <= j; ++i) {
+            double h = 0.0;
+            for (std::size_t row = 0; row < n; ++row) h += Q(row, i) * z[row];
+            H(i, j) = h;
+            for (std::size_t row = 0; row < n; ++row) z[row] -= h * Q(row, i);
+        }
+
+        const double z_norm = norm2(z);
+        H(j + 1, j) = z_norm;
+        ++steps;
+
+        if (z_norm < opts.breakdown_tolerance) {
+            // Lucky breakdown: invariant subspace found.
+            return ArnoldiResult{std::move(Q), std::move(H), steps, true};
+        }
+
+        for (std::size_t i = 0; i < n; ++i) Q(i, j + 1) = z[i] / z_norm;
+    }
+
+    return ArnoldiResult{std::move(Q), std::move(H), steps, false};
+}
+
+CGResult solve_cg(const Matrix& A, const Vector& b, CGOptions opts) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "solve_cg requires a square matrix, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+    const std::size_t n = A.rows();
+    if (b.size() != n) {
+        std::ostringstream oss;
+        oss << "solve_cg: rhs size " << b.size() << " does not match dimension " << n;
+        throw DimensionMismatchError(oss.str());
+    }
+
+    // Check symmetry.
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = i + 1; j < n; ++j)
+            if (std::abs(A(i, j) - A(j, i)) > 1e-10 * (std::abs(A(i, j)) + 1.0))
+                throw LinAlgError("solve_cg: matrix is not symmetric");
+
+    Vector x(n, 0.0);
+    Vector r = b;  // r = b - A*x0, x0 = 0
+    Vector p = r;
+    double rr = dot(r, r);
+
+    if (std::sqrt(rr) < opts.tolerance) {
+        return CGResult{std::move(x), 0, std::sqrt(rr)};
+    }
+
+    for (int iter = 1; iter <= opts.max_iterations; ++iter) {
+        const Vector Ap = A * p;
+        const double pAp = dot(p, Ap);
+        if (std::abs(pAp) == 0.0) break;
+        const double alpha = rr / pAp;
+
+        for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i];
+        for (std::size_t i = 0; i < n; ++i) r[i] -= alpha * Ap[i];
+
+        const double rr_new = dot(r, r);
+        const double res = std::sqrt(rr_new);
+
+        if (res < opts.tolerance) {
+            return CGResult{std::move(x), iter, res};
+        }
+
+        const double beta = rr_new / rr;
+        for (std::size_t i = 0; i < n; ++i) p[i] = r[i] + beta * p[i];
+        rr = rr_new;
+    }
+
+    const double final_res = norm2(A * x - b);
+    std::ostringstream oss;
+    oss << "solve_cg: did not converge in " << opts.max_iterations
+        << " iterations (final residual = " << final_res << ")";
+    throw NonConvergenceError(oss.str());
+}
+
+GMRESResult solve_gmres(const Matrix& A, const Vector& b, GMRESOptions opts) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "solve_gmres requires a square matrix, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+    const std::size_t n = A.rows();
+    if (b.size() != n) {
+        std::ostringstream oss;
+        oss << "solve_gmres: rhs size " << b.size() << " does not match dimension " << n;
+        throw DimensionMismatchError(oss.str());
+    }
+
+    const int m = std::min(opts.restart, static_cast<int>(n));
+
+    Vector x(n, 0.0);
+    int total_iters = 0;
+
+    while (total_iters < opts.max_iterations) {
+        // Compute residual r = b - A*x.
+        const Vector Ax = A * x;
+        Vector r(n);
+        for (std::size_t i = 0; i < n; ++i) r[i] = b[i] - Ax[i];
+        const double beta = norm2(r);
+
+        if (beta < opts.tolerance) {
+            return GMRESResult{std::move(x), total_iters, beta};
+        }
+
+        // Arnoldi to build Krylov basis.
+        const ArnoldiResult ar = arnoldi(A, r, m);
+        const int steps = ar.steps_taken;
+
+        if (steps == 0) break;
+
+        // Solve least-squares problem: min ||beta*e1 - H_hat * y||
+        // where H_hat is (steps+1) x steps.
+        // Apply Givens rotations to reduce H_hat to upper triangular.
+        const auto sz = static_cast<std::size_t>(steps);
+
+        // Work on a copy of the relevant submatrix of H and rhs g.
+        std::vector<std::vector<double>> Hwork(sz + 1, std::vector<double>(sz, 0.0));
+        for (std::size_t i = 0; i <= sz; ++i)
+            for (std::size_t j = 0; j < sz; ++j)
+                Hwork[i][j] = ar.H(i, j);
+
+        std::vector<double> g(sz + 1, 0.0);
+        g[0] = beta;
+
+        // Accumulated Givens rotations.
+        std::vector<double> cs(sz), sn(sz);
+
+        for (std::size_t j = 0; j < sz; ++j) {
+            // Givens to zero H[j+1, j].
+            const double f = Hwork[j][j];
+            const double hh = Hwork[j + 1][j];
+            const double r_val = std::hypot(f, hh);
+            if (r_val == 0.0) { cs[j] = 1.0; sn[j] = 0.0; continue; }
+            cs[j] =  f / r_val;
+            sn[j] = hh / r_val;
+
+            // Apply to column j of H (only the two relevant rows).
+            Hwork[j][j]     =  cs[j] * f + sn[j] * hh;
+            Hwork[j + 1][j] = 0.0;
+
+            // Apply to remaining columns.
+            for (std::size_t l = j + 1; l < sz; ++l) {
+                const double t0 = Hwork[j][l];
+                const double t1 = Hwork[j + 1][l];
+                Hwork[j][l]     =  cs[j] * t0 + sn[j] * t1;
+                Hwork[j + 1][l] = -sn[j] * t0 + cs[j] * t1;
+            }
+
+            // Apply to g.
+            const double g0 = g[j];
+            const double g1 = g[j + 1];
+            g[j]     =  cs[j] * g0 + sn[j] * g1;
+            g[j + 1] = -sn[j] * g0 + cs[j] * g1;
+        }
+
+        // Backward substitution: solve the sz×sz upper triangular system.
+        std::vector<double> y(sz, 0.0);
+        for (int i = static_cast<int>(sz) - 1; i >= 0; --i) {
+            double sum = g[static_cast<std::size_t>(i)];
+            for (std::size_t j = static_cast<std::size_t>(i) + 1; j < sz; ++j)
+                sum -= Hwork[static_cast<std::size_t>(i)][j] * y[j];
+            if (std::abs(Hwork[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)]) > 0.0)
+                y[static_cast<std::size_t>(i)] =
+                    sum / Hwork[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)];
+        }
+
+        // Update x = x + Q[:, 0:steps] * y.
+        for (std::size_t j = 0; j < sz; ++j) {
+            for (std::size_t i = 0; i < n; ++i) {
+                x[i] += y[j] * ar.Q(i, j);
+            }
+        }
+
+        ++total_iters;
+
+        // Check convergence.
+        const double final_res = std::abs(g[sz]);
+        if (final_res < opts.tolerance || ar.breakdown) {
+            return GMRESResult{std::move(x), total_iters, final_res};
+        }
+    }
+
+    const double final_res = norm2(A * x - b);
+    std::ostringstream oss;
+    oss << "solve_gmres: did not converge in " << opts.max_iterations
+        << " restarts (final residual = " << final_res << ")";
+    throw NonConvergenceError(oss.str());
+}
+
+BiCGSTABResult solve_bicgstab(const Matrix& A, const Vector& b, BiCGSTABOptions opts) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "solve_bicgstab requires a square matrix, got "
+            << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+    const std::size_t n = A.rows();
+    if (b.size() != n) {
+        std::ostringstream oss;
+        oss << "solve_bicgstab: rhs size " << b.size() << " does not match dimension " << n;
+        throw DimensionMismatchError(oss.str());
+    }
+
+    Vector x(n, 0.0);
+    Vector r = b;  // r = b - A*x0, x0 = 0
+    Vector r_hat = r;  // shadow residual, fixed throughout
+
+    double rho_old = 1.0, alpha = 1.0, omega = 1.0;
+    Vector v(n, 0.0), p(n, 0.0);
+
+    for (int iter = 1; iter <= opts.max_iterations; ++iter) {
+        const double rho_new = dot(r_hat, r);
+
+        if (std::abs(rho_new) < 1e-300) {
+            std::ostringstream oss;
+            oss << "solve_bicgstab: breakdown (rho near zero) at iteration " << iter;
+            throw NonConvergenceError(oss.str());
+        }
+
+        const double beta = (rho_new / rho_old) * (alpha / omega);
+
+        for (std::size_t i = 0; i < n; ++i)
+            p[i] = r[i] + beta * (p[i] - omega * v[i]);
+
+        v = A * p;
+
+        const double denom_alpha = dot(r_hat, v);
+        if (std::abs(denom_alpha) < 1e-300) {
+            std::ostringstream oss;
+            oss << "solve_bicgstab: breakdown (r_hat·v near zero) at iteration " << iter;
+            throw NonConvergenceError(oss.str());
+        }
+        alpha = rho_new / denom_alpha;
+
+        Vector s(n);
+        for (std::size_t i = 0; i < n; ++i) s[i] = r[i] - alpha * v[i];
+
+        const double s_norm = norm2(s);
+        if (s_norm < opts.tolerance) {
+            for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i];
+            return BiCGSTABResult{std::move(x), iter, s_norm};
+        }
+
+        const Vector t = A * s;
+        const double tt = dot(t, t);
+
+        if (std::abs(tt) < 1e-300) {
+            std::ostringstream oss;
+            oss << "solve_bicgstab: breakdown (t·t near zero) at iteration " << iter;
+            throw NonConvergenceError(oss.str());
+        }
+
+        omega = dot(t, s) / tt;
+
+        if (std::abs(omega) < 1e-300) {
+            std::ostringstream oss;
+            oss << "solve_bicgstab: breakdown (omega near zero) at iteration " << iter;
+            throw NonConvergenceError(oss.str());
+        }
+
+        for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i] + omega * s[i];
+        for (std::size_t i = 0; i < n; ++i) r[i] = s[i] - omega * t[i];
+
+        const double r_norm = norm2(r);
+        if (r_norm < opts.tolerance) {
+            return BiCGSTABResult{std::move(x), iter, r_norm};
+        }
+
+        rho_old = rho_new;
+    }
+
+    const double final_res = norm2(A * x - b);
+    std::ostringstream oss;
+    oss << "solve_bicgstab: did not converge in " << opts.max_iterations
+        << " iterations (final residual = " << final_res << ")";
+    throw NonConvergenceError(oss.str());
+}
+
+}  // namespace linalgebra
diff --git a/src/linalgebra.cpp b/src/linalgebra.cpp
index 18e9cf3..d6d5652 100644
--- a/src/linalgebra.cpp
+++ b/src/linalgebra.cpp
@@ -9,3 +9,8 @@ export import :lu;
 export import :qr;
 export import :qr_iteration;
 export import :cholesky;
+export import :sym_eigen;
+export import :svd;
+export import :iterative;
+export import :precond;
+export import :expm;
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
diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp
index f8a3be5..7f50097 100644
--- a/src/qr_iteration.cpp
+++ b/src/qr_iteration.cpp
@@ -14,9 +14,7 @@ import :qr;
 
 export namespace linalgebra {
 
-// ---------------------------------------------------------------------------
 // Options
-// ---------------------------------------------------------------------------
 
 struct QRIterationOptions {
     double tolerance = 1e-10;
@@ -509,9 +507,7 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opt
     return result;
 }
 
-// ---------------------------------------------------------------------------
 // Francis double-shift QR with implicit bulge chasing
-// ---------------------------------------------------------------------------
 
 QRIterationResult eigenvalues_francis(const Matrix& A, QRIterationOptions opts) {
     require_square(A, "eigenvalues_francis");
diff --git a/src/svd.cpp b/src/svd.cpp
new file mode 100644
index 0000000..8157845
--- /dev/null
+++ b/src/svd.cpp
@@ -0,0 +1,323 @@
+export module linalgebra:svd;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :qr_iteration;
+
+// References used throughout this file:
+//   GVL  — Golub & Van Loan, "Matrix Computations" 4th ed.
+//   T&B  — Trefethen & Bau, "Numerical Linear Algebra"
+
+export namespace linalgebra {
+
+struct SVDResult {
+    Matrix U;      // m × m orthogonal (left singular vectors)
+    Vector sigma;  // min(m,n) singular values, sorted descending
+    Matrix Vt;     // n × n orthogonal (V^T, right singular vectors transposed)
+};
+
+struct SVDOptions {
+    double tolerance     = 1e-12;
+    int    max_iterations = 1000;
+};
+
+// Golub-Kahan bidiagonalization + Golub-Reinsch QR sweeps.
+// Requires rows >= cols; throws DimensionMismatchError otherwise.
+// Reference: GVL §8.6; T&B Lecture 31.
+[[nodiscard]] SVDResult svd(const Matrix& A, SVDOptions opts = {});
+
+}  // namespace linalgebra
+
+namespace linalgebra {
+
+SVDResult svd(const Matrix& A, SVDOptions opts) {
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+
+    if (m < n) {
+        std::ostringstream oss;
+        oss << "svd requires rows >= cols, got " << m << "x" << n;
+        throw DimensionMismatchError(oss.str());
+    }
+    if (n == 0) {
+        return SVDResult{Matrix::identity(m), Vector(0), Matrix::identity(0)};
+    }
+
+    // GK bidiagonalization  A = U_acc * B * Vt_acc
+    // where B is m×n upper bidiagonal (nonzero on diagonal and superdiagonal).
+
+    Matrix work  = A;
+    Matrix U_acc = Matrix::identity(m);
+    Matrix Vt_acc = Matrix::identity(n);
+
+    for (std::size_t k = 0; k < n; ++k) {
+        const std::size_t p_left = m - k;  // column length below row k
+
+        // Left Householder: zero work[k+1:, k].
+        {
+            std::vector<double> u(p_left);
+            for (std::size_t i = 0; i < p_left; ++i) u[i] = work(k + i, k);
+
+            double x_norm = 0.0;
+            for (double v : u) x_norm += v * v;
+            x_norm = std::sqrt(x_norm);
+
+            if (x_norm > 0.0) {
+                const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm;
+                u[0] += sigma;
+                double utu = 0.0;
+                for (double v : u) utu += v * v;
+                const double tau = 2.0 / utu;
+
+                // Apply to work from left: work[k:, k:] -= tau * u * (u^T work[k:, k:])
+                for (std::size_t j = k; j < n; ++j) {
+                    double d = 0.0;
+                    for (std::size_t i = 0; i < p_left; ++i) d += u[i] * work(k + i, j);
+                    const double c = tau * d;
+                    for (std::size_t i = 0; i < p_left; ++i) work(k + i, j) -= c * u[i];
+                }
+
+                // Accumulate into U_acc from right: U_acc[:, k:] -= tau * (U_acc[:, k:] u) * u^T
+                for (std::size_t j = 0; j < m; ++j) {
+                    double d = 0.0;
+                    for (std::size_t i = 0; i < p_left; ++i) d += U_acc(j, k + i) * u[i];
+                    const double c = tau * d;
+                    for (std::size_t i = 0; i < p_left; ++i) U_acc(j, k + i) -= c * u[i];
+                }
+            }
+        }
+
+        // Right Householder: zero work[k, k+2:].
+        if (k + 2 <= n) {
+            const std::size_t p_right = n - k - 1;  // row length after col k+1
+
+            std::vector<double> v(p_right);
+            for (std::size_t j = 0; j < p_right; ++j) v[j] = work(k, k + 1 + j);
+
+            double x_norm = 0.0;
+            for (double val : v) x_norm += val * val;
+            x_norm = std::sqrt(x_norm);
+
+            if (x_norm > 0.0) {
+                const double sigma = (v[0] >= 0.0 ? 1.0 : -1.0) * x_norm;
+                v[0] += sigma;
+                double vtv = 0.0;
+                for (double val : v) vtv += val * val;
+                const double tau = 2.0 / vtv;
+
+                // Apply to work from right: work[:, k+1:] -= tau * (work[:, k+1:] v) * v^T
+                for (std::size_t i = k; i < m; ++i) {
+                    double d = 0.0;
+                    for (std::size_t j = 0; j < p_right; ++j) d += work(i, k + 1 + j) * v[j];
+                    const double c = tau * d;
+                    for (std::size_t j = 0; j < p_right; ++j) work(i, k + 1 + j) -= c * v[j];
+                }
+
+                // Accumulate into Vt_acc from left: Vt_acc[k+1:, :] -= tau * v * (v^T Vt_acc[k+1:, :])
+                for (std::size_t j = 0; j < n; ++j) {
+                    double d = 0.0;
+                    for (std::size_t i = 0; i < p_right; ++i) d += v[i] * Vt_acc(k + 1 + i, j);
+                    const double c = tau * d;
+                    for (std::size_t i = 0; i < p_right; ++i) Vt_acc(k + 1 + i, j) -= c * v[i];
+                }
+            }
+        }
+    }
+
+    // Extract bidiagonal: d[i] = diag, e[i] = superdiag.
+    std::vector<double> d(n), e(n > 1 ? n - 1 : 0, 0.0);
+    for (std::size_t i = 0; i < n; ++i) d[i] = work(i, i);
+    for (std::size_t i = 0; i + 1 < n; ++i) e[i] = work(i, i + 1);
+
+    // GR QR sweeps on the bidiagonal.
+    // Accumulate left rotations into U_f, right rotations into Vt_f.
+
+    Matrix U_f  = Matrix::identity(n);
+    Matrix Vt_f = Matrix::identity(n);
+
+    std::size_t active = n;
+    int total_iters = 0;
+
+    while (active > 1) {
+        if (total_iters >= opts.max_iterations) {
+            std::ostringstream oss;
+            oss << "svd: did not converge in " << opts.max_iterations << " iterations";
+            throw NonConvergenceError(oss.str());
+        }
+        ++total_iters;
+
+        // Deflate small superdiagonals.
+        while (active > 1) {
+            const std::size_t i = active - 2;
+            if (std::abs(e[i]) <= opts.tolerance * (std::abs(d[i]) + std::abs(d[active - 1]))) {
+                e[i] = 0.0;
+                --active;
+            } else {
+                break;
+            }
+        }
+        if (active <= 1) break;
+
+        // Handle zero on diagonal: if d[k] == 0 for k < active-1, chase the
+        // nonzero e[k] to zero using a sequence of left Givens rotations.
+        bool zero_diag = false;
+        for (std::size_t k = 0; k + 1 < active; ++k) {
+            if (std::abs(d[k]) <= opts.tolerance) {
+                zero_diag = true;
+                // Chase e[k] to zero using left Givens rotations in rows k and k+1..active-1.
+                double f = e[k];
+                e[k] = 0.0;
+                for (std::size_t j = k + 1; j < active && f != 0.0; ++j) {
+                    const double g = d[j];
+                    const double r = std::hypot(f, g);
+                    const double c = g / r;
+                    const double s = -f / r;
+                    d[j] = r;
+                    if (j + 1 < active) {
+                        f = s * e[j];
+                        e[j] *= c;
+                    }
+                    // Accumulate into U_f (left rotation on rows k and j).
+                    // Apply a rotation between rows k and j manually on U_f columns.
+                    for (std::size_t col = 0; col < n; ++col) {
+                        const double uk = U_f(k, col);
+                        const double uj = U_f(j, col);
+                        U_f(k, col) =  c * uk - s * uj;
+                        U_f(j, col) =  s * uk + c * uj;
+                    }
+                }
+                break;
+            }
+        }
+        if (zero_diag) continue;
+
+        // Wilkinson shift from bottom 2×2 of B^T B.
+        // B^T B bottom-right 2×2 (indices active-2, active-1):
+        //   [d[a-2]^2 + e[a-3]^2,  d[a-2]*e[a-2]]
+        //   [d[a-2]*e[a-2],          d[a-1]^2 + e[a-2]^2]  (if a>=2)
+        const double a = active >= 2 ? d[active - 2] : 0.0;
+        const double b_val = active >= 2 ? e[active - 2] : 0.0;
+        const double c_val = d[active - 1];
+        const double t11 = a * a + (active >= 3 ? e[active - 3] * e[active - 3] : 0.0);
+        const double t12 = a * b_val;
+        const double t22 = c_val * c_val + b_val * b_val;
+        const double delta = 0.5 * (t11 - t22);
+        const double denom = std::abs(delta) + std::hypot(delta, t12);
+        const double mu = (denom == 0.0) ? t22
+                        : t22 - (delta >= 0.0 ? 1.0 : -1.0) * (t12 * t12) / denom;
+
+        // Golub-Reinsch implicit QR step.
+        double f = d[0] * d[0] - mu;
+        double g = d[0] * e[0];
+
+        for (std::size_t i = 0; i + 1 < active; ++i) {
+            // Right Givens: eliminate g from (f, g) in columns i and i+1.
+            {
+                const double r  = std::hypot(f, g);
+                const double cr = (r == 0.0) ? 1.0 : f / r;
+                const double sr = (r == 0.0) ? 0.0 : g / r;
+
+                if (i > 0) e[i - 1] = r;
+
+                f =  cr * d[i] + sr * e[i];
+                e[i] = -sr * d[i] + cr * e[i];
+                g = sr * d[i + 1];
+                d[i + 1] *= cr;
+
+                // Accumulate right rotation into Vt_f (acts on cols i and i+1 of V, i.e., rows of Vt_f).
+                for (std::size_t row = 0; row < n; ++row) {
+                    const double vi  = Vt_f(i,     row);
+                    const double vi1 = Vt_f(i + 1, row);
+                    Vt_f(i,     row) =  cr * vi + sr * vi1;
+                    Vt_f(i + 1, row) = -sr * vi + cr * vi1;
+                }
+            }
+
+            // Left Givens: eliminate g from (f, g) in rows i and i+1.
+            {
+                const double r  = std::hypot(f, g);
+                const double cl = (r == 0.0) ? 1.0 : f / r;
+                const double sl = (r == 0.0) ? 0.0 : g / r;
+
+                d[i] = r;
+
+                f =  cl * e[i] + sl * d[i + 1];
+                d[i + 1] = -sl * e[i] + cl * d[i + 1];
+                e[i] = f;
+
+                if (i + 2 < active) {
+                    g = sl * e[i + 1];
+                    e[i + 1] *= cl;
+                }
+
+                // Accumulate left rotation into U_f (acts on rows i and i+1).
+                for (std::size_t col = 0; col < n; ++col) {
+                    const double ui  = U_f(i,     col);
+                    const double ui1 = U_f(i + 1, col);
+                    U_f(i,     col) =  cl * ui + sl * ui1;
+                    U_f(i + 1, col) = -sl * ui + cl * ui1;
+                }
+            }
+        }
+        // Set the last diagonal update.
+        e[active - 2] = f;
+    }
+
+    // U = U_acc * U_f^T  (U_f accumulates row operations on B,
+    // which correspond to left singular vectors relative to U_acc).
+    // Vt = Vt_f * Vt_acc  (Vt_f accumulates row ops on Vt_acc).
+    // sigma = |d[i]|, flip signs into U.
+
+    // U_f stores row-wise left rotations applied to the bidiagonal's rows.
+    // The actual left factor is U_acc * U_f^T (since each left Givens G was
+    // applied as B <- G B, meaning U_acc absorbs G^T from the right).
+    const Matrix Uf_t = transpose(U_f);
+    // U_acc is m×m, Uf_t is n×n.  We need m×m U; embed Uf_t into top-left.
+    Matrix U_full = Matrix::zeros(m, m);
+    // Copy U_acc * Uf_t into first n columns; remaining m-n columns of U_acc unchanged.
+    for (std::size_t i = 0; i < m; ++i) {
+        for (std::size_t j = 0; j < n; ++j) {
+            double s = 0.0;
+            for (std::size_t l = 0; l < n; ++l) s += U_acc(i, l) * Uf_t(l, j);
+            U_full(i, j) = s;
+        }
+        for (std::size_t j = n; j < m; ++j) U_full(i, j) = U_acc(i, j);
+    }
+
+    Matrix Vt_full = Vt_f * Vt_acc;
+
+    Vector sigma(n);
+    for (std::size_t i = 0; i < n; ++i) {
+        sigma[i] = std::abs(d[i]);
+        if (d[i] < 0.0) {
+            // Negate corresponding column of U (row of U^T) to keep sigma positive.
+            for (std::size_t j = 0; j < m; ++j) U_full(j, i) = -U_full(j, i);
+        }
+    }
+
+    // Sort singular values descending, applying same permutation to U cols and Vt rows.
+    std::vector<std::size_t> idx(n);
+    std::iota(idx.begin(), idx.end(), std::size_t{0});
+    std::sort(idx.begin(), idx.end(),
+              [&](std::size_t a, std::size_t b) { return sigma[a] > sigma[b]; });
+
+    Vector sigma_sorted(n);
+    Matrix U_sorted(m, m);
+    Matrix Vt_sorted(n, n);
+
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = n; j < m; ++j)
+            U_sorted(i, j) = U_full(i, j);
+
+    for (std::size_t rank = 0; rank < n; ++rank) {
+        const std::size_t src = idx[rank];
+        sigma_sorted[rank] = sigma[src];
+        for (std::size_t i = 0; i < m; ++i) U_sorted(i, rank) = U_full(i, src);
+        for (std::size_t j = 0; j < n; ++j) Vt_sorted(rank, j) = Vt_full(src, j);
+    }
+
+    return SVDResult{std::move(U_sorted), std::move(sigma_sorted), std::move(Vt_sorted)};
+}
+
+}  // namespace linalgebra
diff --git a/src/sym_eigen.cpp b/src/sym_eigen.cpp
new file mode 100644
index 0000000..a38b272
--- /dev/null
+++ b/src/sym_eigen.cpp
@@ -0,0 +1,208 @@
+export module linalgebra:sym_eigen;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :norms;
+import :lu;
+
+// References used throughout this file:
+//   T&B  — Trefethen & Bau, "Numerical Linear Algebra"
+//   GVL  — Golub & Van Loan, "Matrix Computations" 4th ed.
+
+export namespace linalgebra {
+
+// Symmetric tridiagonalization via Householder reflections
+// Reference: GVL §8.3.1; T&B Lecture 26
+//
+// For symmetric A, computes Q and T such that A = Q T Q^T,
+// where T is symmetric tridiagonal and Q is orthogonal.
+
+struct TridiagonalizeResult {
+    Matrix T;  // symmetric tridiagonal (full n×n)
+    Matrix Q;  // orthogonal: A = Q T Q^T
+};
+
+[[nodiscard]] TridiagonalizeResult tridiagonalize(const Matrix& A,
+                                                   double symmetry_tolerance = 1e-12);
+
+// Eigenvectors via inverse iteration
+// Reference: GVL §7.6.1; T&B Lecture 27
+//
+// Given A and a vector of approximate eigenvalues, returns a matrix whose
+// columns are the corresponding (approximate) eigenvectors.
+
+struct InverseIterationOptions {
+    double tolerance     = 1e-10;
+    int    max_iterations = 100;
+};
+
+struct InverseIterationResult {
+    Matrix eigenvectors;  // n × k, column j is eigenvector for eigenvalues[j]
+    Vector residuals;     // ||A*v_j - lambda_j*v_j||_2 for each j
+};
+
+[[nodiscard]] InverseIterationResult eigenvectors_inverse_iteration(
+    const Matrix& A,
+    const Vector& eigenvalues,
+    InverseIterationOptions opts = {});
+
+}  // namespace linalgebra
+
+namespace linalgebra {
+
+TridiagonalizeResult tridiagonalize(const Matrix& A, double symmetry_tolerance) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "tridiagonalize requires a square matrix, got "
+            << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    const std::size_t n = A.rows();
+
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = i + 1; j < n; ++j)
+            if (std::abs(A(i, j) - A(j, i)) > symmetry_tolerance)
+                throw LinAlgError("tridiagonalize requires a symmetric matrix");
+
+    Matrix T = A;
+    Matrix Q = Matrix::identity(n);
+
+    for (std::size_t k = 0; k + 2 <= n; ++k) {
+        const std::size_t p = n - k - 1;  // length of sub-column
+
+        // Form Householder vector u from T[k+1:, k].
+        std::vector<double> u(p);
+        for (std::size_t i = 0; i < p; ++i) u[i] = T(k + 1 + i, k);
+
+        double x_norm = 0.0;
+        for (double v : u) x_norm += v * v;
+        x_norm = std::sqrt(x_norm);
+
+        if (x_norm == 0.0) continue;
+
+        const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm;
+        u[0] += sigma;
+
+        double utu = 0.0;
+        for (double v : u) utu += v * v;
+        const double tau = 2.0 / utu;
+
+        // Apply H = I - tau*u*u^T symmetrically: T <- H T H
+        // Left: T[k+1:, :] -= tau * u * (u^T T[k+1:, :])
+        for (std::size_t j = 0; j < n; ++j) {
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += u[i] * T(k + 1 + i, j);
+            const double coeff = tau * d;
+            for (std::size_t i = 0; i < p; ++i) T(k + 1 + i, j) -= coeff * u[i];
+        }
+
+        // Right: T[:, k+1:] -= tau * (T[:, k+1:] u) * u^T
+        for (std::size_t j = 0; j < n; ++j) {
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += T(j, k + 1 + i) * u[i];
+            const double coeff = tau * d;
+            for (std::size_t i = 0; i < p; ++i) T(j, k + 1 + i) -= coeff * u[i];
+        }
+
+        // Explicitly zero sub-subdiagonal entries for numerical cleanliness.
+        for (std::size_t i = 1; i < p; ++i) {
+            T(k + 1 + i, k) = 0.0;
+            T(k, k + 1 + i) = 0.0;
+        }
+
+        // Accumulate Q: Q[:, k+1:] -= tau * (Q[:, k+1:] u) * u^T
+        for (std::size_t j = 0; j < n; ++j) {
+            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];
+        }
+    }
+
+    return TridiagonalizeResult{std::move(T), std::move(Q)};
+}
+
+InverseIterationResult eigenvectors_inverse_iteration(const Matrix& A,
+                                                       const Vector& eigenvalues,
+                                                       InverseIterationOptions opts) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "eigenvectors_inverse_iteration requires a square matrix, got "
+            << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    const std::size_t n = A.rows();
+    const std::size_t k = eigenvalues.size();
+
+    Matrix evecs(n, k);
+    Vector residuals(k, 0.0);
+
+    // Starting vector: uniform unit vector.
+    Vector v0(n, 1.0 / std::sqrt(static_cast<double>(n)));
+
+    for (std::size_t col = 0; col < k; ++col) {
+        double lambda = eigenvalues[col];
+
+        // Build shifted matrix B = A - lambda*I.
+        Matrix B = A;
+        for (std::size_t i = 0; i < n; ++i) B(i, i) -= lambda;
+
+        // Try to factor; if near-singular (either at factorization or solve time), perturb the shift.
+        auto make_lu = [&]() -> LUResult {
+            try {
+                LUResult result = lu_factor(B, 1e-14);
+                lu_solve(result, v0);  // probe: backward_substitution validates U diagonal
+                return result;
+            } catch (const SingularMatrixError&) {
+                lambda += 1e-7;
+                for (std::size_t i = 0; i < n; ++i) B(i, i) += 1e-7;
+                return lu_factor(B, 1e-14);
+            }
+        };
+        LUResult lu_b = make_lu();
+
+        Vector v = v0;
+
+        bool converged = false;
+        for (int iter = 0; iter < opts.max_iterations; ++iter) {
+            Vector w = lu_solve(lu_b, v);
+
+            // Normalize.
+            double w_norm = norm2(w);
+            if (w_norm == 0.0) break;
+            v = w / w_norm;
+
+            // Residual: ||A v - lambda v||.
+            const Vector Av = A * v;
+            double res = 0.0;
+            for (std::size_t i = 0; i < n; ++i) {
+                const double d = Av[i] - eigenvalues[col] * v[i];
+                res += d * d;
+            }
+            res = std::sqrt(res);
+
+            if (res < opts.tolerance) {
+                converged = true;
+                residuals[col] = res;
+                break;
+            }
+        }
+
+        if (!converged) {
+            std::ostringstream oss;
+            oss << "eigenvectors_inverse_iteration: did not converge for eigenvalue "
+                << col << " (lambda = " << eigenvalues[col] << ") in "
+                << opts.max_iterations << " iterations";
+            throw NonConvergenceError(oss.str());
+        }
+
+        for (std::size_t i = 0; i < n; ++i) evecs(i, col) = v[i];
+    }
+
+    return InverseIterationResult{std::move(evecs), std::move(residuals)};
+}
+
+}  // namespace linalgebra