...
 
Commits (7)
......@@ -4,3 +4,9 @@
[submodule "inastemp"]
path = inastemp
url = https://gitlab.inria.fr/coulaud/inastemp.git
[submodule "experimental/modules/morse_cmake"]
path = experimental/modules/morse_cmake
url = https://gitlab.inria.fr/solverstack/morse_cmake.git
[submodule "experimental/modules/xsimd"]
path = experimental/modules/xsimd
url = https://github.com/QuantStack/xsimd.git
This diff is collapsed.
set(CPACK_PACKAGE_VERSION_PATCH 0)
execute_process(
COMMAND git rev-list HEAD --count
OUTPUT_VARIABLE CPACK_PACKAGE_VERSION_PATCH
RESULT_VARIABLE RET
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE
)
#string_TRIM(PATCH1 PATCH)
set(CPACK_INCLUDE_TOPLEVEL_DIRECTORY "ON")
#
set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE")
set(CPACK_PACKAGE_VERSION_MAJOR "${${CMAKE_PROJECT_NAME}_MAJOR_VERSION}")
set(CPACK_PACKAGE_VERSION_MINOR "${${CMAKE_PROJECT_NAME}_MINOR_VERSION}")
#
set(PACK_PACKAGE_VERSION "${${CMAKE_PROJECT_NAME}_VERSION}.${${CMAKE_PROJECT_NAME}_MINOR_VERSION}-${CPACK_PACKAGE_VERSION_PATCH}")
set(CPACK_SOURCE_GENERATOR "TGZ")
set(CPACK_SOURCE_PACKAGE_FILE_NAME "SCALFMM-${${CMAKE_PROJECT_NAME}_MAJOR_VERSION}.${${CMAKE_PROJECT_NAME}_MINOR_VERSION}-${CPACK_PACKAGE_VERSION_PATCH}")
set(CPACK_SOURCE_IGNORE_FILES "\\\\.git;.DS_Store;.*~;/*.aux;/*.idx;/*.log;/*.out;/*.toc;/*.ilg;CMakeLists.txt.user;/.settings/")
list(APPEND CPACK_SOURCE_IGNORE_FILES "${CMAKE_CURRENT_BINARY_DIR};${CMAKE_CURRENT_SOURCE_DIR}/Utils/;/Build*;/Obsolete/")
#
## This file should be placed in the root directory of your project.
## Then modify the CMakeLists.txt file in the root directory of your
## project to incorporate the testing dashboard.
## # The following are required to uses Dart and the Cdash dashboard
## ENABLE_TESTING()
## INCLUDE(CTest)
set(CTEST_PROJECT_NAME "ScalFMM")
set(CTEST_NIGHTLY_START_TIME "00:00:00 GMT")
set(CTEST_SUBMIT_URL "http://cdash.inria.fr/CDash/submit.php?project=scalfmm")
#--------------------------------------------------------------------
# BUILDNAME variable construction
# This variable will be used to set the build name which will appear
# on the Chameleon dashboard http://cdash.inria.fr/CDash/
#--------------------------------------------------------------------
# Start with the short system name, e.g. "Linux", "FreeBSD" or "Windows"
if(NOT BUILDNAME)
set(BUILDNAME "${CMAKE_SYSTEM_NAME}")
# Add i386 or amd64
if(CMAKE_SIZEOF_VOID_P EQUAL 8)
set(BUILDNAME "${BUILDNAME}-amd64")
else()
set(BUILDNAME "${BUILDNAME}-i386")
endif()
# Add compiler name
get_filename_component(CMAKE_CXX_COMPILER_NAME ${CMAKE_CXX_COMPILER} NAME)
set(BUILDNAME "${BUILDNAME}-${CMAKE_CXX_COMPILER_NAME}")
# Add the build type, e.g. "Debug, Release..."
if(CMAKE_BUILD_TYPE)
set(BUILDNAME "${BUILDNAME}-${CMAKE_BUILD_TYPE}")
endif(CMAKE_BUILD_TYPE)
# Specific options of Scalfmm
if(SCALFMM_USE_SSE)
set(BUILDNAME "${BUILDNAME}-sse")
endif()
if(SCALFMM_USE_BLAS)
set(BUILDNAME "${BUILDNAME}-blas")
endif()
if(SCALFMM_USE_MPI)
set(BUILDNAME "${BUILDNAME}-mpi")
endif()
endif()
# Examples
#---------
# List of source files
set(source_tests_files
sandbox.cpp
test-blas.cpp
test-xtensor.cpp
)
# Add execs - 1 cpp = 1 exec
foreach(exec ${source_tests_files})
set(compile_exec TRUE)
get_filename_component( execname ${exec} NAME_WE )
foreach(fuse_key ${FUSE_DEP_AVAILABLE})
file(STRINGS "${exec}" lines_fuse REGEX "@FUSE_${fuse_key}")
if(lines_fuse AND NOT ${fuse_key} IN_LIST FUSE_LIST)
message( STATUS "This needs ${fuse_key} = ${exec}" )
set(compile_exec FALSE)
endif()
endforeach()
# Dependency are OK
if( compile_exec )
add_executable( ${execname} ${exec})
list(APPEND SCALFMM_EXAMPLES_TARGETS ${execname})
list(APPEND SCALFMM_TESTS_TARGETS ${execname})
set_target_properties(${execname} PROPERTIES ENABLE_EXPORTS TRUE
RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BUILD_TYPE}
)
target_compile_definitions(${execname} PUBLIC -DXTENSOR_USE_XSIMD)
target_link_libraries( ${execname} ${CMAKE_PROJECT_NAME})
# TODO: Remove these include derectories, ust temporary to reuse old components
target_include_directories( ${execname} PRIVATE ${CMAKE_SOURCE_DIR}/../include/)
target_include_directories( ${execname} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
target_compile_definitions( ${execname} PRIVATE BYPASS_CONFIG)
install(TARGETS ${execname} RUNTIME DESTINATION bin)
endif()
endforeach(exec)
add_custom_target(scalfmm_examples ALL DEPENDS ${SCALFMM_EXAMPLES_TARGETS})
#ifndef TRAP_SCALFMM_CONFIG_HPP
#define TRAP_SCALFMM_CONFIG_HPP
#include <scalfmm/config/scalfmm-config.hpp>
#endif
//---------------------
// Experimental example
//---------------------
//#include <iosfwd>
//#include <string>
#include <vector>
#include <iostream>
#include <tuple>
#include <random>
#include <algorithm>
#include <inria/tcli/help_descriptor.hpp>
#include <inria/tcli/tcli.hpp>
#include <range/v3/core.hpp>
#include <range/v3/algorithm/generate.hpp>
#include <range/v3/algorithm/for_each.hpp>
#include <range/v3/algorithm/sort.hpp>
#include <range/v3/view/generate.hpp>
#include <range/v3/range/concepts.hpp>
#include <range/v3/iterator/concepts.hpp>
#include <range/v3/algorithm/result_types.hpp>
#include <range/v3/functional/concepts.hpp>
#include <range/v3/range_fwd.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor/xtensor.hpp>
#include <scalfmm/container/particle_container.hpp>
#include <scalfmm/container/point.hpp>
#include <scalfmm/container/particle.hpp>
#include <scalfmm/interpolation/uniform.hpp>
#include <scalfmm/tree/tree.hpp>
#include <scalfmm/tree/box.hpp>
#include <scalfmm/functional/utils.hpp>
#include <scalfmm/tools/colorized.hpp>
namespace args
{
struct nb_particle : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--n-particles", "-nbp"};
const char* description = "Number of particles to generate.";
using type = std::size_t;
};
struct tree_height : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--tree-height", "-th"};
const char* description = "Tree height (or initial height in case of an adaptive tree).";
using type = std::size_t;
};
struct thread_count : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--threads", "-t"};
const char* description = "Maximum thread count to be used.";
using type = std::size_t;
type def{1};
};
struct order : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--order", "-o"};
const char* description = "Precision order.";
using type = std::size_t;
};
struct file : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--input-file", "-fin"};
const char* description = "Input filename (.fma or .bfma).";
using type = std::string;
};
struct chunk_size
{
inria::tcli::str_vec flags = {"--group-size", "-gs"};
const char* description = "Group tree chunk size.";
using type = std::size_t;
type def{250};
};
auto cli = inria::tcli::make_parser(
nb_particle{},
tree_height{}, thread_count{},
file{}, inria::tcli::help{},
chunk_size{},
order{}
);
} // namespace args
/**
* \brief Store the PerfTest program parameters.
*/
struct params {
std::size_t m_nb_particle{}; ///< Tree height.
std::size_t m_tree_height{}; ///< Tree height.
std::size_t m_nb_threads{}; ///< Maximum thread count (when used).
std::string m_filename{}; ///< Particles file.
std::size_t m_group_size{}; ///< Group tree group size
std::size_t m_order{}; ///< Group tree group size
params() = default;
params(const decltype(args::cli)& cli)
{
m_order = cli.get<args::order>();
m_nb_particle = cli.get<args::nb_particle>();
m_tree_height = cli.get<args::tree_height>();
m_nb_threads = cli.get<args::thread_count>();
m_filename = cli.get<args::file>();
m_group_size = cli.get<args::chunk_size>();
}
};
using namespace scalfmm;
CPP_template(class Iter, class Sent, class Fun)
(requires //ranges::Invocable<Fun&> && ranges::OutputIterator<Iter,ranges::invoke_result_t<Fun&>> && ranges::Sentinel<Sent, Iter>
ranges::Writable<Iter,ranges::invoke_result_t<Fun&>> )
void my_algorithm(Iter first, Sent last, Fun f)
{
// ...
}
// Quick type displayer
template<typename T>
class TD;
int main(int argc, char** argv)
{
args::cli.parse(argc,argv);
params parameters(args::cli);
interpolation::uniform_interpolator<double,4> cheb_interpolator(parameters.m_order ,parameters.m_tree_height);
//variad<int,float> ot{{8,9},{1.3,2.6}};
//variad<int, float> copy{ot};
//std::cout << colors::blue << std::get<0>(copy)[0] << " " << colors::green << std::get<1>(copy)[0] << colors::reset << '\n';
//variad<int, float, double> v_{{1,2}, {3.3,4.4}, {5.5,6.6}};
//variad<int, float, double> v_(10);
//v_.push_back(std::make_tuple(2,1.1,2.2));
//std::cout << std::get<0>(v_).size() << '\n';
//std::cout << colors::blue << std::get<0>(v_)[0] << " " << colors::green << std::get<1>(v_)[0] << " " << colors::yellow << std::get<2>(v_)[0] << colors::reset << '\n';
//const std::size_t tree_height{parameters.m_tree_height};
//using source_t = container::particle<double,3,double,std::size_t>;
//using position_t = container::point<double,3>;
//using container_source_t = container::particle_container<source_t>;
//const std::size_t nb_of_part{parameters.m_nb_particle};
//container_source_t cs{nb_of_part};
//std::random_device rd;
//std::mt19937 gen(rd());
//std::uniform_real_distribution<double> dis(-1.0, 1.0);
//auto random_r = [&dis, &gen](){ return dis(gen); };
//auto make_particle = [&tree_height, &random_r]()
//{
// position_t pos = {random_r(), random_r(), random_r()};
// source_t part( pos
// , random_r()
// , index::get_morton_index(pos, component::box<position_t>{2.0, {0.0,0.0,0.0}}, tree_height)
// );
// return part.as_tuple();
//};
//auto print_particle = [](source_t p)
//{
// std::cout << colors::yellow << "=================" << colors::reset << '\n';
// std::cout << "x = " << std::get<0>(p) << '\n';
// std::cout << "y = " << std::get<1>(p) << '\n';
// std::cout << "z = " << std::get<2>(p) << '\n';
// std::cout << "p = " << std::get<3>(p) << '\n';
// std::cout << "x.p = " << p.position()[0] << '\n';
// std::cout << "norm2 = " << utils::norm2(p.position()) << '\n';
// std::cout << "morton = " << std::get<4>(p) << '\n';
//};
//auto pred = [](auto a, auto b)
//{
// return std::get<4>(a) < std::get<4>(b);
//};
//std::vector<std::tuple<int,float>> v = {{8,3.0},{9,2.1},{10,9.6}};
//my_algorithm(cs.begin(), cs.end(), make_particle);
//std::generate(cs.begin(), cs.end(), make_particle);
//std::sort(cs.begin(), cs.end(), pred);
//std::for_each(cs.begin(), cs.end(), print_particle);
//component::tree t{};
return 0;
}
// --------------------------------
// See LICENCE file at project root
// File : particle_container.cpp
// --------------------------------
#ifndef SCALFMM_CONCEPTS_LIBRARY_HPP
#define SCALFMM_CONCEPTS_LIBRARY_HPP
namespace scalfmm::concepts
{
template<typename R, typename Enabler> struct require_impl;
template<typename R> struct require_impl<R, void> { using type = R; };
template<typename Return, typename... Ts>
struct require_check : require_impl<Return,std::void_t<Ts...>>
{};
template<typename From, typename To>
using Convertible = std::enable_if_t<std::is_convertible_v<From,To>>;
template<typename T>
using Arithmetic = std::enable_if_t<std::is_arithmetic_v<T>>;
template<typename T>
using Integral = std::enable_if_t<std::is_integral_v<T>>;
template<bool Condition>
using If = std::enable_if_t<Condition>;
} // end namespace scalfmm::concepts
// Pseudo require macro
#define requires(...) -> typename ::concept::require_check<__VA_ARGS__>::type
#endif // SCALFMM_CONCEPTS_LIBRARY_HPP
// See LICENCE file at project root
///
/// File generated by cmake
// Do not remove any line
///////////////////////////////////////////////////////
#ifndef SCALFMM_CONFIG_SCALFMM_CONFIG_HPP
#define SCALFMM_CONFIG_SCALFMM_CONFIG_HPP
// Uncomment the next line to use debug mode
#cmakedefine SCALFMM_USE_LOG
///////////////////////////////////////////////////////
// Blas
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_BLAS
#cmakedefine SCALFMM_USE_MKL_AS_BLAS
// Fortran Mangling
#cmakedefine SCALFMM_BLAS_ADD_
#cmakedefine SCALFMM_BLAS_UPCASE
#cmakedefine SCALFMM_BLAS_NOCHANGE
////////////////////////////////////////////////////////
// FFT
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_FFT
#cmakedefine SCALFMM_USE_MKL_AS_FFTW
#cmakedefine SCALFMM_USE_ESSL_AS_FFTW
//////////////////////////////////////////////////////
// MPI
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_MPI
///////////////////////////////////////////////////////
// Memory trace
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_MEM_STATS
///////////////////////////////////////////////////////
// CUDA
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_CUDA
///////////////////////////////////////////////////////
// OPENCL
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_OPENCL
///////////////////////////////////////////////////////
// STARPU
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_USE_STARPU
#cmakedefine SCALFMM_DISABLE_NATIVE_OMP4
///////////////////////////////////////////////////////
// To control starpu config
///////////////////////////////////////////////////////
#cmakedefine SCALFMM_STARPU_USE_COMMUTE
#cmakedefine SCALFMM_STARPU_USE_REDUX
#cmakedefine SCALFMM_STARPU_USE_PRIO
#cmakedefine SCALFMM_STARPU_FORCE_NO_SCHEDULER
#cmakedefine SCALFMM_USE_STARPU_EXTRACT
#endif // CONFIG_H
// --------------------------------
// See LICENCE file at project root
// File : particle.hpp
// --------------------------------
#ifndef SCALFMM_CONTAINER_PARTICLE_HPP
#define SCALFMM_CONTAINER_PARTICLE_HPP
#include <inria/integer_sequence.hpp>
#include <scalfmm/container/point.hpp>
#include <scalfmm/meta/type_pack.hpp>
/**
* \brief Multi-purpose particle implementation
*
* This template implementation of a particle allows simple reuse for several
* use cases. The aim it to provide an interface that is compatible with the
* rest of ScalFMM. It is mainly intended to be used as an interface for the
* particle containers.
*
* The Types parameter pack can accept any type that is to be considered as a
* particle attribute. You can also specify scalfmm::pack type to factorise
* several types.
*
* In the following example, the two specialisations of the class will give the
* same final structure.
*
* ```
* using FReal = double;
* static constexpr std::size_t dimension = 3;
*
* particle<FReal, dimension, int, float, float, float, float>;
* particle<FReal, dimension, int, scalfmm::pack<4, float> >;
* ```
*
* The base of these two classes is
* ```
* std::tuple<double, double, double, int, float, float, float, float>;
* ```
*
* \warning Although the classes will have the same final layout, C++ considers
* these two classes to be different !
*
* ##### Example
*
* ```
* // Define a 3D particle with an int attribute
* using Particle = particle<double, 3, int>;
*
* Particle p;
* p.get<>
* ```
*
*
* \tparam FReal Floating point type
* \tparam dimension Space dimension count
* \tparam Types Attributes type list
*
*/
namespace scalfmm::container
{
template<typename ValueType, std::size_t Dim = 3, typename... Types>
class particle : public scalfmm::meta::pack_expand_tuple< scalfmm::meta::pack<Dim, ValueType>, Types... >
{
public:
/// Storage class : std::tuple<FReal,...(dimension times), Types...>
using tuple_data_t = scalfmm::meta::
pack_expand_tuple< scalfmm::meta::pack<Dim, ValueType>, Types... >;
/// Expand Types list
using types_tuple_t = scalfmm::meta::
pack_expand_tuple< Types... >;
/// Expand std::enable_if if possible
template<bool value>
using sfinae_check = typename std::enable_if<value, bool>::type;
/**
* \brief Check parameter pack size vs. attribute + dimension count at
* compile time
*
* \return true if the parameter pack is shorter than the attribute list
* size + the dimension
*/
template<typename... Ts>
constexpr static bool correct_attribute_list_size()
{
return std::tuple_size<tuple_data_t>::value >= sizeof...(Ts) + dimension;
}
/// 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;
/// Floating point type
using value_type = ValueType;
/// Position type, required by the FVariadicParticleContainer
using position_t = point<value_type, dimension>;
/// Default constructor
particle() = default;
/// Default copy constructor
particle(const particle&) = default;
/// Default copy operator
particle& operator=(const particle&) = default;
/// Default move constructor
particle(particle&&) = default;
/// Default move operator
particle& operator=(particle&&) = default;
/**
* \brief Constructor from position and types
*
* \tparam Ts Attributes parameter pack; if Ts is too long, this will not
* compile.
*
* \param pos Particle position
* \param ts Attributes
*
* \warning There may be less attributes than defined by the particle, in
* that case the missing ones are zero constructed.
*/
template<typename... Ts, sfinae_check<correct_attribute_list_size<Ts...>()> = 0>
particle(const position_t& pos, Ts&&... ts)
: particle(inria::make_index_sequence<dimension>(),
inria::make_index_sequence<std::tuple_size<tuple_data_t>::value-sizeof...(ts)-dimension>(),
pos,
std::forward<Ts>(ts)...
)
{}
/// Constructor from tuple equivalent to #tuple_data_t
template<typename... Ts>
particle(const std::tuple<Ts...>& ts)
: tuple_data_t(ts)
{}
/**
* \brief Position getter
*
* The position is stored in the #_data tuple, to extract it we need to
* recreate a position_t object. This is done by the position_impl() method.
*
* \return A new position_t object
*/
position_t position() const
{
return position_impl(inria::make_index_sequence<dimension>());
}
/**
* \brief Position setter
*
* \parma pos The new position
*/
void set_position(const position_t& pos)
{
return set_position_impl(pos, inria::make_index_sequence<dimension>());
}
/**
* \brief Get a reference to the Ith attribute
*
* \tparam I Index of the attribute to get
*/
template<std::size_t I>
auto attribute() -> decltype(std::get<dimension+I>(*this))
{
return std::get<dimension+I>(*this);
}
/**
* \brief Get a const reference to the Ith attribute
*
* \tparam I Index of the attribute to get
*/
template<std::size_t I>
auto attribute() const -> decltype(std::get<dimension+I>(*this))
{
return std::get<dimension+I>(*this);
}
/**
* \brief Get a tuple filled with copies of the attributes
*/
types_tuple_t attributes() const
{
return attributes_impl(
inria::make_index_sequence<std::tuple_size<types_tuple_t>::value>());
}
/**
* \brief Convert particle to a tuple
*/
tuple_data_t& as_tuple()
{
return *this;
}
/**
* \brief Convert particle to a tuple
*/
const tuple_data_t& as_tuple() const
{
return *this;
}
private:
/**
* \brief Contructor implementation
*
* Builds the particle, zero constructing missing arguuments.
*/
template<typename... Ts, std::size_t... Is, std::size_t... Js>
particle(inria::index_sequence<Is...>, inria::index_sequence<Js...>,
const position_t& pos, Ts&&... ts)
: tuple_data_t(pos[Is]..., std::forward<Ts>(ts)...,
typename std::tuple_element<dimension+Js+sizeof...(ts), tuple_data_t>::type(0)...)
{
//static_assert(sizeof...(Ts) == NbAttributes, "Parameter count is incorrect");
}
/**
* \brief #position method implementation
*
* \tparam Is Index sequence of the position elements in the tuple, deduced
* using arguments
*/
template<std::size_t... Is>
position_t position_impl(inria::index_sequence<Is...>) const
{
return position_t(std::get<Is>(this->as_tuple())...);
}
/**
* \brief #setPosition method implementation
*
* \tparam Is Index sequence of the position elements in the tuple, deduced
* using arguments
*
* \param pos new position
*/
template<std::size_t... Is>
void set_position_impl(const position_t& pos, inria::index_sequence<Is...>)
{
auto l = {std::get<Is>(this->as_tuple()) = pos[Is] ...};
(void)l;
}
/**
* \brief #attributes method implementation
*
* \tparam Is Index sequence of the attributes in the tuple, deduced using
* arguments
*
* \param pos new position
*/
template<std::size_t... Is>
types_tuple_t attributes_impl(inria::index_sequence<Is...>) const
{
return types_tuple_t(std::get<Is+dimension>(*this)...);
}
};
}
#endif // SCALFMM_CONTAINER_PARTICLE_HPP
// --------------------------------
// See LICENCE file at project root
// File : particle_container.hpp
// --------------------------------
#ifndef SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
#define SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
#include <vector>
#include <tuple>
#include <xsimd/config/xsimd_align.hpp>
#include <scalfmm/memory/aligned_allocator.hpp>
#include <scalfmm/container/variadic_container.hpp>
#include <scalfmm/meta/traits.hpp>
namespace scalfmm::container
{
template<typename Particle>
class particle_container
: public variadic_container
< memory::aligned_allocator< XSIMD_DEFAULT_ALIGNMENT
, typename Particle::tuple_data_t
>
, typename Particle::tuple_data_t
>
{
public:
// base type : concatenating the particle tuple with the indexes needed
// in order to allocate the vector
using value_type = typename Particle::tuple_data_t;
using allocator_type = typename memory::aligned_allocator< XSIMD_DEFAULT_ALIGNMENT, value_type>;
using base_type = variadic_container< allocator_type, value_type>;
using size_type = typename base_type::size_type;
using difference_type = typename base_type::difference_type;
using reference = typename base_type::reference_tuple;
using const_reference = typename base_type::const_reference_tuple;
using pointer = typename base_type::pointer_tuple;
using const_pointer = typename base_type::const_pointer_tuple;
using iterator = typename base_type::iterator;
using const_iterator = typename base_type::const_iterator;
using particle_type = Particle;
//Forwardind constructors
using base_type::base_type;
};
}
#endif // SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
// See LICENCE file at project root
//
#ifndef SCALFMM_CONTAINER_POINT_HPP
#define SCALFMM_CONTAINER_POINT_HPP
#include <array>
#include <iterator>
#include <ostream>
#include <type_traits>
#include <cmath>
#include <initializer_list>
#include <scalfmm/meta/traits.hpp>
#include <scalfmm/meta/const_functions.hpp>
namespace scalfmm::container
{
template <typename Arithmetic, std::size_t Dim = 3, typename Enable = void>
struct point
{
static_assert(std::is_arithmetic<Arithmetic>::value, "Point's inner type should be arithmetic!");
};
template <typename Arithmetic, std::size_t Dim>
struct point<Arithmetic, Dim, typename std::enable_if<std::is_arithmetic<Arithmetic>::value>::type>
: public std::array<Arithmetic,Dim>
{
public:
using base_type = std::array<Arithmetic,Dim>;
/// Floating number type
using value_type = Arithmetic;
/// Dimension type
using dimension_type = std::size_t;
/// Space dimension count
constexpr static const std::size_t dimension = Dim;
using base_type::base_type;
template<typename... Types>
point(Types&&... e) : base_type{{std::forward<Types>(e)...}} {}
/** Addition assignment operator */
inline point& operator+=(const point& other)
{
for (int i=0; auto a : other)
{
(*this)[i++] += a;
}
return *this;
}
/** Addition operator */
inline friend point& operator+(point& lhs, const point& rhs)
{
lhs += rhs;
return lhs;
}
inline friend point& operator+(const point lhs, const point rhs)
{
point tmp;
for(int i=0; auto a : tmp) tmp = lhs[i++] + rhs[i++];
return tmp;
}
/** Scalar assignment addition */
inline point& operator+=(const value_type val)
{
for (auto& a : *this)
{
a += val;
}
return *this;
}
/** Scalar addition */
inline friend point& operator+(point& lhs, const value_type val)
{
lhs += val;
return lhs;
}
/** Addition assignment operator */
inline point& operator-=(const point& other)
{
for (int i=0; auto a : other)
{
(*this)[i++] -= a;
}
return *this;
}
/** Substraction operator */
inline friend point& operator-(point& lhs, const point& rhs)
{
lhs -= rhs;
return lhs;
}
inline friend point& operator-(const point lhs, const point rhs)
{
point tmp;
for(int i=0; auto a : tmp) tmp = lhs[i++] - rhs[i++];
return tmp;
}
/** Scalar assignment substraction*/
inline point& operator-=(const value_type val)
{
for (auto& a : *this)
{
a -= val;
}
return *this;
}
/** Scalar substraction */
inline friend point& operator-(point& lhs, const value_type val)
{
lhs -= val;
return lhs;
}
/** Right scalar multiplication assignment operator */
inline point& operator*=(const value_type &val)
{
for (auto& a : *this)
{
a *= val;
}
return *this;
}
/** Right scalar multiplication operator */
inline friend point operator*(point lhs, const value_type val)
{
lhs *= val;
return lhs;
}
/** Left scalar multiplication operator */
inline friend point operator*(const value_type val, point rhs)
{
rhs *= val;
return rhs;
}
///** Data to data division assignment */
inline point& operator/=(const point& other)
{
for (std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] /= other[i];
}
return *this;
}
/** Data to data division operator */
inline friend point operator/(point lhs, const point& rhs)
{
lhs /= rhs;
return lhs;
}
/** Right scalar division assignment operator */
inline point& operator/=(const value_type val)
{
for (auto& a : *this)
{
a /= val;
}
return *this;
}
/** Right scalar division operator */
inline friend point operator/(point lhs, const value_type val)
{
lhs /= val;
return lhs;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
-> std::enable_if_t<meta::is_float<typename T::value_type>::value,bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it) {
if(! meta::feq(*lhs_it, *rhs_it, T{1e-7})) {
return false;
}
}
return true;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
-> std::enable_if_t<meta::is_double<typename T::value_type>::value,bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it) {
if(! meta::feq(*lhs_it, *rhs_it, T{1e-13})) {
return false;
}
}
return true;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<value_type, Dim>& rhs)
-> std::enable_if_t<std::is_integral<typename T::value_type>::value,bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; i++, ++lhs_it, ++rhs_it) {
if( *lhs_it != *rhs_it) {
return false;
}
}
return true;
}
/** Non equality test operator */
inline friend auto operator!=(const point& lhs, const point& rhs)
{
return !(lhs == rhs);
}
/** Formated output stream operator */
inline friend std::ostream &operator<<(std::ostream &os, const point<value_type, Dim> &pos)
{
os << "[";
for (auto it = pos.begin(); it != pos.end() - 1; it++)
os << *it << ", ";
os << pos.back() << "]";
return os;
}
/** Formated input stream operator */
inline friend std::istream &operator>>(std::istream &is, point<value_type, Dim> &pos)
{
char c;
is >> c; // opening '['
for (std::size_t i = 0; i + 1 < Dim; ++i)
{
is >> pos.data()[i];
is >> c; // get coma ','
}
is >> pos.data()[Dim - 1];
is >> c; // closing ']'
return is;
}
};
} // end of namespace scalfmm::container
#endif
// --------------------------------
// See LICENCE file at project root
// File : simple_variadic.cpp
// --------------------------------
#ifndef SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
#define SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
#include <vector>
#include <tuple>
#include <type_traits>
#include <utility>
#include <xsimd/config/xsimd_config.hpp>
namespace scalfmm::container
{
template<typename Types, typename Indices>
class variad_impl;
//template<typename Allocator = memory::aligned_allocator<XSIMD_DEFAULT_ALIGNMENT, Ts>>
template<typename... Ts, std::size_t... Indices>
class variad_impl<std::tuple<Ts...>, std::integer_sequence<std::size_t, Indices...>> : public std::tuple<std::vector<Ts, XSIMD_DEFAULT_ALLOCATOR(Ts)>...>
{
private:
//Discard fold expression results
struct noop_t
{
template<typename... Types>
noop_t(const Types&... ts) {}
};
public:
using base_type = std::tuple< std::vector<Ts, XSIMD_DEFAULT_ALLOCATOR(Ts)>... >;
using value_type = std::tuple< Ts... >;
using allocator_type = std::tuple<XSIMD_DEFAULT_ALLOCATOR(Ts)...>;
using size_type = std::size_t;
using difference_type = std::ptrdiff_t;
using reference = std::tuple<Ts&...>&;
using const_reference = const std::tuple<const Ts&...>&;
using pointer = std::tuple<Ts*...>;
using const_pointer = const std::tuple<const Ts*...>;
private:
allocator_type m_allocator{};
public:
variad_impl() = default;
variad_impl(const variad_impl&) = default;
variad_impl(variad_impl&&) noexcept = default;
variad_impl& operator=(const variad_impl&) = default;
variad_impl& operator=(variad_impl&&) noexcept = default;
~variad_impl() = default;
// TODO
//explicit variad( const allocator_type& alloc){}
explicit variad_impl( size_type count
, const value_type& value
, const allocator_type& alloc = allocator_type() )
: m_allocator(alloc)
{
reserve(count);
for(std::size_t i=0; i<count; ++i)
{
push_back(value);
}
}
explicit variad_impl(size_type count)
{
reserve(count);
for(std::size_t i=0; i<count; ++i)
{
push_back(value_type{});
}
}
inline void reserve(size_type size)
{
noop_t{(std::get<Indices>(*this).reserve(size),0) ...};
}
inline void push_back(const value_type& value)
{
noop_t{(std::get<Indices>(*this).push_back(std::get<Indices>(value)),0) ...};
}
};
template<typename... Ts>
struct variad : public variad_impl<std::tuple<Ts...>, std::make_index_sequence<sizeof...(Ts)>>
{
using base_type = variad_impl<std::tuple<Ts...>, std::make_index_sequence<sizeof...(Ts)>>;
using base_type::base_type;
};
}
#endif // SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
This diff is collapsed.
// See LICENCE file at project root
//
#ifndef SCALFMM_FUNCTIONAL_UTILS_HPP
#define SCALFMM_FUNCTIONAL_UTILS_HPP
#include <cmath>
#include <scalfmm/meta/traits.hpp>
namespace scalfmm::utils
{
// Norm2 on a range
template<typename T>
inline auto norm2(const T& range)
-> std::enable_if_t< meta::has_value_type<T>::value
&& meta::has_range_interface<T>::value
, typename T::value_type>
{
typename T::value_type square_sum{0};
for (auto a : range)
{
square_sum += a*a;
}
return square_sum;
}
//Norm
template<typename T>
inline auto norm(const T& range)
-> std::enable_if_t< meta::has_value_type<T>::value
&& meta::has_range_interface<T>::value
, typename T::value_type>
{ return std::sqrt(norm2(range)); }
}
#endif // SCALFMM_FUNCTIONAL_UTILC_HPP
// --------------------------------
// See LICENCE file at project root
// File : interpolation/builders.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_BUILDERS_HPP
#define SCALFMM_INTERPOLATION_BUILDERS_HPP
#include <scalfmm/meta/const_functions.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor-blas/xlinalg.hpp>
namespace scalfmm::interpolation
{
template<std::size_t, typename T>
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)
{
return get_generator(xt::linspace(std::size_t{0},order-1, order), std::make_index_sequence<dim>{});
}
// TODO for nodes
//auto lin = xt::stack(xt::xtuple(xt::linspace(std::size_t{0},order-1,order)));
//xt::xarray<std::size_t> saved;
//for(std::size_t i=0; i<dimension; ++i)
//{
// auto r = xt::stack(xt::xtuple(xt::ones<std::size_t>({meta::pow(order,(dimension-1-i))})));
// auto tmp = xt::linalg::kron(lin,r);
// auto l = xt::stack(xt::xtuple(xt::ones<std::size_t>({meta::pow(order,i)})));
// auto reshaped = xt::linalg::kron(l,tmp);
// reshaped.reshape({1,-1});
// auto res = xt::stack(xt::xtuple(reshaped,saved));
// saved = res;
//}
// std::cout << "ids: " << saved << '\n';
template<std::size_t dim> struct nodes{ static_assert(dim<4, "Dimension for interpolation node not supported."); };
template<>
struct nodes<1>
{
inline static auto get(std::size_t order)
{
return xt::linspace(std::size_t{0}, order-1, order);
}
};
template<>
struct nodes<2>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,2)};
const xt::xarray<std::size_t>::shape_type shape{2,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/order)%order;
}
return ns;
}
};
template<>
struct nodes<3>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,3)};
const xt::xarray<std::size_t>::shape_type shape{3,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/order)%order;
ns(2,i) = i/(order*order);
}
return ns;
}
};
template<>
struct nodes<4>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,4)};
const xt::xarray<std::size_t>::shape_type shape{4,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/(order))%order;
ns(2,i) = (i/(order*order))%order;
ns(3,i) = i/(order*order*order);
}
return ns;
}
};
// Builders for relative center
}
#endif // SCALFMM_INTERPOLATION_BUILDERS_HPP
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebyshev_interpolator.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#include <array>
#include <cmath>
#include <type_traits>
#include <scalfmm/interpolation/chebyshev/chebyshev_root_factory.hpp>
#include <scalfmm/meta/const_functions.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xio.hpp>
namespace scalfmm::interpolation
{
template<std::size_t, typename T>
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 node_ids(const std::size_t order)
{
return get_generator(xt::linspace(std::size_t{0},order-1, order), std::make_index_sequence<dim>{});
}
template<typename ValueType, std::size_t Dimension>
struct chebyshev_interpolator
{
public:
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
chebyshev_interpolator() = default;
chebyshev_interpolator(chebyshev_interpolator const &) = default;
chebyshev_interpolator(chebyshev_interpolator &&) = default;
chebyshev_interpolator& operator=(chebyshev_interpolator const&) = default;
chebyshev_interpolator& operator=(chebyshev_interpolator &&) = default;
explicit chebyshev_interpolator(const std::size_t order, std::size_t level)
: m_order(order), m_level(level), m_nnodes(meta::pow(order,dimension))
{
//constexpr unsigned int size_of_coef = std::pow(m_order, dimension);
//for(auto a : m_shape) { a = size_of_coef; }
//m_permutations.reshape(m_shape);
m_ts.resize({m_order,m_order});
//m_flat_ts.resize({(m_order)*m_order});
// Initializing roots
chebyshev_root_factory<value_type> cheb_fact{};
auto cheb_roots = cheb_fact.get(m_order);
// Calculating Ts
for( std::size_t o = 1; o < m_order; ++o )
for( std::size_t j = 0; j < m_order; ++j )
{ m_ts(o,j) = chebyshev::ts(o,cheb_roots(j)); }
auto lin = node_ids<dimension>(order);
std::cout << lin.size() << '\n';
std::cout << lin.dimension() << '\n';
std::cout << lin << '\n';
}
private:
const std::size_t m_order{};
const std::size_t m_level{};
const std::size_t m_nnodes{};
const std::array<std::size_t, dimension> m_shape{};
xt::xtensor<unsigned int, dimension> m_permutations{m_shape};
xt::xtensor<value_type, 2> m_ts{};
xt::xtensor<std::size_t, dimension> m_node_ids{};
//xt::xarray<ValueType,Dimension>
};
}
#endif // SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebychev/chebyshev_root_factory.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_ROOT_FACTORY_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_ROOT_FACTORY_HPP
#include <vector>
#include <exception>
#include <cmath>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
namespace scalfmm::interpolation
{
namespace chebyshev
{
template<typename T>
constexpr T ts(const unsigned int n, T x)
{
assert(std::abs(x)-1.<10.*std::numeric_limits<T>::epsilon());
if (std::abs(x)>1.)
{
x = (x > T( 1.) ? T( 1.) : x);
x = (x < T(-1.) ? T(-1.) : x);
}
return T(cos(n * acos(x)));
}
template<typename T>
constexpr T us(const unsigned int n, T x)
{
assert(std::abs(x)-1.<10.*std::numeric_limits<T>::epsilon());
if (std::abs(x)>1.)
{
x = (x > T( 1.) ? T( 1.) : x);
x = (x < T(-1.) ? T(-1.) : x);
}
return T(n * (sin(n * acos(x))) / sqrt(1 - x*x));
}
} // chebyshev
template<typename ValueType, typename Enable = void>
struct chebyshev_root_factory
{
static_assert( std::is_floating_point<ValueType>::value
, "Can't construct root elements with non floating types!");
};
template<typename ValueType>
struct chebyshev_root_factory<ValueType, std::enable_if_t<std::is_floating_point<ValueType>::value>>
{
using value_type = ValueType;
using array_type = xt::xtensor<value_type,1>;
auto get(const unsigned int order) const
{
if(order > 14 || order < 0) throw std::out_of_range("Order out of range, max order is 14.");
return xt::view(m_roots, xt::range(indexes[order],indexes[order]+order));
}
const std::array<std::ptrdiff_t,16> indexes{-1,-1,0,2,5,9,14,20,27,35,44,54,65,77,90,104};
const array_type m_roots{
value_type(-0.707106781186548),
value_type( 0.707106781186547),
value_type(-8.66025403784439e-01),
value_type( 0.0),
value_type( 8.66025403784438e-01),
value_type(-0.923879532511287),
value_type(-0.382683432365090),
value_type( 0.382683432365090),
value_type( 0.923879532511287),
value_type(-9.51056516295154e-01),
value_type(-5.87785252292473e-01),
value_type( 0.0),
value_type( 5.87785252292473e-01),
value_type( 9.51056516295154e-01),
value_type(-0.965925826289068),
value_type(-0.707106781186548),
value_type(-0.258819045102521),
value_type( 0.258819045102521),
value_type( 0.707106781186547),
value_type( 0.965925826289068),
value_type(-9.74927912181824e-01),
value_type(-7.81831482468030e-01),
value_type(-4.33883739117558e-01),
value_type( 0.0 ),
value_type( 4.33883739117558e-01),
value_type( 7.81831482468030e-01),
value_type( 9.74927912181824e-01),
value_type(-0.980785280403230),
value_type(-0.831469612302545),
value_type(-0.555570233019602),
value_type(-0.195090322016128),
value_type( 0.195090322016128),
value_type( 0.555570233019602),
value_type( 0.831469612302545),
value_type( 0.980785280403230),
value_type(-9.84807753012208e-01),
value_type(-8.66025403784439e-01),
value_type(-6.42787609686539e-01),
value_type(-3.42020143325669e-01),
value_type( 0.0 ),
value_type( 3.42020143325669e-01),
value_type( 6.42787609686539e-01),
value_type( 8.66025403784438e-01),
value_type( 9.84807753012208e-01),
value_type(-0.987688340595138),
value_type(-0.891006524188368),
value_type(-0.707106781186548),
value_type(-0.453990499739547),
value_type(-0.156434465040231),
value_type( 0.156434465040231),
value_type( 0.453990499739547),
value_type( 0.707106781186547),
value_type( 0.891006524188368),
value_type( 0.987688340595138),
value_type(-9.89821441880933e-01),
value_type(-9.09631995354518e-01),
value_type(-7.55749574354258e-01),
value_type(-5.40640817455598e-01),
value_type(-2.81732556841430e-01),
value_type( 0.0 ),
value_type( 2.81732556841430e-01),
value_type( 5.40640817455597e-01),
value_type( 7.55749574354258e-01),
value_type( 9.09631995354518e-01),
value_type( 9.89821441880933e-01),
value_type(-0.991444861373810),
value_type(-0.923879532511287),
value_type(-0.793353340291235),
value_type(-0.608761429008721),