Commit 38a32ca9 authored by ESTERIE Pierre's avatar ESTERIE Pierre

Some correction for simd p2m

parent ead11a88
......@@ -24,7 +24,7 @@ int main()
constexpr size_type particle_dim{3};
//constexpr size_type potential_dim{3};
constexpr size_type order{4};
constexpr size_type order{6};
using particle_type = scalfmm::container::particle<double, particle_dim, double, double, double>;
//using tuple_type = std::tuple<double,double,double>;
......@@ -35,7 +35,7 @@ int main()
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, particle_dim>;
const size_type nb_particles{5000};
const size_type nb_particles{10000};
const size_type tree_height{5};
const double width{3.723};
......@@ -58,24 +58,26 @@ int main()
};
std::generate(particles_source.begin(), particles_source.end(), make_particle);
center = {0, 0., 0.};
std::generate(particles_target.begin(), particles_target.end(), make_particle);
auto start = std::chrono::high_resolution_clock::now();
interpolator_type s(order, tree_height);
auto end = std::chrono::high_resolution_clock::now();
std::cout << colors::green << "Interpolator construction time : "
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "ms\n"
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "us\n"
<< colors::reset;
// Construct leaf
leaf_type l(std::move(particles_source), center, width);
multipole_type multipoles{scalfmm::meta::pow(order,particle_dim)};
start = std::chrono::high_resolution_clock::now();
scalfmm::operators::apply_p2m(s, order, l, multipoles);
end = std::chrono::high_resolution_clock::now();
std::cout << colors::yellow << "p2m took... : "
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "us\n"
<< colors::reset;
return 0;
}
......@@ -3,14 +3,13 @@
#ifndef SCALFMM_CONTAINER_POINT_HPP
#define SCALFMM_CONTAINER_POINT_HPP
#include <algorithm>
#include <array>
#include <cmath>
#include <initializer_list>
#include <iterator>
#include <ostream>
#include <type_traits>
#include <initializer_list>
#include <algorithm>
#include <xsimd/xsimd.hpp>
......@@ -38,20 +37,17 @@ namespace scalfmm::container
/// Space dimension count
constexpr static const std::size_t dimension = Dim;
point(std::initializer_list<value_type> l)
{
std::copy(l.begin(), l.end(), this->begin());
}
point(std::initializer_list<value_type> l) { std::copy(l.begin(), l.end(), this->begin()); }
point() = default;
point(point const&) = default;
point(point &&) noexcept = default;
point(point&&) noexcept = default;
point& operator=(point const&) = default;
point& operator=(point &&) = default;
point& operator=(point&&) = default;
point(value_type to_splat)
{
for(std::size_t i=0; i<dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] = to_splat;
}
......@@ -60,7 +56,7 @@ namespace scalfmm::container
template<typename U>
point(point<U> const& other)
{
for(std::size_t i=0; i < dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] = value_type(other[i]);
}
......@@ -69,7 +65,7 @@ namespace scalfmm::container
/** Addition assignment operator */
inline point& operator+=(const point& other)
{
for(std::size_t i=0; i<dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] += other[i];
}
......@@ -86,9 +82,9 @@ namespace scalfmm::container
inline friend point operator+(point lhs, point rhs)
{
point tmp;
for(int i = 0; i<dimension; ++i)
for(int i = 0; i < dimension; ++i)
{
tmp = lhs[i] + rhs[i];
tmp[i] = lhs[i] + rhs[i];
}
return tmp;
}
......@@ -131,9 +127,9 @@ namespace scalfmm::container
inline friend point operator-(point lhs, point rhs)
{
point tmp;
for(int i = 0; i<dimension; ++i)
for(int i = 0; i < dimension; ++i)
{
tmp = lhs[i] - rhs[i];
tmp[i] = lhs[i] - rhs[i];
}
return tmp;
}
......@@ -166,7 +162,7 @@ namespace scalfmm::container
}
/** Data to data multiplication operator */
inline friend point operator*(point & lhs, const point& rhs)
inline friend point operator*(point& lhs, const point& rhs)
{
lhs *= rhs;
return lhs;
......@@ -175,7 +171,7 @@ namespace scalfmm::container
/** Right scalar multiplication assignment operator */
inline point& operator*=(const value_type& val)
{
for(std::size_t i=0; i<dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] *= val;
}
......@@ -216,9 +212,9 @@ namespace scalfmm::container
/** Right scalar division assignment operator */
inline point& operator/=(const value_type val)
{
for(auto& a: *this)
for(std::size_t i = 0; i < dimension; ++i)
{
a /= val;
(*this)[i] /= val;
}
return *this;
}
......@@ -231,8 +227,8 @@ namespace scalfmm::container
}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// template<class T>
// inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_float<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
......@@ -249,8 +245,8 @@ namespace scalfmm::container
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// template<class T>
// inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_double<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
......@@ -267,8 +263,8 @@ namespace scalfmm::container
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// template<class T>
// inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<std::is_integral<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
......@@ -285,8 +281,8 @@ namespace scalfmm::container
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// template<class T>
// inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_simd<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
......
......@@ -132,13 +132,11 @@ namespace scalfmm::operators
simd_position_type local_position{};
mapping_part_position(part_position, local_position);
if(part == 0) std::cout << "part position = " << part_position << '\n';
if(part == 0) std::cout << "local position = " << local_position << '\n';
for(std::size_t o = 0; auto p: poly_of_part)
for(std::size_t o = 0; auto& p: poly_of_part)
{
polynomials_f_type polynomials{};
p = simd::apply_f<simd_position_type::dimension>(polynomials, part_position, order, ++o);
p = simd::apply_f<simd_position_type::dimension>(polynomials, local_position, order, ++o);
}
// Assembling S
......@@ -146,8 +144,6 @@ namespace scalfmm::operators
auto construct_s = [&S, &poly_of_part, &idx, inc](auto&... current_indices) {
auto s_simd = simd::generate_s<position_type::dimension>(poly_of_part, {current_indices...});
//std::cout << s_simd << '\n';
//std::cout << "idx = " << idx << '\n';
xsimd::store_aligned(&S[idx], s_simd);
idx += inc;
};
......@@ -156,8 +152,10 @@ namespace scalfmm::operators
stops.fill(order);
meta::looper<position_type::dimension>()(construct_s, stops);
poly_of_part.clear();
S.clear();
particle_iterator++;
}
}
......
......@@ -5,7 +5,6 @@
#ifndef SCALFMM_INTERPOLATION_MAPPING_HPP
#define SCALFMM_INTERPOLATION_MAPPING_HPP
namespace scalfmm::interpolation
{
template<typename T>
......@@ -16,8 +15,8 @@ namespace scalfmm::interpolation
using inner_type = typename T::value_type;
explicit mapper(const value_type& center, const value_type& width)
: m_a(center - width / inner_type(2.0))
, m_b(center + width / inner_type(2.0))
: m_a(center - (width / inner_type(2.0)))
, m_b(center + (width / inner_type(2.0)))
{
}
......@@ -66,8 +65,7 @@ namespace scalfmm::interpolation
inline void operator()(const value_type& glob_pos, value_type& loc_pos) const override
{
std::cout << inner_type(2.0)*glob_pos << '\n';
loc_pos = (inner_type(2.0)*glob_pos-m_b-m_a) / (m_b-m_a);
loc_pos = (inner_type(2.0) * glob_pos - m_b - m_a) / (m_b - m_a);
}
private:
......@@ -75,6 +73,6 @@ namespace scalfmm::interpolation
using base_type::m_b;
};
}
#endif // SCALFMM_INTERPOLATION_MAPPING_HPP
} // namespace scalfmm::interpolation
#endif // SCALFMM_INTERPOLATION_MAPPING_HPP
......@@ -7,15 +7,47 @@
#include <array>
#include <cassert>
#include <cmath>
#include <limits>
#include <type_traits>
#include <scalfmm/interpolation/interpolator.hpp>
#include <scalfmm/utils/math.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xtensor.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.hpp>
namespace scalfmm::interpolation
{
namespace impl
{
template<typename T, std::size_t N>
inline void compare(xsimd::batch<T,N>& x)
{
using batch_type = xsimd::batch<T,N>;
auto booleans = xsimd::abs(x) > batch_type(1.);
if(xsimd::any(booleans))
{
x = xsimd::select(x>batch_type(1.0), batch_type(1.), x);
x = xsimd::select(x<batch_type(-1.0), batch_type(-1.), x);
}
}
template<typename T>
inline auto compare(T& x)
-> std::enable_if_t<std::is_scalar<T>::value>
{
if (std::abs(x)>1.)
{
x = (x > T( 1.) ? T( 1.) : x);
x = (x < T(-1.) ? T(-1.) : x);
}
}
}
template<typename T>
struct uniform_roots : roots_callable<T>
{
......@@ -30,24 +62,17 @@ namespace scalfmm::interpolation
inline value_type operator()(value_type x, const std::size_t order, const std::size_t n) override
{
//using bool_type = typename xsimd::simd_batch_traits<value_type>::batch_bool_type;
//assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
impl::compare(x);
//bool_type booleans = xsimd::abs(x) > value_type(1.);
//if(xsimd::any(booleans))
//{
// xsimd::select(x>value_type(1.0), value_type(1.),x);
// xsimd::select(x<value_type(-1.0), value_type(-1.),x);
//}
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
value_type scale{};
const int omn{int(order) - int(n) - 1};
if(omn % 2)
if(omn % 2 == 0)
{
scale = value_type(-1.); // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
}
......@@ -71,12 +96,11 @@ namespace scalfmm::interpolation
// new version (reducing round-off)
// regular grid on [-1,1] (h simplifies, only the size of the domain and a remains i.e. 2. and -1.)
L *= (value_type(int(order) - 1) * (x + value_type(1.)) - value_type(2.) * value_type(m));
L *= (value_type(int(order) - 1) * (x + value_type(1.)) - (value_type(2.) * value_type(m)));
}
}
L *= scale;
return L;
}
};
......
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