From 907604afdbd7405327b8ef4184bf05bcad3f0eec Mon Sep 17 00:00:00 2001
From: Olivier Coulaud <olivier.coulaud@inria.fr>
Date: Fri, 21 Mar 2025 14:59:56 +0100
Subject: [PATCH] define LOW_RANK_NO_WEIGHTS added
---
include/scalfmm/utils/low_rank.hpp | 26 ++++++++++++++++----------
1 file changed, 16 insertions(+), 10 deletions(-)
diff --git a/include/scalfmm/utils/low_rank.hpp b/include/scalfmm/utils/low_rank.hpp
index a54fd737..4f176f1f 100644
--- a/include/scalfmm/utils/low_rank.hpp
+++ b/include/scalfmm/utils/low_rank.hpp
@@ -64,7 +64,12 @@ namespace scalfmm::low_rank
{
for(std::size_t i = Ib; i < Ie; ++i)
{
+#ifdef LOW_RANK_NO_WEIGHTS
+ view.at(idx) = mk.evaluate(X[i], Y[j]).at(kn * MatrixKernelType::km + km);
+#else
view.at(idx) = weights[i] * weights[j] * mk.evaluate(X[i], Y[j]).at(kn * MatrixKernelType::km + km);
+
+#endif
++idx;
}
}
@@ -164,7 +169,7 @@ namespace scalfmm::low_rank
template<typename MatrixType, typename ValueType>
void check_reconstruction(MatrixType const& K, const MatrixType& U, const MatrixType& V, const ValueType epsilon,
- bool stop = false)
+ std::string help, bool stop = false)
{
auto nnodes = K.shape()[0];
auto normK{xt::linalg::norm(K)};
@@ -177,8 +182,8 @@ namespace scalfmm::low_rank
if(error_recom > epsilon)
{
- std::cout << cpp_tools::colors::red << "error rank " << U.shape()[1] << " rel error " << error_recom
- << " epsilon " << epsilon << " ratio " << error_recom / epsilon << std::endl
+ std::cout << cpp_tools::colors::red << help << " error rank " << U.shape()[1] << " rel error "
+ << error_recom << " epsilon " << epsilon << " ratio " << error_recom / epsilon << std::endl
<< cpp_tools::colors::reset;
if(stop)
{
@@ -231,8 +236,13 @@ namespace scalfmm::low_rank
{
for(std::size_t j{0}; j < nnodes; ++j)
{
- // normW += weights[i] * weights[i];
+// normW += weights[i] * weights[i];
+#ifdef LOW_RANK_NO_WEIGHTS
+ K.at(i, j) = mk.evaluate(X[i], Y[j]).at(kn * MatrixKernelType::km + km);
+#else
K.at(i, j) = weights[i] * weights[j] * mk.evaluate(X[i], Y[j]).at(kn * MatrixKernelType::km + km);
+
+#endif
}
}
// std::cout << " normW " << normW << " epsilon " << epsilon << " new eps/normW2 " << epsilon / normW << std::endl;
@@ -257,10 +267,6 @@ namespace scalfmm::low_rank
U.at(i, j) *= S.at(j);
}
}
-
- // std::cout << " K with weights \n";
- // check_reconstruction(K, U, xt::eval(xt::transpose(V)), epsilon, false);
- //
return std::make_tuple(U, xt::eval(xt::transpose(V)));
}
@@ -353,7 +359,7 @@ namespace scalfmm::low_rank
xt::xtensor<ValueType, 2> U_row(U);
xt::xtensor<ValueType, 2> V_row(V);
-#ifndef LOW_RANK_CHECK
+#ifdef LOW_RANK_CHECK
xt::xtensor<ValueType, 2, xt::layout_type::column_major> K;
K.resize({nnodes, nnodes});
for(std::size_t i{0}; i < nnodes; ++i)
@@ -363,7 +369,7 @@ namespace scalfmm::low_rank
K.at(i, j) = mk.evaluate(X[i], Y[j]).at(kn * MatrixKernelType::km + km);
}
}
- check_reconstruction(K, U, V, epsilon);
+ check_reconstruction(K, U, V, epsilon, "kn=" + std::to_string(kn) + "km=" + std::to_string(km));
#endif
return std::make_tuple(U_row, V_row);
}
--
GitLab