Commit 145524c3 authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

Merge branch 'experimental' into 'simd-perf-fix'

# Conflicts:
#   experimental/include/scalfmm/algorithms/full_direct.hpp
#   experimental/include/scalfmm/matrix_kernels/laplace.hpp
parents 6062b14d d4654444
......@@ -5,10 +5,13 @@
#ifndef SCALFMM_ALGORITHMS_FULL_DIRECT_HPP
#define SCALFMM_ALGORITHMS_FULL_DIRECT_HPP
#include <array>
#include <iostream>
#include <iterator>
//
#include "scalfmm/meta/traits.hpp"
#include "scalfmm/meta/utils.hpp"
#include <scalfmm/container/iterator.hpp>
#include <scalfmm/container/point.hpp>
#include <scalfmm/matrix_kernels/laplace.hpp>
......@@ -30,9 +33,9 @@ namespace scalfmm::algorithms
{
// using particle_type = Particle;
using iterator_type = decltype(begin);
using value_type = typename scalfmm::container::iterator_traits<iterator_type>::inputs_type;
using range_outputs_type = typename scalfmm::container::iterator_traits<iterator_type>::range_outputs_type;
using particle_type = typename scalfmm::container::iterator_traits<iterator_type>::particle_type;
using value_type = typename particle_type::position_type::value_type;
static_assert(MatrixKernel::km == scalfmm::container::iterator_traits<iterator_type>::inputs_size,
"Different input size between Matrix kernel and container!");
static_assert(MatrixKernel::kn == scalfmm::container::iterator_traits<iterator_type>::outputs_size,
......@@ -46,7 +49,7 @@ namespace scalfmm::algorithms
// Get particle position
auto pt_x = particle_type(*it_p).position();
// val the array of outputs
std::array<double, MatrixKernel::kn> val{};
std::array<value_type, MatrixKernel::kn> val{};
auto compute = [&pt_x, &val, &matrix_kernel](ContainerIterator begin, ContainerIterator end) {
for(auto it_p2 = begin; it_p2 < end; ++it_p2)
{
......@@ -89,6 +92,70 @@ namespace scalfmm::algorithms
// ++idx;
}
}
template<typename ContainerParticles, typename MatrixKernel>
inline auto full_direct(ContainerParticles particles, const MatrixKernel& matrix_kernel)
{
// using particle_type = Particle;
using particle_type = typename ContainerParticles::value_type;
using value_type = typename particle_type::position_type::value_type;
using range_outputs_type = typename particle_type::range_outputs_type;
static_assert(MatrixKernel::km == particle_type::inputs_size,
"Different input size between Matrix kernel and container!");
static_assert(MatrixKernel::kn == particle_type::outputs_size,
"Different output size between Matrix kernel and container!");
// int idx{};
const int dimension = particle_type::position_type::dimension;
//#pragma omp parallel for shared(matrix_kernel)
for(std::size_t idx = 0; idx < particles.size(); ++idx)
{
// Get particle position
auto pt_x = particles.at(idx).position();
// val the array of outputs
std::array<value_type, MatrixKernel::kn> val{};
auto compute = [&pt_x, &val, &matrix_kernel, &particles](std::size_t start, std::size_t end) {
for(std::size_t idx_2 = start; idx_2 < end; ++idx_2)
{
auto q = particles.at(idx_2).inputs();
auto pt_y = particles.at(idx_2).position();
auto val_mat = matrix_kernel.evaluate(pt_x, pt_y);
for(int j = 0; j < MatrixKernel::kn; ++j)
{
for(int i = 0; i < MatrixKernel::km; ++i)
{
val[j] += val_mat.at(j, i) * q[i];
}
}
}
};
compute(0, idx);
compute(idx + 1, particles.size());
//
auto out = particles.at(idx).outputs();
for(int j = 0; j < MatrixKernel::km; ++j)
{
out[j] = val[j];
}
// std::cout << idx << " " << std::distance(begin, it_p) << " "
// << val[0]<<std::endl;
// auto outputs_tuple = meta::sub_tuple(particles.at(idx), range_outputs_type{});
// outputs_tuple = meta::to_tuple(val);
//
// const auto& p = particle_type(*it_p);
// std::cout << idx << " pos: " << pt_x << " i ";
// for (auto& e : p.inputs()) {
// std::cout << " o " << e;
// }
// for (auto& e : particle_type(*it_p).outputs()) {
// std::cout << " " << e;
// }
// std::cout << std::endl;
// ++idx;
}
}
} // namespace scalfmm::algorithms
......
......@@ -16,6 +16,7 @@
#include <xsimd/math/xsimd_scalar.hpp>
#include <xtensor/xmath.hpp>
#include <xtensor/xfixed.hpp>
#include <xtensor/xmath.hpp>
#include "scalfmm/meta/utils.hpp"
#include <scalfmm/container/point.hpp>
......@@ -76,7 +77,8 @@ namespace scalfmm::matrix_kernels::laplace
std::index_sequence<Is...> is) const noexcept
{
using std::sqrt;
return matrix_type<ValueType>{ValueType(1.0) / sqrt((((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...))};
return matrix_type<ValueType>{ValueType(1.0) /
sqrt((((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...))};
}
const std::size_t separation_criterion{1};
};
......@@ -288,7 +290,7 @@ namespace scalfmm::matrix_kernels::laplace
static constexpr auto homogeneity_tag{homogeneity::homogenous};
static constexpr auto symmetry_tag{symmetry::symmetric};
static constexpr std::size_t km{1};
static constexpr std::size_t kn{1+Dim};
static constexpr std::size_t kn{1 + Dim};
template<typename ValueType>
using matrix_type = std::array<ValueType, kn*km>;
template<typename ValueType>
......@@ -346,8 +348,7 @@ namespace scalfmm::matrix_kernels::laplace
///
/// The kernel \f$k(x,y): R^{1} -> R^{kn}\f$ with\f$ kn = km = 1\f$
/// \f$k(x,y) = \ln| x - y |\f$
/// The kernel is homogeneous\f$ k(ax,ay) = 1/a k(x,y)\f$
/// scale factor does not exist for non_homogenous kernel
/// The kernel is non homogeneous
///
struct ln_2d
{
......@@ -391,6 +392,57 @@ namespace scalfmm::matrix_kernels::laplace
const std::size_t separation_criterion{1};
};
///////
/// \brief The grad ln_2d struct corresponds to the \f$ ld/dr log(r) \f$ kernel
///
/// The kernel \f$k(x,y): R^{1} -> R^{kn}\f$ with\f$ kn =2 ; km = 1\f$
/// \f$k(x,y) = (x-y)/| x - y |^2\f$
/// The kernel is homogeneous
/// \f$k(ax,ay) = 1/a k(x,y)\f$
///
#ifdef TO_COMPILE
struct grad_ln_2d
{
static constexpr auto homogeneity_tag{homogeneity::homogenous};
static constexpr auto symmetry_tag{symmetry::symmetric};
static constexpr std::size_t km{1};
static constexpr std::size_t kn{2};
template<typename ValueType>
using matrix_type = xt::xtensor_fixed<ValueType, xt::xshape<kn, km>>;
template<typename ValueType>
using vector_type = xt::xtensor_fixed<ValueType, xt::xshape<kn>>;
template<typename ValueType>
[[nodiscard]] inline constexpr auto mutual_coefficient() const
{
return vector_type<ValueType>({ValueType(1.)});
}
template<typename ValueType>
[[nodiscard]] inline auto evaluate(container::point<ValueType, 2> const& x,
container::point<ValueType, 2> const& y) const noexcept
{
ValueType tmp = ValueType(1.0) / ((((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...));
return matrix_type<ValueType>({{(tmp * (ys.at(Is) - xs.at(Is)))}...}); // r2 * (xs - ys);
}
template<typename ValueType>
[[nodiscard]] inline auto scale_factor(ValueType cell_width) const noexcept
{
auto tmp{ValueType(2.) / cell_width};
vector_type<ValueType> sf;
std::fill(std::begin(sf), std::end(sf), tmp);
return sf;
}
template<typename ValueType>
[[nodiscard]] inline auto scale_factor(ValueType cell_width, std::size_t tree_level) const noexcept
{
return scale_factor(cell_width / math::pow(ValueType(2.), tree_level));
}
const std::size_t separation_criterion{1};
};
#endif
} // namespace scalfmm::matrix_kernels::laplace
#endif // SCALFMM_MATRIX_KERNELS_LAPLACE_HPP
......@@ -56,7 +56,7 @@ namespace scalfmm::component
std::swap(m_c1[i], m_c2[i]);
}
}
m_center = (m_c1 + m_c2) / 2;
m_center = (m_c1 + m_c2) * 0.5;
}
public:
......@@ -107,7 +107,7 @@ namespace scalfmm::component
d += side_length;
}
m_center = (m_c1 + m_c2) / 2;
m_center = (m_c1 + m_c2) * 0.5;
}
/** Builds a cube using its center and width
......@@ -186,7 +186,7 @@ namespace scalfmm::component
*/
[[nodiscard]] auto contains(const position_type& p) const -> bool
{
for(std::size_t i = 0; i < dimension; i++)
for(std::size_t i = 0; i < dimension; ++i)
{
if(p[i] < m_c1[i] || p[i] > m_c2[i])
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment