From 4602b36e9d5ea08656e3222846a1f161bbb1cec1 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Mon, 27 Apr 2026 07:24:41 +0300 Subject: Module refactor --- tests/test_lu.cpp | 109 +++++++++++++++--------------------------------------- 1 file changed, 29 insertions(+), 80 deletions(-) (limited to 'tests/test_lu.cpp') diff --git a/tests/test_lu.cpp b/tests/test_lu.cpp index dfb76c4..c4f9c5d 100644 --- a/tests/test_lu.cpp +++ b/tests/test_lu.cpp @@ -1,8 +1,4 @@ -#include "linalg_error.hpp" -#include "lu.hpp" -#include "matrix.hpp" -#include "norms.hpp" -#include "vector.hpp" +import linalgebra; #include #include @@ -11,31 +7,23 @@ #include #include -using linalg::DimensionMismatchError; -using linalg::Matrix; -using linalg::SingularMatrixError; -using linalg::Vector; -using linalg::LUResult; - -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::SingularMatrixError; +using linalgebra::Vector; +using linalgebra::LUResult; namespace { -// ||PA - LU||_F (Frobenius, computed element-wise via ||vec||_2) double reconstruction_error(const Matrix& A, const LUResult& lu) { const std::size_t n = A.rows(); - // Build PA by permuting rows of A. Matrix PA(n, n); for (std::size_t i = 0; i < n; ++i) { for (std::size_t j = 0; j < n; ++j) { PA(i, j) = A(lu.perm[i], j); } } - // Compute LU product. const Matrix LU = lu.L * lu.U; - // Compute Frobenius norm of (PA - LU). double err = 0.0; for (std::size_t i = 0; i < n; ++i) { for (std::size_t j = 0; j < n; ++j) { @@ -46,12 +34,10 @@ double reconstruction_error(const Matrix& A, const LUResult& lu) { return std::sqrt(err); } -// ||Ax - b||_2 double solve_residual(const Matrix& A, const Vector& x, const Vector& b) { - return linalg::norm2(A * x - b); + return linalgebra::norm2(A * x - b); } -// Generate a reproducible random nonsingular n x n matrix. Matrix random_matrix(std::size_t n, unsigned seed = 42) { std::mt19937 rng(seed); std::uniform_real_distribution dist(-10.0, 10.0); @@ -66,17 +52,13 @@ Matrix random_matrix(std::size_t n, unsigned seed = 42) { } // namespace -// --------------------------------------------------------------------------- -// Factorization correctness -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: 3x3 known system", "[lu]") { const Matrix A{ {2.0, 1.0, -1.0}, {-3.0, -1.0, 2.0}, {-2.0, 1.0, 2.0} }; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); REQUIRE(lu.L.rows() == 3); REQUIRE(lu.U.rows() == 3); @@ -86,41 +68,33 @@ TEST_CASE("LU factorization: 3x3 known system", "[lu]") { CHECK(lu.L(i, i) == Catch::Approx(1.0)); } - // ||PA - LU|| must be near zero. CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-12)); } TEST_CASE("LU factorization: identity matrix", "[lu]") { const Matrix I = Matrix::identity(4); - const LUResult lu = linalg::lu_factor(I); + const LUResult lu = linalgebra::lu_factor(I); CHECK(reconstruction_error(I, lu) == Catch::Approx(0.0).margin(1e-14)); - // U should equal I (up to row ordering already handled by PA=LU). for (std::size_t i = 0; i < 4; ++i) { CHECK(lu.U(i, i) == Catch::Approx(1.0)); } } TEST_CASE("LU factorization: matrix requiring row swaps", "[lu]") { - // First column entry is zero. No-pivot LU would immediately fail. const Matrix A{ {0.0, 1.0, 2.0}, {3.0, 4.0, 5.0}, {6.0, 7.0, 8.0} }; - // A is singular (rows are in AP), but check that partial pivoting still - // proceeds and detects singularity correctly. - // Row 3 - row 2 = row 2 - row 1, so rank < 3. - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: first-column zero, nonsingular", "[lu]") { - // [[0, 1], [1, 0]] — requires a swap at step 0. const Matrix A{{0.0, 1.0}, {1.0, 0.0}}; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-14)); - // Solving Ax = b: A swaps components. const Vector b{3.0, 7.0}; - const Vector x = linalg::lu_solve(lu, b); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-12)); CHECK(x[0] == Catch::Approx(7.0)); CHECK(x[1] == Catch::Approx(3.0)); @@ -129,17 +103,12 @@ TEST_CASE("LU factorization: first-column zero, nonsingular", "[lu]") { TEST_CASE("LU factorization: random nonsingular matrices", "[lu]") { for (std::size_t n : {5u, 10u, 20u}) { const Matrix A = random_matrix(n, 123u + static_cast(n)); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-10)); } } -// --------------------------------------------------------------------------- -// Solve correctness -// --------------------------------------------------------------------------- - TEST_CASE("LU solve: known 3x3 system", "[lu]") { - // From Cramer / textbook: solution is x = (2, 3, -1). const Matrix A{ {2.0, 1.0, -1.0}, {-3.0, -1.0, 2.0}, @@ -148,8 +117,8 @@ TEST_CASE("LU solve: known 3x3 system", "[lu]") { const Vector expected{2.0, 3.0, -1.0}; const Vector b = A * expected; - const LUResult lu = linalg::lu_factor(A); - const Vector x = linalg::lu_solve(lu, b); + const LUResult lu = linalgebra::lu_factor(A); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(x[0] == Catch::Approx(expected[0]).epsilon(1e-12)); CHECK(x[1] == Catch::Approx(expected[1]).epsilon(1e-12)); @@ -160,92 +129,77 @@ TEST_CASE("LU solve: known 3x3 system", "[lu]") { TEST_CASE("LU solve: random nonsingular systems", "[lu]") { for (std::size_t n : {5u, 15u, 30u}) { const Matrix A = random_matrix(n, 7u * static_cast(n)); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); - // Random rhs. std::mt19937 rng(n); std::uniform_real_distribution dist(-5.0, 5.0); Vector b(n); for (std::size_t i = 0; i < n; ++i) b[i] = dist(rng); - const Vector x = linalg::lu_solve(lu, b); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-9)); } } TEST_CASE("LU solve: diagonal system", "[lu]") { - // D = diag(2, 3, 4), b = (2, 9, 8), solution = (1, 3, 2). const Matrix D{ {2.0, 0.0, 0.0}, {0.0, 3.0, 0.0}, {0.0, 0.0, 4.0} }; const Vector b{2.0, 9.0, 8.0}; - const LUResult lu = linalg::lu_factor(D); - const Vector x = linalg::lu_solve(lu, b); + const LUResult lu = linalgebra::lu_factor(D); + const Vector x = linalgebra::lu_solve(lu, b); CHECK(x[0] == Catch::Approx(1.0)); CHECK(x[1] == Catch::Approx(3.0)); CHECK(x[2] == Catch::Approx(2.0)); } -// --------------------------------------------------------------------------- -// Failure cases -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: non-square matrix throws", "[lu]") { const Matrix A(3, 4); - CHECK_THROWS_AS(linalg::lu_factor(A), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), DimensionMismatchError); } TEST_CASE("LU factorization: exactly singular matrix throws", "[lu]") { - // Zero row → singular. const Matrix A{ {1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {0.0, 0.0, 0.0} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: rank-deficient matrix throws", "[lu]") { - // Row 2 is a linear combination of rows 0 and 1. const Matrix A{ {1.0, 2.0}, {2.0, 4.0} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU factorization: near-singular matrix throws at default tolerance", "[lu]") { - // Pivot reduced to ~1e-16, should trip the singularity check. const Matrix A{ {1.0, 1.0}, {1.0, 1.0 + 1e-16} }; - CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError); + CHECK_THROWS_AS(linalgebra::lu_factor(A), SingularMatrixError); } TEST_CASE("LU solve: mismatched rhs throws", "[lu]") { const Matrix A = Matrix::identity(3); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const Vector b(5, 1.0); - CHECK_THROWS_AS(linalg::lu_solve(lu, b), DimensionMismatchError); + CHECK_THROWS_AS(linalgebra::lu_solve(lu, b), DimensionMismatchError); } -// --------------------------------------------------------------------------- -// L and U structure -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: L is unit lower triangular", "[lu]") { const Matrix A = random_matrix(6, 999u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const std::size_t n = A.rows(); for (std::size_t i = 0; i < n; ++i) { - // Unit diagonal. CHECK(lu.L(i, i) == Catch::Approx(1.0)); - // Strict upper triangle is zero. for (std::size_t j = i + 1; j < n; ++j) { CHECK(lu.L(i, j) == Catch::Approx(0.0).margin(1e-15)); } @@ -254,7 +208,7 @@ TEST_CASE("LU factorization: L is unit lower triangular", "[lu]") { TEST_CASE("LU factorization: U is upper triangular", "[lu]") { const Matrix A = random_matrix(6, 777u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); const std::size_t n = A.rows(); for (std::size_t i = 1; i < n; ++i) { @@ -264,20 +218,15 @@ TEST_CASE("LU factorization: U is upper triangular", "[lu]") { } } -// --------------------------------------------------------------------------- -// Permutation sign and determinant -// --------------------------------------------------------------------------- - TEST_CASE("LU factorization: sign of permutation is ±1", "[lu]") { const Matrix A = random_matrix(5, 321u); - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); CHECK((lu.sign == 1 || lu.sign == -1)); } TEST_CASE("LU factorization: determinant via sign * prod(diag(U))", "[lu]") { - // det([[3,1],[2,4]]) = 12 - 2 = 10 const Matrix A{{3.0, 1.0}, {2.0, 4.0}}; - const LUResult lu = linalg::lu_factor(A); + const LUResult lu = linalgebra::lu_factor(A); double det = static_cast(lu.sign); for (std::size_t i = 0; i < A.rows(); ++i) { det *= lu.U(i, i); -- cgit v1.2.3