Commit 86b0d024 authored by ESTERIE Pierre's avatar ESTERIE Pierre

Road to p2m...

parent b2b14ded
......@@ -29,7 +29,7 @@ set(${CMAKE_PROJECT_NAME}_VERSION "${${CMAKE_PROJECT_NAME}_MAJOR_VERSION}.${${C
# Creating main lib
# -----------------
add_library(${CMAKE_PROJECT_NAME} INTERFACE)
target_compile_features(${CMAKE_PROJECT_NAME} INTERFACE cxx_std_20)
target_compile_features(${CMAKE_PROJECT_NAME} INTERFACE cxx_std_17)
# Set library includes
# --------------------
target_include_directories(${CMAKE_PROJECT_NAME} INTERFACE
......@@ -107,8 +107,8 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
#
# Boost
# -----
find_package(Boost 1.69.0 COMPONENTS serialization REQUIRED)
message(STATUS "${Boost_LIBRARIES}")
#find_package(Boost 1.69.0 COMPONENTS serialization REQUIRED)
#message(STATUS "${Boost_LIBRARIES}")
#
# OpenMP
......@@ -130,24 +130,24 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
# XTL, XSIMD & XTENSOR
# --------------------
# Find XSIMD properly
find_package(xsimd CONFIG REQUIRED)
find_package(xsimd REQUIRED)
if(xsimd_FOUND)
message(STATUS "XSIMD found")
target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xsimd)
else(xsimd_FOUND)
message(FATAL_ERROR "Can't find XSIMD !")
endif(xsimd_FOUND)
find_package(xtl CONFIG REQUIRED)
find_package(xtl REQUIRED)
if(xtl_FOUND)
message(STATUS "XTL found")
target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xtl)
else(xtl_FOUND)
message(FATAL_ERROR "Can't find XTL !")
endif(xtl_FOUND)
find_package(xtensor CONFIG REQUIRED)
find_package(xtensor REQUIRED)
if(xtensor_FOUND)
if(xsimd_FOUND)
list(APPEND SCALFMM_COMPILE_DEFINITIONS XTENSOR_USE_XSIMD)
set(SCALFMM_COMPILE_DEFINITIONS XTENSOR_USE_XSIMD)
endif(xsimd_FOUND)
message(STATUS "XTENSOR found")
target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xtensor)
......
......@@ -7,9 +7,6 @@
#include <utility>
#include <random>
#include <xtensor/xarray.hpp>
#include <xtensor/xio.hpp>
#include <scalfmm/container/particle.hpp>
#include <scalfmm/container/particle_container.hpp>
#include <scalfmm/interpolation/uniform.hpp>
......@@ -42,7 +39,6 @@ int main()
container_type particles_source{nb_particles};
container_type particles_target{nb_particles};
//std::random_device rd;
std::mt19937 gen(33);
std::uniform_real_distribution<double> dis(-1.0, 1.0);
......@@ -79,5 +75,6 @@ int main()
std::cout << colors::yellow << "p2m took... : "
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "us\n"
<< colors::reset;
return 0;
}
......@@ -89,7 +89,7 @@ namespace scalfmm::container
/// Space dimensions
constexpr static std::size_t dimension = Dim;
/// Size of #tuple_data_t tuple
constexpr static std::size_t NbAttributes = std::tuple_size<tuple_data_t>::value - dimension;
constexpr static std::size_t attributes_size = std::tuple_size<tuple_data_t>::value - dimension;
/// Floating point type
using value_type = ValueType;
......
......@@ -82,7 +82,7 @@ namespace scalfmm::container
inline friend point operator+(point lhs, point rhs)
{
point tmp;
for(int i = 0; i < dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
tmp[i] = lhs[i] + rhs[i];
}
......@@ -92,9 +92,9 @@ namespace scalfmm::container
/** Scalar assignment addition */
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;
}
......@@ -109,9 +109,9 @@ namespace scalfmm::container
/** Addition assignment operator */
inline point& operator-=(const point& other)
{
for(int i = 0; auto a: other)
for(std::size_t i = 0; i<dimension; ++i)
{
(*this)[++i] -= a;
(*this)[i] -= other[i];
}
return *this;
......@@ -127,7 +127,7 @@ namespace scalfmm::container
inline friend point operator-(point lhs, point rhs)
{
point tmp;
for(int i = 0; i < dimension; ++i)
for(std::size_t i = 0; i < dimension; ++i)
{
tmp[i] = lhs[i] - rhs[i];
}
......@@ -137,9 +137,9 @@ namespace scalfmm::container
/** Scalar assignment substraction*/
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;
}
......@@ -203,14 +203,25 @@ namespace scalfmm::container
}
/** Data to data division operator */
inline friend point operator/(point lhs, const point& rhs)
inline friend point& operator/(point& lhs, const point& rhs)
{
lhs /= rhs;
return lhs;
}
/** Data to data division operator */
inline friend point operator/(point lhs, point rhs)
{
point tmp;
for(std::size_t i = 0; i < dimension; ++i)
{
tmp[i] = lhs[i] / rhs[i];
}
return tmp;
}
/** Right scalar division assignment operator */
inline point& operator/=(const value_type val)
inline point& operator/=(value_type val)
{
for(std::size_t i = 0; i < dimension; ++i)
{
......
......@@ -6,7 +6,10 @@
#define SCALFMM_CORE_OPERATORS_HPP
#include <tuple>
#include <algorithm>
#include <fstream>
#include <utility>
#include <string>
#include <xsimd/xsimd.hpp>
#include <xtensor/xarray.hpp>
......@@ -40,6 +43,13 @@ namespace scalfmm::simd
noop_f{(xsimd::store_aligned(std::get<Is>(particle_ptrs), std::get<Is>(src)), 0)...};
}
template<typename Particle, typename TuplePtr, std::size_t... Is>
constexpr inline auto load_attributes(TuplePtr const& particle_ptrs, std::index_sequence<Is...> s)
{
return std::make_tuple(xsimd::load_aligned(std::get<Particle::dimension+Is>(particle_ptrs))...);
}
template<typename F, typename Container, typename... Types, std::size_t... Is>
constexpr inline auto apply_f(std::index_sequence<Is...> s, F&& f, Container&& input, Types&&... ts)
{
......@@ -51,7 +61,7 @@ namespace scalfmm::simd
constexpr inline auto generate_s(std::index_sequence<ISs...> s, Container const& cont,
std::array<std::size_t, sizeof...(ISs)> const& index)
{
return meta::multiply(std::get<ISs>(cont[index[ISs]])...);
return meta::multiply(std::get<ISs>(cont[index[ISs]])...);
}
} // namespace impl
......@@ -61,10 +71,10 @@ namespace scalfmm::simd
return impl::load_position<Position>(particle_ptrs, std::make_index_sequence<Position::dimension>{});
}
template<typename Particle, std::size_t Index, typename TuplePtr>
constexpr inline auto load_attribute(TuplePtr const& particle_ptrs)
template<typename Particle, typename TuplePtr>
constexpr inline auto load_attributes(TuplePtr const& particle_ptrs)
{
return xsimd::load_aligned(std::get<Particle::dimension + Index>(particle_ptrs));
return impl::load_attributes<Particle>(particle_ptrs, std::make_index_sequence<Particle::attributes_size>{});
}
template<typename Particle, typename Position, typename TuplePtr>
......@@ -102,9 +112,9 @@ namespace scalfmm::operators
using particle_type = typename leaf_type::particle_type;
using position_type = typename leaf_type::position_type;
using value_type = typename Leaf::value_type;
using array_type = xt::xarray<value_type>;
using simd_type = xsimd::simd_type<value_type>;
using simd_position_type = container::point<simd_type>;
using simd_potentials_type = container::point<simd_type>;
using polynomials_f_type = typename Interpolator::template polynomials_f_type<simd_type>;
......@@ -114,16 +124,15 @@ namespace scalfmm::operators
// Vectorizable particles
const std::size_t vec_size = leaf_size - leaf_size % inc;
typename array_type::shape_type pol_shape = {position_type::dimension, order};
array_type pol_at_local_pos(pol_shape);
const interpolation::map_glob_loc<simd_position_type> mapping_part_position(
leaf.center(), simd_position_type(simd_type(leaf.width())));
auto particle_iterator = leaf.begin();
constexpr std::size_t potential_size{particle_type::attributes_size};
std::vector<simd_position_type> poly_of_part(order);
std::vector<value_type> S(meta::pow(order, position_type::dimension) * inc);
std::vector<value_type, XSIMD_DEFAULT_ALLOCATOR(value_type)> S(meta::pow(order, position_type::dimension) * inc);
for(std::size_t part = 0; part < vec_size; part += inc)
{
......@@ -131,19 +140,26 @@ namespace scalfmm::operators
simd_position_type part_position(simd::load_position<simd_position_type>(particle_iterator));
simd_position_type local_position{};
mapping_part_position(part_position, local_position);
//if(part == 0) std::cout << local_position << '\n';
std::size_t o = 0;
for(std::size_t o = 0; auto& p: poly_of_part)
for(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, local_position, order, o);
//if(part == 0) std::cout << poly_of_part[o] << '\n';
o++;
}
// Assembling S
std::size_t idx{0};
std::size_t idx = 0;
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...});
auto construct_s = [/*&part,*/ &S, &poly_of_part, &idx, inc](auto&... current_indices) {
//meta::print(std::cout, current_indices...);
//std::cout << '\n';
auto s_simd = simd::generate_s<position_type::dimension>(poly_of_part, {{current_indices...}});
//if(part == 0) std::cout << s_simd[0] << '\n';
xsimd::store_aligned(&S[idx], s_simd);
idx += inc;
};
......@@ -153,6 +169,9 @@ namespace scalfmm::operators
meta::looper<position_type::dimension>()(construct_s, stops);
auto simd_potentials = simd::load_attributes<particle_type>(particle_iterator);
//std::cout << std::get<0>(simd_potentials) << '\n';
poly_of_part.clear();
S.clear();
particle_iterator++;
......
......@@ -16,13 +16,13 @@ namespace scalfmm::interpolation
constexpr auto id(T t)
{
return std::forward<T>(t);
};
}
template<typename Gen, std::size_t... I>
constexpr auto get_generator(Gen&& gen, std::index_sequence<I...>)
{
return std::move(xt::stack(std::make_tuple(id<I>(gen)...)));
};
}
template<std::size_t dim>
constexpr auto get(const std::size_t order)
......
......@@ -65,7 +65,7 @@ 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);
loc_pos = (inner_type(2.0) * glob_pos - m_b - m_a) / (m_b - m_a);
}
private:
......
......@@ -72,7 +72,7 @@ namespace scalfmm::interpolation
value_type scale{};
const int omn{int(order) - int(n) - 1};
if(omn % 2 == 0)
if(omn % 2 != 0)
{
scale = value_type(-1.); // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
}
......
......@@ -10,6 +10,12 @@
namespace scalfmm::meta
{
template<typename... Args>
void print(std::ostream& out, Args&&... args)
{
((out << ' ' << std::forward<Args>(args)),...);
}
template<typename... Args>
constexpr inline auto multiply(Args&... args)
{
......@@ -21,12 +27,10 @@ namespace scalfmm::meta
{
template<typename F, std::size_t Dimension, typename... Is>
constexpr inline void operator()(F&& f, std::array<std::size_t,Dimension> const& stops, Is const& ... is)
constexpr inline void operator()(F&& f, std::array<std::size_t,Dimension> const& stops, Is& ... is)
{
//std::cout << "stops " << N << " : " << std::get<N-1>(stops) << '\n';
for(std::size_t i=0; i<std::get<N-1>(stops); ++i)
{
//std::cout << N << " : " << i << '\n';
looper<N-1>()(std::forward<F>(f), stops, i, is...);
}
}
......@@ -36,13 +40,11 @@ namespace scalfmm::meta
struct looper<1>
{
template<typename F, std::size_t Dimension, typename... Is>
constexpr inline void operator()(F&& f, std::array<std::size_t,Dimension> const& stops, Is const& ... is)
constexpr inline void operator()(F&& f, std::array<std::size_t,Dimension> const& stops, Is& ... is)
{
//std::cout << "stops 1" << " : " << std::get<0>(stops) << '\n';
for(std::size_t i=0; i<std::get<0>(stops); ++i)
{
//std::cout << '1' << " : " << i << '\n';
f(is...);
f(i,is...);
}
}
};
......
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