Finish TODOs

7 files changed, 805 insertions, 4 deletions

diff --git a/tests/test_expm.cpp b/tests/test_expm.cpp
new file mode 100644
index 0000000..a0d3ea0
--- /dev/null
+++ b/tests/test_expm.cpp
@@ -0,0 +1,144 @@
+import linalgebra;
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <numbers>
+#include <random>
+
+using linalgebra::DimensionMismatchError;
+using linalgebra::Matrix;
+
+namespace {
+
+double frobenius_diff(const Matrix& A, const Matrix& B) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < A.rows(); ++i)
+        for (std::size_t j = 0; j < A.cols(); ++j) {
+            const double d = A(i, j) - B(i, j);
+            s += d * d;
+        }
+    return std::sqrt(s);
+}
+
+// Taylor series
+Matrix expm_taylor5(const Matrix& A) {
+    const std::size_t n = A.rows();
+    Matrix result = Matrix::identity(n);
+    Matrix power  = Matrix::identity(n);
+    double fact = 1.0;
+    for (int k = 1; k <= 5; ++k) {
+        power = power * A;
+        fact *= static_cast<double>(k);
+        const double inv_fact = 1.0 / fact;
+        for (std::size_t i = 0; i < n; ++i)
+            for (std::size_t j = 0; j < n; ++j)
+                result(i, j) += inv_fact * power(i, j);
+    }
+    return result;
+}
+
+}  // namespace
+
+TEST_CASE("expm: zero matrix gives identity", "[expm]") {
+    Matrix Z = Matrix::zeros(3, 3);
+    auto E = linalgebra::expm(Z);
+    auto I = Matrix::identity(3);
+    REQUIRE(frobenius_diff(E, I) < 1e-12);
+}
+
+TEST_CASE("expm: scalar multiple of identity", "[expm]") {
+    // expm(s*I) = exp(s)*I
+    const double s = 2.0;
+    Matrix A = Matrix::zeros(3, 3);
+    A(0, 0) = s; A(1, 1) = s; A(2, 2) = s;
+    auto E = linalgebra::expm(A);
+
+    const double expected = std::exp(s);
+    for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+            REQUIRE(E(i, j) == Catch::Approx(i == j ? expected : 0.0).margin(1e-10));
+}
+
+TEST_CASE("expm: 2x2 nilpotent", "[expm]") {
+    // A = [[0,1],[0,0]], expm(A) = [[1,1],[0,1]] exactly.
+    Matrix A{{0.0, 1.0}, {0.0, 0.0}};
+    auto E = linalgebra::expm(A);
+    REQUIRE(E(0, 0) == Catch::Approx(1.0).margin(1e-12));
+    REQUIRE(E(0, 1) == Catch::Approx(1.0).margin(1e-12));
+    REQUIRE(E(1, 0) == Catch::Approx(0.0).margin(1e-12));
+    REQUIRE(E(1, 1) == Catch::Approx(1.0).margin(1e-12));
+}
+
+TEST_CASE("expm: 2x2 rotation generator", "[expm]") {
+    // A = [[0,-t],[t,0]], expm(A) = [[cos(t), -sin(t)],[sin(t), cos(t)]].
+    const double t = std::numbers::pi / 4.0;
+    Matrix A{{0.0, -t}, {t, 0.0}};
+    auto E = linalgebra::expm(A);
+
+    REQUIRE(E(0, 0) == Catch::Approx(std::cos(t)).epsilon(1e-10));
+    REQUIRE(E(0, 1) == Catch::Approx(-std::sin(t)).epsilon(1e-10));
+    REQUIRE(E(1, 0) == Catch::Approx(std::sin(t)).epsilon(1e-10));
+    REQUIRE(E(1, 1) == Catch::Approx(std::cos(t)).epsilon(1e-10));
+}
+
+TEST_CASE("expm: diagonal matrix", "[expm]") {
+    // A = diag(1, 2), expm(A) = diag(e, e^2).
+    Matrix A{{1.0, 0.0}, {0.0, 2.0}};
+    auto E = linalgebra::expm(A);
+    REQUIRE(E(0, 0) == Catch::Approx(std::exp(1.0)).epsilon(1e-10));
+    REQUIRE(E(1, 1) == Catch::Approx(std::exp(2.0)).epsilon(1e-10));
+    REQUIRE(E(0, 1) == Catch::Approx(0.0).margin(1e-12));
+    REQUIRE(E(1, 0) == Catch::Approx(0.0).margin(1e-12));
+}
+
+TEST_CASE("expm: comparison with Taylor series (small-norm A)", "[expm]") {
+    // For small ||A||, expm(A) ≈ Taylor series to order 5.
+    std::mt19937 rng(42);
+    std::uniform_real_distribution<double> dist(-0.01, 0.01);
+    Matrix A(4, 4);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            A(i, j) = dist(rng);
+
+    auto E = linalgebra::expm(A);
+    auto E_taylor = expm_taylor5(A);
+    REQUIRE(frobenius_diff(E, E_taylor) < 1e-10);
+}
+
+TEST_CASE("expm: large norm requires scaling (10*I)", "[expm]") {
+    // expm(10*I) = exp(10)*I.
+    Matrix A = Matrix::zeros(3, 3);
+    A(0, 0) = 10.0; A(1, 1) = 10.0; A(2, 2) = 10.0;
+    auto E = linalgebra::expm(A);
+
+    const double e10 = std::exp(10.0);  // ≈ 22026.47
+    for (std::size_t i = 0; i < 3; ++i)
+        REQUIRE(E(i, i) == Catch::Approx(e10).epsilon(1e-8));
+}
+
+TEST_CASE("expm: non-square throws", "[expm]") {
+    Matrix A(2, 3);
+    REQUIRE_THROWS_AS(linalgebra::expm(A), DimensionMismatchError);
+}
+
+TEST_CASE("expm: symmetric A gives symmetric expm(A)", "[expm]") {
+    std::mt19937 rng(7);
+    std::uniform_real_distribution<double> dist(-1.0, 1.0);
+    Matrix B(4, 4);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            B(i, j) = dist(rng);
+    // Make symmetric: A = B + B^T.
+    Matrix A(4, 4, 0.0);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            A(i, j) = B(i, j) + B(j, i);
+
+    auto E = linalgebra::expm(A);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            REQUIRE(E(i, j) == Catch::Approx(E(j, i)).margin(1e-9));
+}
diff --git a/tests/test_iterative.cpp b/tests/test_iterative.cpp
new file mode 100644
index 0000000..64ee85d
--- /dev/null
+++ b/tests/test_iterative.cpp
@@ -0,0 +1,211 @@
+import linalgebra;
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <random>
+
+using linalgebra::DimensionMismatchError;
+using linalgebra::LinAlgError;
+using linalgebra::Matrix;
+using linalgebra::Vector;
+
+namespace {
+
+double frobenius_diff(const Matrix& A, const Matrix& B) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < A.rows(); ++i)
+        for (std::size_t j = 0; j < A.cols(); ++j) {
+            const double d = A(i, j) - B(i, j);
+            s += d * d;
+        }
+    return std::sqrt(s);
+}
+
+Matrix make_dd(std::size_t n, std::mt19937& rng) {
+    std::uniform_real_distribution<double> dist(0.0, 1.0);
+    Matrix A(n, n);
+    for (std::size_t i = 0; i < n; ++i) {
+        double row_sum = 0.0;
+        for (std::size_t j = 0; j < n; ++j) {
+            A(i, j) = dist(rng);
+            if (i != j) row_sum += std::abs(A(i, j));
+        }
+        A(i, i) = row_sum + 1.0;  // strictly diagonally dominant
+    }
+    return A;
+}
+
+Matrix make_spd(std::size_t n, std::mt19937& rng) {
+    std::uniform_real_distribution<double> dist(-1.0, 1.0);
+    Matrix B(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            B(i, j) = dist(rng);
+    Matrix A = linalgebra::transpose(B) * B;
+    for (std::size_t i = 0; i < n; ++i) A(i, i) += static_cast<double>(n);
+    return A;
+}
+
+}  // namespace
+
+// arnoldi
+
+TEST_CASE("arnoldi: orthonormality and AQ=QH relation", "[iterative][arnoldi]") {
+    Matrix A{{2.0, 1.0, 0.0},
+             {1.0, 3.0, 1.0},
+             {0.0, 1.0, 4.0}};
+    Vector b{1.0, 0.0, 0.0};
+    const int k = 2;
+    auto res = linalgebra::arnoldi(A, b, k);
+
+    const std::size_t n = 3;
+    const auto steps = static_cast<std::size_t>(res.steps_taken);
+
+    // Q columns must be orthonormal.
+    for (std::size_t i = 0; i <= steps; ++i) {
+        for (std::size_t j = 0; j <= steps; ++j) {
+            double dot = 0.0;
+            for (std::size_t r = 0; r < n; ++r) dot += res.Q(r, i) * res.Q(r, j);
+            const double expected = (i == j) ? 1.0 : 0.0;
+            REQUIRE(dot == Catch::Approx(expected).margin(1e-10));
+        }
+    }
+
+    for (std::size_t j = 0; j < steps; ++j) {
+        Vector qj(n);
+        for (std::size_t i = 0; i < n; ++i) qj[i] = res.Q(i, j);
+        const Vector Aqj = A * qj;
+
+        for (std::size_t i = 0; i <= steps; ++i) {
+            double qh = 0.0;
+            for (std::size_t r = 0; r < n; ++r) qh += res.Q(r, i) * res.H(i, j);
+        }
+
+
+        // Direct check: ||A*qj - Q*H[:,j]||
+        Vector Hcol(steps + 1);
+        for (std::size_t i = 0; i <= steps; ++i) Hcol[i] = res.H(i, j);
+        double resid = 0.0;
+        for (std::size_t row = 0; row < n; ++row) {
+            double qh_row = 0.0;
+            for (std::size_t i = 0; i <= steps; ++i) qh_row += res.Q(row, i) * Hcol[i];
+            const double d = Aqj[row] - qh_row;
+            resid += d * d;
+        }
+        REQUIRE(std::sqrt(resid) < 1e-10);
+    }
+}
+
+TEST_CASE("arnoldi: breakdown on scaled identity", "[iterative][arnoldi]") {
+    // A = 2*I → Krylov space is one-dimensional.
+    Matrix A(3, 3, 0.0);
+    A(0, 0) = 2.0; A(1, 1) = 2.0; A(2, 2) = 2.0;
+    Vector b{1.0, 0.0, 0.0};
+    auto res = linalgebra::arnoldi(A, b, 3);
+    REQUIRE(res.breakdown == true);
+    REQUIRE(res.steps_taken <= 3);
+}
+
+TEST_CASE("arnoldi: dimension mismatch throws", "[iterative][arnoldi]") {
+    Matrix A(3, 3, 0.0);
+    Vector b(4, 0.0);
+    REQUIRE_THROWS_AS(linalgebra::arnoldi(A, b, 2), DimensionMismatchError);
+}
+
+// solve_cg
+
+TEST_CASE("solve_cg: 2x2 SPD", "[iterative][cg]") {
+    Matrix A{{4.0, 1.0}, {1.0, 3.0}};
+    Vector b{1.0, 2.0};
+    auto res = linalgebra::solve_cg(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-10);
+}
+
+TEST_CASE("solve_cg: identity system converges in 1 iteration", "[iterative][cg]") {
+    auto I = Matrix::identity(5);
+    Vector b{1.0, 2.0, 3.0, 4.0, 5.0};
+    auto res = linalgebra::solve_cg(I, b);
+    for (std::size_t i = 0; i < 5; ++i)
+        REQUIRE(res.x[i] == Catch::Approx(b[i]).margin(1e-10));
+}
+
+TEST_CASE("solve_cg: 5x5 random SPD", "[iterative][cg]") {
+    std::mt19937 rng(99);
+    Matrix A = make_spd(5, rng);
+    std::uniform_real_distribution<double> dist(-2.0, 2.0);
+    Vector b(5);
+    for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng);
+
+    auto res = linalgebra::solve_cg(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8);
+
+    // Cross-check vs LU.
+    auto lu  = linalgebra::lu_factor(A);
+    auto x_lu = linalgebra::lu_solve(lu, b);
+    REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8);
+}
+
+TEST_CASE("solve_cg: non-symmetric throws", "[iterative][cg]") {
+    Matrix A{{1.0, 2.0}, {0.0, 1.0}};
+    Vector b{1.0, 1.0};
+    REQUIRE_THROWS_AS(linalgebra::solve_cg(A, b), LinAlgError);
+}
+
+// solve_gmres
+
+TEST_CASE("solve_gmres: 2x2 non-symmetric", "[iterative][gmres]") {
+    Matrix A{{2.0, 1.0}, {1.0, 3.0}};
+    Vector b{5.0, 7.0};
+    auto res = linalgebra::solve_gmres(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-9);
+}
+
+TEST_CASE("solve_gmres: 5x5 diagonally dominant", "[iterative][gmres]") {
+    std::mt19937 rng(11);
+    Matrix A = make_dd(5, rng);
+    std::uniform_real_distribution<double> dist(-3.0, 3.0);
+    Vector b(5);
+    for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng);
+
+    auto res = linalgebra::solve_gmres(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8);
+
+    auto lu = linalgebra::lu_factor(A);
+    auto x_lu = linalgebra::lu_solve(lu, b);
+    REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8);
+}
+
+TEST_CASE("solve_gmres: zero rhs gives zero solution", "[iterative][gmres]") {
+    Matrix A{{2.0, 1.0}, {1.0, 3.0}};
+    Vector b(2, 0.0);
+    auto res = linalgebra::solve_gmres(A, b);
+    REQUIRE(linalgebra::norm2(res.x) < 1e-12);
+}
+
+
+// solve_bicgstab
+
+TEST_CASE("solve_bicgstab: 2x2 non-symmetric", "[iterative][bicgstab]") {
+    Matrix A{{2.0, 1.0}, {1.0, 3.0}};
+    Vector b{5.0, 7.0};
+    auto res = linalgebra::solve_bicgstab(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-9);
+}
+
+TEST_CASE("solve_bicgstab: 5x5 diagonally dominant", "[iterative][bicgstab]") {
+    std::mt19937 rng(22);
+    Matrix A = make_dd(5, rng);
+    std::uniform_real_distribution<double> dist(-3.0, 3.0);
+    Vector b(5);
+    for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng);
+
+    auto res = linalgebra::solve_bicgstab(A, b);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8);
+
+    auto lu = linalgebra::lu_factor(A);
+    auto x_lu = linalgebra::lu_solve(lu, b);
+    REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8);
+}
diff --git a/tests/test_precond.cpp b/tests/test_precond.cpp
new file mode 100644
index 0000000..160badf
--- /dev/null
+++ b/tests/test_precond.cpp
@@ -0,0 +1,142 @@
+import linalgebra;
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <cstddef>
+
+using linalgebra::DimensionMismatchError;
+using linalgebra::Matrix;
+using linalgebra::SingularMatrixError;
+using linalgebra::Vector;
+
+// condition_number_1norm
+
+TEST_CASE("condition_number_1norm: identity", "[precond][condition]") {
+    for (std::size_t n : {1, 2, 5}) {
+        auto I = Matrix::identity(n);
+        REQUIRE(linalgebra::condition_number_1norm(I) == Catch::Approx(1.0).epsilon(1e-10));
+    }
+}
+
+TEST_CASE("condition_number_1norm: diagonal matrix", "[precond][condition]") {
+    // diag(1, 10, 100): ||A||_1 = 100, ||A^{-1}||_1 = 1, so cond = 100.
+    Matrix A(3, 3, 0.0);
+    A(0, 0) = 1.0; A(1, 1) = 10.0; A(2, 2) = 100.0;
+    const double c = linalgebra::condition_number_1norm(A);
+    // Power-iteration estimator gives a lower bound; for simple diagonal it should be exact.
+    REQUIRE(c == Catch::Approx(100.0).epsilon(1e-8));
+}
+
+TEST_CASE("condition_number_1norm: ill-conditioned Hilbert 4x4", "[precond][condition]") {
+    // Hilbert matrix H(i,j) = 1/(i+j+1).
+    Matrix H(4, 4);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            H(i, j) = 1.0 / static_cast<double>(i + j + 1);
+    const double c = linalgebra::condition_number_1norm(H);
+    REQUIRE(c > 1e3);  // Hilbert matrices are notoriously ill-conditioned
+}
+
+TEST_CASE("condition_number_1norm: non-square throws", "[precond][condition]") {
+    Matrix A(2, 3);
+    REQUIRE_THROWS_AS(linalgebra::condition_number_1norm(A), DimensionMismatchError);
+}
+
+// precond_jacobi
+
+TEST_CASE("precond_jacobi: diagonal matrix", "[precond][jacobi]") {
+    Matrix A(3, 3, 0.0);
+    A(0, 0) = 2.0; A(1, 1) = 4.0; A(2, 2) = 8.0;
+    auto P = linalgebra::precond_jacobi(A);
+
+    REQUIRE(P.inv_diag[0] == Catch::Approx(0.5));
+    REQUIRE(P.inv_diag[1] == Catch::Approx(0.25));
+    REQUIRE(P.inv_diag[2] == Catch::Approx(0.125));
+
+    Vector x{1.0, 1.0, 1.0};
+    auto y = linalgebra::apply(P, x);
+    REQUIRE(y[0] == Catch::Approx(0.5));
+    REQUIRE(y[1] == Catch::Approx(0.25));
+    REQUIRE(y[2] == Catch::Approx(0.125));
+}
+
+TEST_CASE("precond_jacobi: zero diagonal throws", "[precond][jacobi]") {
+    Matrix A{{1.0, 0.0}, {0.0, 0.0}};
+    REQUIRE_THROWS_AS(linalgebra::precond_jacobi(A), SingularMatrixError);
+}
+
+TEST_CASE("precond_jacobi: non-square throws", "[precond][jacobi]") {
+    Matrix A(2, 3);
+    REQUIRE_THROWS_AS(linalgebra::precond_jacobi(A), DimensionMismatchError);
+}
+
+// precond_ilu0
+
+TEST_CASE("precond_ilu0: identity", "[precond][ilu0]") {
+    auto I = Matrix::identity(3);
+    auto P = linalgebra::precond_ilu0(I);
+    Vector b{3.0, 1.0, 4.0};
+    auto y = linalgebra::apply(P, b);
+    for (std::size_t i = 0; i < 3; ++i)
+        REQUIRE(y[i] == Catch::Approx(b[i]).margin(1e-12));
+}
+
+TEST_CASE("precond_ilu0: 3x3 SPD", "[precond][ilu0]") {
+    // Dense ILU0 = exact LU without pivoting, so for any nonsingular A,
+    // apply(P, A*x) should recover x.
+    Matrix A{{4.0, 2.0, 0.0},
+             {2.0, 3.0, 1.0},
+             {0.0, 1.0, 2.0}};
+    auto P = linalgebra::precond_ilu0(A);
+
+    Vector x_true{1.0, 2.0, 3.0};
+    Vector rhs = A * x_true;
+    auto x_rec = linalgebra::apply(P, rhs);
+
+    for (std::size_t i = 0; i < 3; ++i)
+        REQUIRE(x_rec[i] == Catch::Approx(x_true[i]).margin(1e-10));
+}
+
+TEST_CASE("precond_ilu0: near-singular pivot throws", "[precond][ilu0]") {
+    Matrix A{{0.0, 1.0}, {1.0, 1.0}};
+    REQUIRE_THROWS_AS(linalgebra::precond_ilu0(A), SingularMatrixError);
+}
+
+// lstsq
+
+TEST_CASE("lstsq: square full-rank system", "[precond][lstsq]") {
+    Matrix A{{2.0, 1.0}, {1.0, 3.0}};
+    Vector b{5.0, 7.0};
+    auto res = linalgebra::lstsq(A, b);
+    REQUIRE(res.rank == 2);
+    REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-10);
+}
+
+TEST_CASE("lstsq: overdetermined full-rank", "[precond][lstsq]") {
+    // A is 3x2, consistent overdetermined system.
+    Matrix A{{1.0, 1.0}, {2.0, 1.0}, {3.0, 1.0}};
+    Vector b{6.0, 5.0, 7.0};
+    auto res = linalgebra::lstsq(A, b);
+    REQUIRE(res.rank == 2);
+    // Residual should be the minimum achievable (verify normal equations: A^T A x = A^T b).
+    Vector AtAx = linalgebra::transpose(A) * (A * res.x);
+    Vector Atb  = linalgebra::transpose(A) * b;
+    for (std::size_t i = 0; i < 2; ++i)
+        REQUIRE(AtAx[i] == Catch::Approx(Atb[i]).margin(1e-8));
+}
+
+TEST_CASE("lstsq: rank-deficient", "[precond][lstsq]") {
+    // col2 = 2*col1 → rank 1.
+    Matrix A{{1.0, 2.0}, {2.0, 4.0}, {3.0, 6.0}};
+    Vector b{1.0, 2.0, 3.0};
+    auto res = linalgebra::lstsq(A, b);
+    REQUIRE(res.rank < 2);
+    // The residual should be the minimum achievable.
+    REQUIRE(res.residual_norm < 1e-8);
+}
+
+TEST_CASE("lstsq: non-tall matrix throws", "[precond][lstsq]") {
+    Matrix A(2, 3);
+    Vector b(2, 0.0);
+    REQUIRE_THROWS_AS(linalgebra::lstsq(A, b), DimensionMismatchError);
+}
diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp
index d78fbb1..cbde901 100644
--- a/tests/test_qr.cpp
+++ b/tests/test_qr.cpp
@@ -224,9 +224,7 @@ TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") {
     CHECK_NOTHROW(linalgebra::qr_householder(A));
 }
 
-// ---------------------------------------------------------------------------
 // qr_colpiv tests
-// ---------------------------------------------------------------------------
 
 TEST_CASE("QR ColPiv: full rank reconstruction", "[qr][colpiv]") {
     const Matrix A = random_matrix(6, 4, 100u);
diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp
index 64a5dd6..95bb044 100644
--- a/tests/test_qr_iteration.cpp
+++ b/tests/test_qr_iteration.cpp
@@ -423,9 +423,7 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n",
     std::cout << "=============================================\n";
 }
 
-// ---------------------------------------------------------------------------
 // Francis double-shift QR tests
-// ---------------------------------------------------------------------------
 
 TEST_CASE("Francis QR: 2x2 real eigenvalues", "[qr_iteration][francis]") {
     Matrix A{{3.0, 1.0}, {0.0, 2.0}};
diff --git a/tests/test_svd.cpp b/tests/test_svd.cpp
new file mode 100644
index 0000000..d2121d6
--- /dev/null
+++ b/tests/test_svd.cpp
@@ -0,0 +1,140 @@
+import linalgebra;
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <random>
+
+using linalgebra::DimensionMismatchError;
+using linalgebra::Matrix;
+using linalgebra::SVDResult;
+
+namespace {
+
+double frobenius_diff(const Matrix& A, const Matrix& B) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < A.rows(); ++i)
+        for (std::size_t j = 0; j < A.cols(); ++j) {
+            const double d = A(i, j) - B(i, j);
+            s += d * d;
+        }
+    return std::sqrt(s);
+}
+
+// Reconstruct A = U * diag(sigma) * Vt and measure error
+double reconstruction_error(const Matrix& A, const SVDResult& svd) {
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+    const std::size_t p = svd.sigma.size();
+
+    // Build U_thin (m x p) — first p columns of U.
+    Matrix U_thin(m, p);
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < p; ++j)
+            U_thin(i, j) = svd.U(i, j);
+
+    // Build Sigma_diag (p x p).
+    Matrix S(p, p, 0.0);
+    for (std::size_t i = 0; i < p; ++i) S(i, i) = svd.sigma[i];
+
+    // Build Vt_thin (p x n) — first p rows of Vt.
+    Matrix Vt_thin(p, n);
+    for (std::size_t i = 0; i < p; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            Vt_thin(i, j) = svd.Vt(i, j);
+
+    const Matrix recon = U_thin * S * Vt_thin;
+    return frobenius_diff(A, recon);
+}
+
+// Check orthogonality of the first k columns of M
+double col_ortho_error(const Matrix& M, std::size_t k) {
+    const std::size_t m = M.rows();
+    double s = 0.0;
+    for (std::size_t i = 0; i < k; ++i) {
+        for (std::size_t j = 0; j < k; ++j) {
+            double dot = 0.0;
+            for (std::size_t r = 0; r < m; ++r) dot += M(r, i) * M(r, j);
+            const double expected = (i == j) ? 1.0 : 0.0;
+            const double d = dot - expected;
+            s += d * d;
+        }
+    }
+    return std::sqrt(s);
+}
+
+}  // namespace
+
+TEST_CASE("svd: 2x2 diagonal", "[svd]") {
+    Matrix A{{3.0, 0.0}, {0.0, -2.0}};
+    auto res = linalgebra::svd(A);
+
+    REQUIRE(res.sigma[0] == Catch::Approx(3.0).epsilon(1e-10));
+    REQUIRE(res.sigma[1] == Catch::Approx(2.0).epsilon(1e-10));
+    REQUIRE(reconstruction_error(A, res) < 1e-10);
+}
+
+TEST_CASE("svd: 3x2 tall matrix", "[svd]") {
+    Matrix A{{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
+    auto res = linalgebra::svd(A);
+
+    REQUIRE(reconstruction_error(A, res) < 1e-10);
+
+    REQUIRE(res.sigma[0] >= res.sigma[1]);
+    REQUIRE(res.sigma[1] >= 0.0);
+
+    REQUIRE(col_ortho_error(res.U, 2) < 1e-10);
+
+    REQUIRE(col_ortho_error(linalgebra::transpose(res.Vt), 2) < 1e-10);
+}
+
+TEST_CASE("svd: rank-1 matrix", "[svd]") {
+    // A = u * v^T for u = [1,2,3], v = [2,1].
+    Matrix A{{2.0, 1.0}, {4.0, 2.0}, {6.0, 3.0}};
+    auto res = linalgebra::svd(A);
+
+    REQUIRE(reconstruction_error(A, res) < 1e-10);
+    // Second singular value should be near zero.
+    REQUIRE(res.sigma[0] > 1e-10);
+    REQUIRE(res.sigma[1] < 1e-8);
+}
+
+TEST_CASE("svd: identity 4x4", "[svd]") {
+    auto I = Matrix::identity(4);
+    auto res = linalgebra::svd(I);
+
+    for (std::size_t i = 0; i < 4; ++i)
+        REQUIRE(res.sigma[i] == Catch::Approx(1.0).epsilon(1e-10));
+    REQUIRE(reconstruction_error(I, res) < 1e-10);
+}
+
+TEST_CASE("svd: 4x4 random", "[svd]") {
+    std::mt19937 rng(42);
+    std::uniform_real_distribution<double> dist(-2.0, 2.0);
+    Matrix A(4, 4);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            A(i, j) = dist(rng);
+
+    auto res = linalgebra::svd(A);
+    REQUIRE(reconstruction_error(A, res) < 1e-9);
+    REQUIRE(col_ortho_error(res.U, 4) < 1e-10);
+    REQUIRE(col_ortho_error(linalgebra::transpose(res.Vt), 4) < 1e-10);
+}
+
+TEST_CASE("svd: Hilbert 4x4", "[svd]") {
+    Matrix H(4, 4);
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+            H(i, j) = 1.0 / static_cast<double>(i + j + 1);
+    auto res = linalgebra::svd(H);
+    REQUIRE(reconstruction_error(H, res) < 1e-10);
+    for (std::size_t i = 0; i < 4; ++i) REQUIRE(res.sigma[i] >= 0.0);
+}
+
+TEST_CASE("svd: wide matrix throws", "[svd]") {
+    Matrix A(2, 4);
+    REQUIRE_THROWS_AS(linalgebra::svd(A), DimensionMismatchError);
+}
diff --git a/tests/test_sym_eigen.cpp b/tests/test_sym_eigen.cpp
new file mode 100644
index 0000000..3906785
--- /dev/null
+++ b/tests/test_sym_eigen.cpp
@@ -0,0 +1,168 @@
+import linalgebra;
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <random>
+
+using linalgebra::DimensionMismatchError;
+using linalgebra::LinAlgError;
+using linalgebra::Matrix;
+using linalgebra::Vector;
+
+namespace {
+
+double frobenius_diff(const Matrix& A, const Matrix& B) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < A.rows(); ++i)
+        for (std::size_t j = 0; j < A.cols(); ++j) {
+            const double d = A(i, j) - B(i, j);
+            s += d * d;
+        }
+    return std::sqrt(s);
+}
+
+// Check orthogonality ||Q^T Q - I||_F < tol.
+double ortho_error(const Matrix& Q) {
+    const std::size_t n = Q.rows();
+    const std::size_t m = Q.cols();
+    const Matrix Qt = linalgebra::transpose(Q);
+    const Matrix QtQ = Qt * Q;
+    const Matrix I = Matrix::identity(m);
+    return frobenius_diff(QtQ, I);
+}
+
+// Build a symmetric matrix from B^T B + shift * I.
+Matrix make_sym(std::size_t n, double shift, std::mt19937& rng) {
+    std::uniform_real_distribution<double> dist(-1.0, 1.0);
+    Matrix B(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            B(i, j) = dist(rng);
+    Matrix A = linalgebra::transpose(B) * B;
+    for (std::size_t i = 0; i < n; ++i) A(i, i) += shift;
+    return A;
+}
+
+}  // namespace
+
+// tridiagonalize
+
+TEST_CASE("tridiagonalize: 2x2 symmetric", "[sym_eigen][tridiagonalize]") {
+    Matrix A{{4.0, 2.0}, {2.0, 3.0}};
+    auto res = linalgebra::tridiagonalize(A);
+
+    // Reconstruction: Q T Q^T == A.
+    const Matrix QTQT = res.Q * res.T * linalgebra::transpose(res.Q);
+    REQUIRE(frobenius_diff(A, QTQT) < 1e-12);
+
+    // Q must be orthogonal.
+    REQUIRE(ortho_error(res.Q) < 1e-12);
+
+    // T must be tridiagonal: T(i,j) == 0 for |i-j| > 1.
+    const std::size_t n = res.T.rows();
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            if (i > j + 1 || j > i + 1)
+                REQUIRE(std::abs(res.T(i, j)) < 1e-12);
+}
+
+TEST_CASE("tridiagonalize: 4x4 symmetric", "[sym_eigen][tridiagonalize]") {
+    Matrix A{{6.0, 2.0, 1.0, 0.0},
+             {2.0, 5.0, 3.0, 1.0},
+             {1.0, 3.0, 4.0, 2.0},
+             {0.0, 1.0, 2.0, 3.0}};
+    auto res = linalgebra::tridiagonalize(A);
+
+    const Matrix QTQT = res.Q * res.T * linalgebra::transpose(res.Q);
+    REQUIRE(frobenius_diff(A, QTQT) < 1e-10);
+    REQUIRE(ortho_error(res.Q) < 1e-12);
+
+    const std::size_t n = res.T.rows();
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            if (i > j + 1 || j > i + 1)
+                REQUIRE(std::abs(res.T(i, j)) < 1e-10);
+}
+
+TEST_CASE("tridiagonalize: identity", "[sym_eigen][tridiagonalize]") {
+    auto I = Matrix::identity(5);
+    auto res = linalgebra::tridiagonalize(I);
+    REQUIRE(frobenius_diff(I, res.Q * res.T * linalgebra::transpose(res.Q)) < 1e-12);
+}
+
+TEST_CASE("tridiagonalize: random symmetric", "[sym_eigen][tridiagonalize]") {
+    std::mt19937 rng(77);
+    Matrix A = make_sym(8, 5.0, rng);
+    auto res = linalgebra::tridiagonalize(A);
+    REQUIRE(frobenius_diff(A, res.Q * res.T * linalgebra::transpose(res.Q)) < 1e-9);
+    REQUIRE(ortho_error(res.Q) < 1e-11);
+}
+
+TEST_CASE("tridiagonalize: non-symmetric throws", "[sym_eigen][tridiagonalize]") {
+    Matrix A{{1.0, 2.0}, {3.0, 4.0}};
+    REQUIRE_THROWS_AS(linalgebra::tridiagonalize(A), LinAlgError);
+}
+
+TEST_CASE("tridiagonalize: non-square throws", "[sym_eigen][tridiagonalize]") {
+    Matrix A(2, 3);
+    REQUIRE_THROWS_AS(linalgebra::tridiagonalize(A), DimensionMismatchError);
+}
+
+// eigenvectors_inverse_iteration
+
+TEST_CASE("eigenvectors_inverse_iteration: 2x2 diagonal", "[sym_eigen][inverse_iter]") {
+    Matrix A(2, 2, 0.0);
+    A(0, 0) = 3.0; A(1, 1) = 7.0;
+    Vector lambdas{3.0, 7.0};
+
+    auto res = linalgebra::eigenvectors_inverse_iteration(A, lambdas);
+
+    // Each column should satisfy A*v ≈ lambda*v.
+    for (std::size_t col = 0; col < 2; ++col) {
+        Vector v(2);
+        v[0] = res.eigenvectors(0, col);
+        v[1] = res.eigenvectors(1, col);
+        const Vector Av = A * v;
+        const double lam = lambdas[col];
+        double resid = 0.0;
+        for (std::size_t i = 0; i < 2; ++i) {
+            const double d = Av[i] - lam * v[i];
+            resid += d * d;
+        }
+        REQUIRE(std::sqrt(resid) < 1e-8);
+    }
+}
+
+TEST_CASE("eigenvectors_inverse_iteration: 3x3 symmetric", "[sym_eigen][inverse_iter]") {
+    // Known symmetric matrix: compute eigenvalues with francis, then eigenvectors.
+    Matrix A{{6.0, 2.0, 1.0},
+             {2.0, 3.0, 1.0},
+             {1.0, 1.0, 1.0}};
+
+    auto eig = linalgebra::eigenvalues_francis(A);
+    // Use real eigenvalues only.
+    auto res = linalgebra::eigenvectors_inverse_iteration(A, eig.eigenvalues_real);
+
+    for (std::size_t col = 0; col < 3; ++col) {
+        Vector v(3);
+        for (std::size_t i = 0; i < 3; ++i) v[i] = res.eigenvectors(i, col);
+        const Vector Av = A * v;
+        const double lam = eig.eigenvalues_real[col];
+        double resid = 0.0;
+        for (std::size_t i = 0; i < 3; ++i) {
+            const double d = Av[i] - lam * v[i];
+            resid += d * d;
+        }
+        REQUIRE(std::sqrt(resid) < 1e-6);
+    }
+}
+
+TEST_CASE("eigenvectors_inverse_iteration: non-square throws", "[sym_eigen][inverse_iter]") {
+    Matrix A(2, 3);
+    Vector lambdas{1.0};
+    REQUIRE_THROWS_AS(linalgebra::eigenvectors_inverse_iteration(A, lambdas),
+                      DimensionMismatchError);
+}