Commit 9e69e054 authored by ESTERIE Pierre's avatar ESTERIE Pierre

Revert "Some correction for simd p2m"

This reverts commit cd4b2a0d.
parent cd4b2a0d
......@@ -24,7 +24,7 @@ int main()
constexpr size_type particle_dim{3};
//constexpr size_type potential_dim{3};
constexpr size_type order{6};
constexpr size_type order{4};
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{10000};
const size_type nb_particles{5000};
const size_type tree_height{5};
const double width{3.723};
......@@ -58,26 +58,24 @@ 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() << "us\n"
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "ms\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,13 +3,14 @@
#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>
......@@ -37,17 +38,20 @@ 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;
}
......@@ -56,7 +60,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]);
}
......@@ -65,7 +69,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];
}
......@@ -82,9 +86,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[i] = lhs[i] + rhs[i];
tmp = lhs[i] + rhs[i];
}
return tmp;
}
......@@ -127,9 +131,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[i] = lhs[i] - rhs[i];
tmp = lhs[i] - rhs[i];
}
return tmp;
}
......@@ -162,7 +166,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;
......@@ -171,7 +175,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;
}
......@@ -212,9 +216,9 @@ namespace scalfmm::container
/** Right scalar division assignment operator */
inline point& operator/=(const value_type val)
{
for(std::size_t i = 0; i < dimension; ++i)
for(auto& a: *this)
{
(*this)[i] /= val;
a /= val;
}
return *this;
}
......@@ -227,8 +231,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();
......@@ -245,8 +249,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();
......@@ -263,8 +267,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();
......@@ -281,8 +285,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,11 +132,13 @@ 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, local_position, order, ++o);
p = simd::apply_f<simd_position_type::dimension>(polynomials, part_position, order, ++o);
}
// Assembling S
......@@ -144,6 +146,8 @@ 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;
};
......@@ -152,10 +156,8 @@ namespace scalfmm::operators
stops.fill(order);
meta::looper<position_type::dimension>()(construct_s, stops);
poly_of_part.clear();
S.clear();
particle_iterator++;
}
}
......
......@@ -5,6 +5,7 @@
#ifndef SCALFMM_INTERPOLATION_MAPPING_HPP
#define SCALFMM_INTERPOLATION_MAPPING_HPP
namespace scalfmm::interpolation
{
template<typename T>
......@@ -15,8 +16,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))
{
}
......@@ -65,7 +66,8 @@ namespace scalfmm::interpolation
inline void operator()(const value_type& glob_pos, value_type& loc_pos) const override
{
loc_pos = (inner_type(2.0) * glob_pos - m_b - m_a) / (m_b - m_a);
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);
}
private:
......@@ -73,6 +75,6 @@ namespace scalfmm::interpolation
using base_type::m_b;
};
} // namespace scalfmm::interpolation
#endif // SCALFMM_INTERPOLATION_MAPPING_HPP
}
#endif // SCALFMM_INTERPOLATION_MAPPING_HPP
......@@ -7,45 +7,16 @@
#include <cmath>
#include <limits>
#include <type_traits>
#include <scalfmm/interpolation/interpolator.hpp>
#include <scalfmm/utils/math.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xtensor.hpp>
#include <xsimd/xsimd.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>
{
......@@ -60,17 +31,24 @@ namespace scalfmm::interpolation
inline value_type operator()(value_type x, const std::size_t order, const std::size_t n) override
{
//assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
//using bool_type = typename xsimd::simd_batch_traits<value_type>::batch_bool_type;
impl::compare(x);
assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
//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)
value_type scale{};
const int omn{int(order) - int(n) - 1};
if(omn % 2 == 0)
if(omn % 2)
{
scale = value_type(-1.); // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
}
......@@ -94,11 +72,12 @@ 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