1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
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);
}
|