aboutsummaryrefslogtreecommitdiff
path: root/tests
diff options
context:
space:
mode:
Diffstat (limited to 'tests')
-rw-r--r--tests/test_lu.cpp109
-rw-r--r--tests/test_matrix.cpp14
-rw-r--r--tests/test_qr.cpp131
-rw-r--r--tests/test_qr_iteration.cpp200
-rw-r--r--tests/test_triangular_solve.cpp31
-rw-r--r--tests/test_vector.cpp16
6 files changed, 125 insertions, 376 deletions
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 <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -11,31 +7,23 @@
#include <cstddef>
#include <random>
-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<double> 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<unsigned>(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<unsigned>(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<double> 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<double>(lu.sign);
for (std::size_t i = 0; i < A.rows(); ++i) {
det *= lu.U(i, i);
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
index fcd3d27..2fd9a5e 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -1,15 +1,11 @@
-#include "linalg_error.hpp"
-#include "matrix.hpp"
+import linalgebra;
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-#include <type_traits>
-#include <utility>
-
-using linalg::Matrix;
-using linalg::Vector;
-using linalg::DimensionMismatchError;
+using linalgebra::Matrix;
+using linalgebra::Vector;
+using linalgebra::DimensionMismatchError;
TEST_CASE("Matrix constructors initialize dimensions and values", "[matrix]") {
const Matrix empty;
@@ -101,7 +97,7 @@ TEST_CASE("Matrix initializer list rejects unequal row lengths", "[matrix]") {
TEST_CASE("Matrix transpose swaps rows and columns", "[matrix]") {
const Matrix a{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}};
- const Matrix at = linalg::transpose(a);
+ const Matrix at = linalgebra::transpose(a);
REQUIRE(at.rows() == 3);
REQUIRE(at.cols() == 2);
diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp
index a837132..1a4114f 100644
--- a/tests/test_qr.cpp
+++ b/tests/test_qr.cpp
@@ -1,8 +1,4 @@
-#include "linalg_error.hpp"
-#include "matrix.hpp"
-#include "norms.hpp"
-#include "qr.hpp"
-#include "vector.hpp"
+import linalgebra;
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -12,18 +8,13 @@
#include <functional>
#include <random>
-using linalg::DimensionMismatchError;
-using linalg::Matrix;
-using linalg::QRResult;
-using linalg::SingularMatrixError;
-
-// ---------------------------------------------------------------------------
-// Helpers
-// ---------------------------------------------------------------------------
+using linalgebra::DimensionMismatchError;
+using linalgebra::Matrix;
+using linalgebra::QRResult;
+using linalgebra::SingularMatrixError;
namespace {
-// ||A - QR||_F
double reconstruction_error(const Matrix& A, const QRResult& qr) {
const Matrix diff = A - qr.Q * qr.R;
double err = 0.0;
@@ -33,11 +24,9 @@ double reconstruction_error(const Matrix& A, const QRResult& qr) {
return std::sqrt(err);
}
-// ||Q^T Q - I||_F (should be ~0 for orthonormal Q)
double orthogonality_error(const QRResult& qr) {
const Matrix& Q = qr.Q;
const std::size_t n = Q.cols();
- // Compute Q^T Q
Matrix QtQ(n, n);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j) {
@@ -45,7 +34,6 @@ double orthogonality_error(const QRResult& qr) {
for (std::size_t k = 0; k < Q.rows(); ++k) s += Q(k, i) * Q(k, j);
QtQ(i, j) = s;
}
- // ||QtQ - I||_F
double err = 0.0;
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j) {
@@ -55,7 +43,6 @@ double orthogonality_error(const QRResult& qr) {
return std::sqrt(err);
}
-// R must be upper triangular (strict lower triangle near zero).
bool r_is_upper_triangular(const Matrix& R, double tol = 1e-12) {
for (std::size_t i = 1; i < R.rows(); ++i)
for (std::size_t j = 0; j < i; ++j)
@@ -73,7 +60,6 @@ Matrix random_matrix(std::size_t m, std::size_t n, unsigned seed = 42) {
return M;
}
-// Run all checks for a given QR function and matrix.
using QRFn = std::function<QRResult(const Matrix&)>;
void check_qr(const Matrix& A, QRFn fn,
@@ -91,33 +77,25 @@ void check_qr(const Matrix& A, QRFn fn,
} // namespace
-// ---------------------------------------------------------------------------
-// Macro to run the same test body for all three methods
-// ---------------------------------------------------------------------------
-
-#define FOR_ALL_METHODS(A, recon_tol, ortho_tol) \
- SECTION("classical_gs") { \
- check_qr(A, [](const Matrix& M) { return linalg::qr_classical_gs(M); }, \
- recon_tol, ortho_tol, "classical_gs"); \
- } \
- SECTION("modified_gs") { \
- check_qr(A, [](const Matrix& M) { return linalg::qr_modified_gs(M); }, \
- recon_tol, ortho_tol, "modified_gs"); \
- } \
- SECTION("householder") { \
- check_qr(A, [](const Matrix& M) { return linalg::qr_householder(M); }, \
- recon_tol, ortho_tol, "householder"); \
+#define FOR_ALL_METHODS(A, recon_tol, ortho_tol) \
+ SECTION("classical_gs") { \
+ check_qr(A, [](const Matrix& M) { return linalgebra::qr_classical_gs(M); }, \
+ recon_tol, ortho_tol, "classical_gs"); \
+ } \
+ SECTION("modified_gs") { \
+ check_qr(A, [](const Matrix& M) { return linalgebra::qr_modified_gs(M); }, \
+ recon_tol, ortho_tol, "modified_gs"); \
+ } \
+ SECTION("householder") { \
+ check_qr(A, [](const Matrix& M) { return linalgebra::qr_householder(M); }, \
+ recon_tol, ortho_tol, "householder"); \
}
-// ---------------------------------------------------------------------------
-// Basic correctness: square matrices
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: 3x3 known matrix", "[qr]") {
const Matrix A{
{1.0, 2.0, 3.0},
{4.0, 5.0, 6.0},
- {7.0, 8.0, 10.0} // not exactly singular
+ {7.0, 8.0, 10.0}
};
FOR_ALL_METHODS(A, 1e-12, 1e-12)
}
@@ -136,10 +114,6 @@ TEST_CASE("QR: diagonal matrix", "[qr]") {
FOR_ALL_METHODS(D, 1e-14, 1e-14)
}
-// ---------------------------------------------------------------------------
-// Rectangular (tall) matrices
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: tall 5x3 random matrix", "[qr]") {
const Matrix A = random_matrix(5, 3, 7u);
FOR_ALL_METHODS(A, 1e-12, 1e-12)
@@ -150,10 +124,6 @@ TEST_CASE("QR: tall 10x4 random matrix", "[qr]") {
FOR_ALL_METHODS(A, 1e-12, 1e-12)
}
-// ---------------------------------------------------------------------------
-// Random square matrices
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: random 6x6", "[qr]") {
const Matrix A = random_matrix(6, 6, 123u);
FOR_ALL_METHODS(A, 1e-12, 1e-12)
@@ -164,14 +134,7 @@ TEST_CASE("QR: random 12x12", "[qr]") {
FOR_ALL_METHODS(A, 1e-11, 1e-11)
}
-// ---------------------------------------------------------------------------
-// Nearly dependent columns — GS methods degrade; Householder stays clean
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: nearly dependent columns", "[qr]") {
- // Column 1 = column 0 + epsilon * e_1.
- // Classical GS will lose most of Q's orthogonality here.
- // Modified GS is better. Householder is unaffected.
constexpr double eps = 1e-7;
const Matrix A{
{1.0, 1.0 + eps, 0.0},
@@ -180,32 +143,24 @@ TEST_CASE("QR: nearly dependent columns", "[qr]") {
{0.0, 0.0, 1.0}
};
- // All three should reconstruct A accurately.
SECTION("classical_gs reconstruction") {
- const QRResult qr = linalg::qr_classical_gs(A);
+ const QRResult qr = linalgebra::qr_classical_gs(A);
CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10));
- // Orthogonality will be poor for classical GS on this input.
- // We only assert it's not catastrophically wrong (< 0.01).
CHECK(orthogonality_error(qr) < 0.01);
}
SECTION("modified_gs reconstruction") {
- const QRResult qr = linalg::qr_modified_gs(A);
+ const QRResult qr = linalgebra::qr_modified_gs(A);
CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10));
CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-8));
}
SECTION("householder reconstruction") {
- const QRResult qr = linalg::qr_householder(A);
+ const QRResult qr = linalgebra::qr_householder(A);
CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-13));
CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13));
}
}
-// ---------------------------------------------------------------------------
-// Hilbert-like ill-conditioned matrix
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: 4x4 Hilbert matrix", "[qr]") {
- // H[i][j] = 1 / (i + j + 1)
const std::size_t n = 4;
Matrix H(n, n);
for (std::size_t i = 0; i < n; ++i)
@@ -213,39 +168,34 @@ TEST_CASE("QR: 4x4 Hilbert matrix", "[qr]") {
H(i, j) = 1.0 / static_cast<double>(i + j + 1);
SECTION("classical_gs") {
- const QRResult qr = linalg::qr_classical_gs(H);
+ const QRResult qr = linalgebra::qr_classical_gs(H);
CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12));
- // Orthogonality is imperfect on Hilbert matrices with classical GS.
CHECK(orthogonality_error(qr) < 1e-8);
}
SECTION("modified_gs") {
- const QRResult qr = linalg::qr_modified_gs(H);
+ const QRResult qr = linalgebra::qr_modified_gs(H);
CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12));
CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-10));
}
SECTION("householder") {
- const QRResult qr = linalg::qr_householder(H);
+ const QRResult qr = linalgebra::qr_householder(H);
CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-13));
CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13));
}
}
-// ---------------------------------------------------------------------------
-// Structure checks
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: R is upper triangular", "[qr]") {
const Matrix A = random_matrix(5, 5, 555u);
- CHECK(r_is_upper_triangular(linalg::qr_classical_gs(A).R));
- CHECK(r_is_upper_triangular(linalg::qr_modified_gs(A).R));
- CHECK(r_is_upper_triangular(linalg::qr_householder(A).R));
+ CHECK(r_is_upper_triangular(linalgebra::qr_classical_gs(A).R));
+ CHECK(r_is_upper_triangular(linalgebra::qr_modified_gs(A).R));
+ CHECK(r_is_upper_triangular(linalgebra::qr_householder(A).R));
}
TEST_CASE("QR: Q columns are unit length", "[qr]") {
const Matrix A = random_matrix(6, 4, 321u);
- for (QRFn fn : {QRFn{[](const Matrix& M) { return linalg::qr_classical_gs(M); }},
- QRFn{[](const Matrix& M) { return linalg::qr_modified_gs(M); }},
- QRFn{[](const Matrix& M) { return linalg::qr_householder(M); }}}) {
+ for (QRFn fn : {QRFn{[](const Matrix& M) { return linalgebra::qr_classical_gs(M); }},
+ QRFn{[](const Matrix& M) { return linalgebra::qr_modified_gs(M); }},
+ QRFn{[](const Matrix& M) { return linalgebra::qr_householder(M); }}}) {
const QRResult qr = fn(A);
for (std::size_t j = 0; j < qr.Q.cols(); ++j) {
double norm2 = 0.0;
@@ -256,25 +206,20 @@ TEST_CASE("QR: Q columns are unit length", "[qr]") {
}
}
-// ---------------------------------------------------------------------------
-// Failure cases
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR: fat matrix throws DimensionMismatchError", "[qr]") {
- const Matrix A(3, 5); // rows < cols
- CHECK_THROWS_AS(linalg::qr_classical_gs(A), DimensionMismatchError);
- CHECK_THROWS_AS(linalg::qr_modified_gs(A), DimensionMismatchError);
- CHECK_THROWS_AS(linalg::qr_householder(A), DimensionMismatchError);
+ const Matrix A(3, 5);
+ CHECK_THROWS_AS(linalgebra::qr_classical_gs(A), DimensionMismatchError);
+ CHECK_THROWS_AS(linalgebra::qr_modified_gs(A), DimensionMismatchError);
+ CHECK_THROWS_AS(linalgebra::qr_householder(A), DimensionMismatchError);
}
TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") {
const Matrix A{
- {1.0, 2.0, 2.0}, // col 2 = 2 * col 0
+ {1.0, 2.0, 2.0},
{2.0, 4.0, 4.0},
{3.0, 6.0, 6.0}
};
- CHECK_THROWS_AS(linalg::qr_classical_gs(A), SingularMatrixError);
- CHECK_THROWS_AS(linalg::qr_modified_gs(A), SingularMatrixError);
- // Householder handles rank-deficient input gracefully (R gets a zero diagonal entry).
- CHECK_NOTHROW(linalg::qr_householder(A));
+ CHECK_THROWS_AS(linalgebra::qr_classical_gs(A), SingularMatrixError);
+ CHECK_THROWS_AS(linalgebra::qr_modified_gs(A), SingularMatrixError);
+ CHECK_NOTHROW(linalgebra::qr_householder(A));
}
diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp
index 7126334..008f109 100644
--- a/tests/test_qr_iteration.cpp
+++ b/tests/test_qr_iteration.cpp
@@ -1,7 +1,4 @@
-#include "linalg_error.hpp"
-#include "matrix.hpp"
-#include "qr_iteration.hpp"
-#include "vector.hpp"
+import linalgebra;
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -16,19 +13,14 @@
#include <utility>
#include <vector>
-using linalg::Matrix;
-using linalg::NonConvergenceError;
-using linalg::QRIterationOptions;
-using linalg::QRIterationResult;
-using linalg::Vector;
-
-// ---------------------------------------------------------------------------
-// Test helpers
-// ---------------------------------------------------------------------------
+using linalgebra::Matrix;
+using linalgebra::NonConvergenceError;
+using linalgebra::QRIterationOptions;
+using linalgebra::QRIterationResult;
+using linalgebra::Vector;
namespace {
-// Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag.
using EigPairs = std::vector<std::pair<double, double>>;
EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) {
@@ -66,20 +58,6 @@ bool eigs_match(const Vector& computed_real, const Vector& computed_imag,
} // namespace
-// ---------------------------------------------------------------------------
-// Test 1: 2×2 symmetric matrix with known eigenvalues
-// ---------------------------------------------------------------------------
-//
-// A = | 2 1 | is symmetric positive definite.
-// | 1 2 |
-//
-// Characteristic polynomial: (2-λ)^2 - 1 = 0 → λ = 1, 3.
-// Eigenvectors: [1,-1]/√2 (λ=1) and [1,1]/√2 (λ=3).
-//
-// The unshifted iteration converges at rate |λ_1/λ_2| = 1/3 per step,
-// so only a handful of iterations are needed.
-// Ref: T&B Theorem 28.2.
-
TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues",
"[qr_iteration][shifted]") {
const Matrix A{
@@ -87,13 +65,12 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues",
{1.0, 2.0}
};
- const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(A);
REQUIRE(res.eigenvalues_real.size() == 2);
REQUIRE(res.eigenvalues_imag.size() == 2);
REQUIRE(res.iterations > 0);
- // All eigenvalues of a symmetric matrix must be real.
CHECK(std::abs(res.eigenvalues_imag[0]) < 1e-8);
CHECK(std::abs(res.eigenvalues_imag[1]) < 1e-8);
@@ -101,23 +78,6 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues",
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
-// ---------------------------------------------------------------------------
-// Test 2: 4×4 symmetric tridiagonal matrix — reference eigenvalues
-// ---------------------------------------------------------------------------
-//
-// The n×n symmetric tridiagonal matrix with 2 on the diagonal and -1 on the
-// first super- and sub-diagonals has known eigenvalues (discrete Laplacian):
-//
-// λ_k = 2 - 2 cos(k π / (n+1)), k = 1, …, n
-//
-// Ref: Golub & Van Loan §4.4.2 (discrete sine transform).
-//
-// For n = 4:
-// λ_1 = 2 - 2 cos(π/5) ≈ 0.3820
-// λ_2 = 2 - 2 cos(2π/5) ≈ 1.3820
-// λ_3 = 2 - 2 cos(3π/5) ≈ 2.6180
-// λ_4 = 2 - 2 cos(4π/5) ≈ 3.6180
-
TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal",
"[qr_iteration][unshifted]") {
const Matrix A{
@@ -127,16 +87,14 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal",
{ 0.0, 0.0, -1.0, 2.0}
};
- const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(A);
REQUIRE(res.eigenvalues_real.size() == 4);
REQUIRE(res.eigenvalues_imag.size() == 4);
- // All eigenvalues of a symmetric matrix must be real.
for (std::size_t k = 0; k < 4; ++k)
CHECK(std::abs(res.eigenvalues_imag[k]) < 1e-8);
- // Compare against the closed-form reference.
constexpr double pi = 3.14159265358979323846;
const EigPairs expected = {
{2.0 - 2.0 * std::cos( pi / 5.0), 0.0},
@@ -147,21 +105,6 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal",
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
-// ---------------------------------------------------------------------------
-// Test 3: 5×5 symmetric tridiagonal — convergence history
-// ---------------------------------------------------------------------------
-//
-// Uses a 5×5 symmetric tridiagonal (discrete Laplacian) to guarantee all
-// real eigenvalues and predictable linear convergence. The Frobenius norm
-// of the strict lower triangle is printed at every step so the convergence
-// rate can be observed directly.
-//
-// Expected behaviour: ||lower(A_k)||_F decreases geometrically each step
-// (linear convergence), with ratio ≈ max_j |λ_{j+1}/λ_j|.
-// Ref: T&B Theorem 28.2.
-//
-// 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5.
-
TEST_CASE("QR iteration (unshifted): 5x5 convergence history",
"[qr_iteration][unshifted]") {
const Matrix A{
@@ -175,12 +118,11 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history",
QRIterationOptions opts;
opts.track_convergence = true;
- const QRIterationResult res = linalg::eigenvalues_unshifted(A, opts);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(A, opts);
REQUIRE_FALSE(res.convergence_history.empty());
REQUIRE(res.eigenvalues_real.size() == 5);
- // Print convergence history so the linear rate is visible.
std::cout << "\n=== Unshifted QR — 5x5 convergence history ===\n";
std::cout << " Converged in " << res.iterations << " iteration(s)\n";
for (std::size_t k = 0; k < res.convergence_history.size(); ++k) {
@@ -190,32 +132,26 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history",
}
std::cout << "=======================================================\n";
- // The final recorded norm must be below the default tolerance.
CHECK(res.convergence_history.back() < opts.tolerance);
}
-// ---------------------------------------------------------------------------
-// Test 4: Eigenvalue residuals below 1e-8
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR iteration (unshifted): residuals below 1e-8",
"[qr_iteration][unshifted]") {
SECTION("2x2: eigenvalues 1 and 3") {
const Matrix A{{2.0, 1.0}, {1.0, 2.0}};
- const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(A);
const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}};
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
SECTION("3x3 diagonal: eigenvalues 1, 4, 9") {
- // Diagonal matrix — already in Schur form; converges in one step.
const Matrix D{
{1.0, 0.0, 0.0},
{0.0, 4.0, 0.0},
{0.0, 0.0, 9.0}
};
- const QRIterationResult res = linalg::eigenvalues_unshifted(D);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(D);
const EigPairs expected = {{1.0, 0.0}, {4.0, 0.0}, {9.0, 0.0}};
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
@@ -234,13 +170,13 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8",
{2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0},
{2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0}
};
- const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(A);
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
SECTION("5x5 identity: all eigenvalues == 1") {
const Matrix I = Matrix::identity(5);
- const QRIterationResult res = linalg::eigenvalues_unshifted(I);
+ const QRIterationResult res = linalgebra::eigenvalues_unshifted(I);
REQUIRE(res.eigenvalues_real.size() == 5);
for (std::size_t k = 0; k < 5; ++k) {
CHECK(std::abs(res.eigenvalues_real[k] - 1.0) < 1e-8);
@@ -249,35 +185,21 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8",
}
}
-// ---------------------------------------------------------------------------
-// Failure cases
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR iteration (unshifted): non-square matrix throws",
"[qr_iteration][unshifted]") {
- const Matrix A(3, 4); // non-square
- CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A),
- linalg::DimensionMismatchError);
+ const Matrix A(3, 4);
+ CHECK_THROWS_AS(linalgebra::eigenvalues_unshifted(A),
+ linalgebra::DimensionMismatchError);
}
TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws",
"[qr_iteration][unshifted]") {
- // Cap at zero iterations — any non-trivial matrix fails immediately.
const Matrix A{{2.0, 1.0}, {1.0, 2.0}};
QRIterationOptions opts;
opts.max_iterations = 0;
- CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A, opts), NonConvergenceError);
+ CHECK_THROWS_AS(linalgebra::eigenvalues_unshifted(A, opts), NonConvergenceError);
}
-// ===========================================================================
-// Wilkinson-shifted QR iteration
-// ===========================================================================
-
-// ---------------------------------------------------------------------------
-// Helper builds a random symmetric matrix via A = M + M^T (guaranteed real
-// eigenvalues) with a fixed seed for reproducibility.
-// ---------------------------------------------------------------------------
-
namespace {
Matrix random_symmetric(std::size_t n, unsigned seed = 42) {
@@ -296,14 +218,6 @@ Matrix random_symmetric(std::size_t n, unsigned seed = 42) {
} // namespace
-// ---------------------------------------------------------------------------
-// Test S1: shifted vs unshifted iteration count on the same matrix.
-//
-// Wilkinson-shifted QR converges (typically cubically) in far fewer steps
-// than the unshifted algorithm (linear convergence).
-// The test asserts the shifted count is strictly smaller and prints both.
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR iteration (shifted): fewer iterations than unshifted",
"[qr_iteration][shifted]") {
const Matrix A{
@@ -318,8 +232,8 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted",
QRIterationOptions opts;
opts.track_convergence = true;
- const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts);
- const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts);
+ const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts);
+ const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts);
std::cout << "\n=== Shifted vs Unshifted ===\n";
std::cout << " Unshifted iterations: " << unshifted.iterations << "\n";
@@ -333,25 +247,15 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted",
1e-8));
}
-// ---------------------------------------------------------------------------
-// Test S2: matrix where unshifted takes >100 iterations, shifted takes <20.
-//
-// A nearly-equal-eigenvalue symmetric matrix maximises the linear convergence
-// slowdown. Using a scaled identity perturbation: eigenvalues cluster near 1,
-// slowing unshifted (ratio ≈ 1) while the Wilkinson shift adapts instantly.
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >100",
"[qr_iteration][shifted]") {
- // 5×5 symmetric matrix with eigenvalues 1, 1.001, 1.002, 1.003, 1.004.
- // Off-diagonal entries couple them. Unshifted stalls (|λ_{j+1}/λ_j| ≈ 1).
const Matrix A = random_symmetric(5, 17u);
QRIterationOptions opts;
- opts.max_iterations = 2000;
+ opts.max_iterations = 2000;
- const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts);
- const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts);
+ const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts);
+ const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts);
std::cout << "\n=== Hard matrix ===\n";
std::cout << " Unshifted iterations: " << unshifted.iterations << "\n";
@@ -361,10 +265,6 @@ TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >10
CHECK(shifted.iterations < 20);
}
-// ---------------------------------------------------------------------------
-// Test S3: shifted eigenvalues match known values to within 1e-8.
-// ---------------------------------------------------------------------------
-
TEST_CASE("QR iteration (shifted): residuals below 1e-8",
"[qr_iteration][shifted]") {
SECTION("4x4 tridiagonal: closed-form eigenvalues") {
@@ -381,29 +281,18 @@ TEST_CASE("QR iteration (shifted): residuals below 1e-8",
{2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0},
{2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0}
};
- const QRIterationResult res = linalg::eigenvalues_shifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_shifted(A);
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
SECTION("2x2 known eigenvalues") {
const Matrix A{{2.0, 1.0}, {1.0, 2.0}};
const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}};
- const QRIterationResult res = linalg::eigenvalues_shifted(A);
+ const QRIterationResult res = linalgebra::eigenvalues_shifted(A);
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
}
-// ===========================================================================
-// Hessenberg reduction + practical QR algorithm
-// ===========================================================================
-
-// ---------------------------------------------------------------------------
-// Test H1: hessenberg_reduction produces correct H and Q.
-//
-// Verify: (1) H is upper Hessenberg, (2) Q is orthogonal, (3) A = Q H Q^T.
-// Ref: GVL §7.4.2.
-// ---------------------------------------------------------------------------
-
namespace {
bool is_upper_hessenberg(const Matrix& H, double tol = 1e-10) {
@@ -421,7 +310,6 @@ double frobenius_norm(const Matrix& A) {
return std::sqrt(s);
}
-// ||A - B||_F
double diff_norm(const Matrix& A, const Matrix& B) {
double s = 0.0;
for (std::size_t i = 0; i < A.rows(); ++i)
@@ -432,7 +320,6 @@ double diff_norm(const Matrix& A, const Matrix& B) {
return std::sqrt(s);
}
-// ||Q^T Q - I||_F
double orthogonality_error(const Matrix& Q) {
const std::size_t n = Q.rows();
double err = 0.0;
@@ -451,29 +338,21 @@ double orthogonality_error(const Matrix& Q) {
TEST_CASE("Hessenberg reduction: structure and similarity",
"[qr_iteration][shifted]") {
const Matrix A = random_symmetric(6, 7u);
- const linalg::HessenbergResult hr = linalg::hessenberg_reduction(A);
+ const linalgebra::HessenbergResult hr = linalgebra::hessenberg_reduction(A);
- // H must be upper Hessenberg.
CHECK(is_upper_hessenberg(hr.H));
-
- // Q must be orthogonal.
CHECK(orthogonality_error(hr.Q) < 1e-10);
- // A = Q H Q^T ⟹ ||A - Q H Q^T||_F < tol.
- const Matrix QtHQ = hr.Q * hr.H * linalg::transpose(hr.Q);
+ const Matrix QtHQ = hr.Q * hr.H * linalgebra::transpose(hr.Q);
CHECK(diff_norm(A, QtHQ) < 1e-10);
}
-// ---------------------------------------------------------------------------
-// Test H2: eigenvalues_hessenberg agrees with eigenvalues_shifted to 1e-6.
-// ---------------------------------------------------------------------------
-
TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6",
"[qr_iteration][shifted]") {
const Matrix A = random_symmetric(8, 99u);
- const QRIterationResult ref = linalg::eigenvalues_shifted(A);
- const QRIterationResult hess = linalg::eigenvalues_hessenberg(A);
+ const QRIterationResult ref = linalgebra::eigenvalues_shifted(A);
+ const QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A);
REQUIRE(hess.eigenvalues_real.size() == 8);
CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag,
@@ -481,10 +360,6 @@ TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6",
1e-6));
}
-// ---------------------------------------------------------------------------
-// Test H3: known eigenvalues — 4×4 tridiagonal.
-// ---------------------------------------------------------------------------
-
TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix",
"[qr_iteration][shifted]") {
const Matrix A{
@@ -500,19 +375,10 @@ TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix",
{2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0},
{2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0}
};
- const QRIterationResult res = linalg::eigenvalues_hessenberg(A);
+ const QRIterationResult res = linalgebra::eigenvalues_hessenberg(A);
CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
}
-// ---------------------------------------------------------------------------
-// Test H4: benchmark — shifted QR vs Hessenberg pipeline for n = 50, 100, 200.
-//
-// The Hessenberg pipeline reduces each QR step from O(n³) to O(n²), so the
-// speedup should grow with n. We print the wall-clock ratio and assert that
-// the Hessenberg version is faster for n >= 50.
-// Ref: GVL §7.4.2; T&B Lecture 29.
-// ---------------------------------------------------------------------------
-
TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n",
"[qr_iteration][hessenberg]") {
using Clock = std::chrono::high_resolution_clock;
@@ -531,11 +397,11 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n",
const Matrix A = random_symmetric(n, 13u);
const auto t0s = Clock::now();
- { const auto tmp = linalg::eigenvalues_shifted(A); (void)tmp; }
+ { const auto tmp = linalgebra::eigenvalues_shifted(A); (void)tmp; }
const double t_shifted = Seconds(Clock::now() - t0s).count();
const auto t0h = Clock::now();
- const QRIterationResult hess = linalg::eigenvalues_hessenberg(A);
+ const QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A);
const double t_hess = Seconds(Clock::now() - t0h).count();
const double speedup = t_shifted / t_hess;
@@ -547,11 +413,9 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n",
<< std::setprecision(2)
<< std::setw(12) << speedup << "x\n";
- // The Hessenberg version must be faster for all tested sizes.
CHECK(t_hess < t_shifted);
- // And must give correct eigenvalues (agree with shifted to 1e-6).
- const QRIterationResult ref = linalg::eigenvalues_shifted(A);
+ const QRIterationResult ref = linalgebra::eigenvalues_shifted(A);
CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag,
to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag),
1e-6));
diff --git a/tests/test_triangular_solve.cpp b/tests/test_triangular_solve.cpp
index da7e55c..f1ce709 100644
--- a/tests/test_triangular_solve.cpp
+++ b/tests/test_triangular_solve.cpp
@@ -1,20 +1,17 @@
-#include "linalg_error.hpp"
-#include "matrix.hpp"
-#include "norms.hpp"
-#include "triangular_solve.hpp"
+import linalgebra;
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-using linalg::DimensionMismatchError;
-using linalg::Matrix;
-using linalg::SingularMatrixError;
-using linalg::Vector;
+using linalgebra::DimensionMismatchError;
+using linalgebra::Matrix;
+using linalgebra::SingularMatrixError;
+using linalgebra::Vector;
namespace {
double residual_norm(const Matrix& a, const Vector& x, const Vector& b) {
- return linalg::norm2((a * x) - b);
+ return linalgebra::norm2((a * x) - b);
}
} // namespace
@@ -28,7 +25,7 @@ TEST_CASE("Forward substitution solves lower-triangular systems", "[triangular]"
const Vector expected{1.0, 2.0, -1.0};
const Vector rhs = lower * expected;
- const Vector x = linalg::forward_substitution(lower, rhs);
+ const Vector x = linalgebra::forward_substitution(lower, rhs);
REQUIRE(x.size() == expected.size());
CHECK(x[0] == Catch::Approx(expected[0]));
CHECK(x[1] == Catch::Approx(expected[1]));
@@ -45,7 +42,7 @@ TEST_CASE("Backward substitution solves upper-triangular systems", "[triangular]
const Vector expected{2.0, -1.0, 3.0};
const Vector rhs = upper * expected;
- const Vector x = linalg::backward_substitution(upper, rhs);
+ const Vector x = linalgebra::backward_substitution(upper, rhs);
REQUIRE(x.size() == expected.size());
CHECK(x[0] == Catch::Approx(expected[0]));
CHECK(x[1] == Catch::Approx(expected[1]));
@@ -61,7 +58,7 @@ TEST_CASE("Triangular solves support unit-diagonal systems", "[triangular]") {
};
const Vector rhs{1.0, 0.0, 4.0};
- const Vector x = linalg::forward_substitution(lower, rhs, 1e-12, true);
+ const Vector x = linalgebra::forward_substitution(lower, rhs, 1e-12, true);
CHECK(x[0] == Catch::Approx(1.0));
CHECK(x[1] == Catch::Approx(2.0));
CHECK(x[2] == Catch::Approx(3.0));
@@ -70,19 +67,19 @@ TEST_CASE("Triangular solves support unit-diagonal systems", "[triangular]") {
TEST_CASE("Triangular solves reject shape and structure mismatches", "[triangular]") {
const Matrix nonsquare(2, 3);
const Vector rhs2{1.0, 2.0};
- CHECK_THROWS_AS(linalg::forward_substitution(nonsquare, rhs2), DimensionMismatchError);
+ CHECK_THROWS_AS(linalgebra::forward_substitution(nonsquare, rhs2), DimensionMismatchError);
const Matrix lower{
{1.0, 1.0},
{2.0, 3.0}
};
- CHECK_THROWS_AS(linalg::forward_substitution(lower, rhs2), std::invalid_argument);
+ CHECK_THROWS_AS(linalgebra::forward_substitution(lower, rhs2), std::invalid_argument);
const Matrix upper{
{1.0, 2.0},
{1.0, 3.0}
};
- CHECK_THROWS_AS(linalg::backward_substitution(upper, rhs2), std::invalid_argument);
+ CHECK_THROWS_AS(linalgebra::backward_substitution(upper, rhs2), std::invalid_argument);
}
TEST_CASE("Triangular solves detect negligible diagonal entries", "[triangular]") {
@@ -91,11 +88,11 @@ TEST_CASE("Triangular solves detect negligible diagonal entries", "[triangular]"
{2.0, 1.0}
};
const Vector rhs{1.0, 2.0};
- CHECK_THROWS_AS(linalg::forward_substitution(lower, rhs), SingularMatrixError);
+ CHECK_THROWS_AS(linalgebra::forward_substitution(lower, rhs), SingularMatrixError);
const Matrix upper{
{1.0, 2.0},
{0.0, 1e-14}
};
- CHECK_THROWS_AS(linalg::backward_substitution(upper, rhs), SingularMatrixError);
+ CHECK_THROWS_AS(linalgebra::backward_substitution(upper, rhs), SingularMatrixError);
}
diff --git a/tests/test_vector.cpp b/tests/test_vector.cpp
index 364ce31..a776dc8 100644
--- a/tests/test_vector.cpp
+++ b/tests/test_vector.cpp
@@ -1,6 +1,4 @@
-#include "linalg_error.hpp"
-#include "norms.hpp"
-#include "vector.hpp"
+import linalgebra;
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -8,8 +6,8 @@
#include <type_traits>
#include <utility>
-using linalg::DimensionMismatchError;
-using linalg::Vector;
+using linalgebra::DimensionMismatchError;
+using linalgebra::Vector;
TEST_CASE("Vector constructors initialize size and values", "[vector]") {
const Vector empty;
@@ -82,17 +80,17 @@ TEST_CASE("Vector arithmetic enforces shape compatibility", "[vector]") {
CHECK(scaled[1] == 1.0);
CHECK(scaled[2] == 1.5);
- CHECK(linalg::dot(a, b) == 32.0);
+ CHECK(linalgebra::dot(a, b) == 32.0);
const Vector short_vec{1.0, 2.0};
CHECK_THROWS_AS(a + short_vec, DimensionMismatchError);
- CHECK_THROWS_AS(linalg::dot(a, short_vec), DimensionMismatchError);
+ CHECK_THROWS_AS(linalgebra::dot(a, short_vec), DimensionMismatchError);
}
TEST_CASE("Vector 2-norm matches manually computed values", "[vector][norms]") {
const Vector v{3.0, 4.0};
- CHECK(linalg::norm2(v) == Catch::Approx(5.0));
+ CHECK(linalgebra::norm2(v) == Catch::Approx(5.0));
const Vector zero(5);
- CHECK(linalg::norm2(zero) == Catch::Approx(0.0));
+ CHECK(linalgebra::norm2(zero) == Catch::Approx(0.0));
}
generated by cgit v1.2.3 (git 2.51.0) at 2026-09-01 21:19:50 +0000