lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add solving/inverting lower & upper triangular systems/matrices and tests

Browse source
This commit is contained in:
Nikita Lisitsa 2023-12-28 19:01:14 +03:00
parent 05b7f2d560
commit 7611d375dc
2 changed files with 168 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

78
libs/geom/include/psemek/geom/gauss.hpp
View file

@ -243,4 +243,82 @@ namespace psemek::geom
}
}
template <typename T, std::size_t N, typename ... RHS>
bool solve_lower_triangular(matrix<T, N, N> const & m, RHS & ... rhs)
{
for (std::size_t i = 0; i < N; ++i)
{
if (m[i][i] == T{}) return false;
for (std::size_t j = 0; j < i; ++j)
{
detail::for_each([&](auto & rhs){
detail::gauss_helper h{rhs};
for (std::size_t c = 0; c < h.columns(); ++c)
{
h.at(i, c) -= h.at(j, c) * m[i][j];
}
}, rhs...);
}
detail::for_each([&](auto & rhs){
detail::gauss_helper h{rhs};
for (std::size_t c = 0; c < h.columns(); ++c)
{
h.at(i, c) /= m[i][i];
}
}, rhs...);
}
return true;
}
template <typename T, std::size_t N, typename ... RHS>
bool solve_upper_triangular(matrix<T, N, N> const & m, RHS & ... rhs)
{
for (std::size_t i = N; i --> 0;)
{
if (m[i][i] == T{}) return false;
for (std::size_t j = i + 1; j < N; ++j)
{
detail::for_each([&](auto & rhs){
detail::gauss_helper h{rhs};
for (std::size_t c = 0; c < h.columns(); ++c)
{
h.at(i, c) -= h.at(j, c) * m[i][j];
}
}, rhs...);
}
detail::for_each([&](auto & rhs){
detail::gauss_helper h{rhs};
for (std::size_t c = 0; c < h.columns(); ++c)
{
h.at(i, c) /= m[i][i];
}
}, rhs...);
}
return true;
}
template <typename T, std::size_t N>
std::optional<matrix<T, N, N>> inverse_lower_triangular(matrix<T, N, N> const & m)
{
auto result = matrix<T, N, N>::identity();
if (!solve_lower_triangular(m, result))
return std::nullopt;
return result;
}
template <typename T, std::size_t N>
std::optional<matrix<T, N, N>> inverse_upper_triangular(matrix<T, N, N> const & m)
{
auto result = matrix<T, N, N>::identity();
if (!solve_upper_triangular(m, result))
return std::nullopt;
return result;
}
}

90
libs/geom/tests/matrix.cpp
View file

@ -196,6 +196,96 @@ test_case(geom_matrix_inverse)
expect_small(frobenius_norm(*i - by_rows(vector{-2.f, 1.f}, vector{1.5f, -0.5f})), 1e-6);
}
test_case(geom_matrix_solve__lower__triangular)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 1024; ++iteration)
{
matrix<float, 8, 8> m = matrix<float, 8, 8>::zero();
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j <= i; ++j)
m[i][j] = (i == j) ? (1.f + std::abs(normal(rng))) : normal(rng);
vector<float, 8> v;
for (std::size_t i = 0; i < v.dimension(); ++i)
v[i] = normal(rng);
auto u = v;
solve_lower_triangular(m, u);
expect_small(length(m * u - v), 1e-5);
}
}
test_case(geom_matrix_solve__upper__triangular)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 1024; ++iteration)
{
matrix<float, 8, 8> m = matrix<float, 8, 8>::zero();
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = i; j < m.columns(); ++j)
m[i][j] = (i == j) ? (1.f + std::abs(normal(rng))) : normal(rng);
vector<float, 8> v;
for (std::size_t i = 0; i < v.dimension(); ++i)
v[i] = normal(rng);
auto u = v;
solve_upper_triangular(m, u);
expect_small(length(m * u - v), 1e-5);
}
}
test_case(geom_matrix_inverse__lower__triangular)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 1024; ++iteration)
{
matrix<float, 8, 8> m = matrix<float, 8, 8>::zero();
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j <= i; ++j)
m[i][j] = (i == j) ? (1.f + std::abs(normal(rng))) : normal(rng);
auto inv = inverse_lower_triangular(m);
expect(inv);
expect_small(frobenius_norm(m.identity() - m * *inv), 1e-4);
}
}
test_case(geom_matrix_inverse__upper__triangular)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 1024; ++iteration)
{
matrix<float, 8, 8> m = matrix<float, 8, 8>::zero();
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = i; j < m.columns(); ++j)
m[i][j] = (i == j) ? (1.f + std::abs(normal(rng))) : normal(rng);
auto inv = inverse_upper_triangular(m);
expect(inv);
expect_small(frobenius_norm(m.identity() - m * *inv), 1e-4);
}
}
test_case(geom_matrix_qr)
{
random::generator rng;