From f50a9732bc40b4da1ec1d40552610f937fb5b3fc Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr> Date: Tue, 16 Feb 2016 18:12:12 +0100 Subject: [PATCH] Ajout des sources et tests c++ et modification du CMake. --- CMake/FindFFTW.cmake | 9 +- CMake/FindPythonFull.cmake | 5 +- CMake/GoogleTestHelper.cmake | 47 ++ HySoP/CMakeLists.txt | 15 +- HySoP/setup.py.in | 2 +- HySoP/src/CMakeLists.txt | 49 +- HySoP/src/hysop++/main/diffSolver.cpp | 216 +++++ HySoP/src/hysop++/main/poissonSolver.cpp | 182 +++++ HySoP/src/hysop++/src/data/accumulatorIndex.h | 65 ++ HySoP/src/hysop++/src/data/index.h | 125 +++ .../hysop++/src/data/memory/fftwAllocator.h | 53 ++ .../src/data/memory/minimalAllocator.h | 52 ++ HySoP/src/hysop++/src/data/multi_array/< | 352 +++++++++ .../data/multi_array/const_multi_array_ref.h | 53 ++ .../data/multi_array/const_multi_array_view.h | 53 ++ .../src/data/multi_array/multi_array.h | 51 ++ .../src/data/multi_array/multi_array_clean.h | 47 ++ .../data/multi_array/multi_array_defines.h | 374 +++++++++ .../src/data/multi_array/multi_array_ext.h | 247 ++++++ .../src/data/multi_array/multi_array_impl.h | 192 +++++ .../src/data/multi_array/multi_array_ref.h | 70 ++ .../src/data/multi_array/multi_array_view.h | 66 ++ HySoP/src/hysop++/src/data/multi_array/tags | 6 + HySoP/src/hysop++/src/detail/index_seq.h | 43 + HySoP/src/hysop++/src/domain/boundary.cpp | 26 + HySoP/src/hysop++/src/domain/boundary.h | 23 + HySoP/src/hysop++/src/domain/domain.h | 154 ++++ .../hysop++/src/domain/domainConfiguration.h | 98 +++ HySoP/src/hysop++/src/fft/extension.cpp | 24 + HySoP/src/hysop++/src/fft/extension.h | 23 + .../hysop++/src/fft/fftDomainConfiguration.h | 114 +++ HySoP/src/hysop++/src/fft/fftw3.h | 440 +++++++++++ HySoP/src/hysop++/src/fft/fftwComplex.h | 69 ++ HySoP/src/hysop++/src/fft/planner.h | 745 ++++++++++++++++++ HySoP/src/hysop++/src/fft/transform.cpp | 14 + HySoP/src/hysop++/src/fft/transform.h | 188 +++++ HySoP/src/hysop++/src/maths/polynomial.h | 427 ++++++++++ HySoP/src/hysop++/src/solver/diffSolver.h | 29 + HySoP/src/hysop++/src/solver/fftDiffSolver.h | 191 +++++ .../src/hysop++/src/solver/fftPoissonSolver.h | 205 +++++ HySoP/src/hysop++/src/solver/poissonSolver.h | 27 + HySoP/src/hysop++/src/utils/constants.h | 16 + HySoP/src/hysop++/src/utils/default.h | 20 + HySoP/src/hysop++/src/utils/defines.h | 46 ++ HySoP/src/hysop++/src/utils/types.h | 59 ++ HySoP/src/hysop++/src/utils/utils.cpp | 17 + HySoP/src/hysop++/src/utils/utils.h | 196 +++++ HySoP/src/hysop++/tests/CMakeLists.txt | 7 + .../tests/testDiffSolver/CMakeLists.txt | 12 + .../src/hysop++/tests/testDiffSolver/main.cpp | 10 + .../tests/testDiffSolver/testDiffSolver.cpp | 221 ++++++ .../tests/testDiffSolver/testDiffSolver.h | 14 + .../hysop++/tests/testPlanner/CMakeLists.txt | 12 + HySoP/src/hysop++/tests/testPlanner/main.cpp | 10 + .../hysop++/tests/testPlanner/testPlanner.cpp | 180 +++++ .../hysop++/tests/testPlanner/testPlanner.h | 20 + .../tests/testPoissonSolver/CMakeLists.txt | 12 + .../hysop++/tests/testPoissonSolver/main.cpp | 10 + .../testPoissonSolver/testPoissonSolver.cpp | 186 +++++ .../testPoissonSolver/testPoissonSolver.h | 14 + .../hysop++/tests/testPolynoms/CMakeLists.txt | 12 + HySoP/src/hysop++/tests/testPolynoms/main.cpp | 10 + .../tests/testPolynoms/testPolynoms.cpp | 103 +++ .../hysop++/tests/testPolynoms/testPolynoms.h | 13 + 64 files changed, 6364 insertions(+), 7 deletions(-) create mode 100644 CMake/GoogleTestHelper.cmake create mode 100644 HySoP/src/hysop++/main/diffSolver.cpp create mode 100644 HySoP/src/hysop++/main/poissonSolver.cpp create mode 100644 HySoP/src/hysop++/src/data/accumulatorIndex.h create mode 100644 HySoP/src/hysop++/src/data/index.h create mode 100644 HySoP/src/hysop++/src/data/memory/fftwAllocator.h create mode 100644 HySoP/src/hysop++/src/data/memory/minimalAllocator.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/< create mode 100644 HySoP/src/hysop++/src/data/multi_array/const_multi_array_ref.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/const_multi_array_view.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_clean.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_ext.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_ref.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/multi_array_view.h create mode 100644 HySoP/src/hysop++/src/data/multi_array/tags create mode 100644 HySoP/src/hysop++/src/detail/index_seq.h create mode 100644 HySoP/src/hysop++/src/domain/boundary.cpp create mode 100644 HySoP/src/hysop++/src/domain/boundary.h create mode 100644 HySoP/src/hysop++/src/domain/domain.h create mode 100644 HySoP/src/hysop++/src/domain/domainConfiguration.h create mode 100644 HySoP/src/hysop++/src/fft/extension.cpp create mode 100644 HySoP/src/hysop++/src/fft/extension.h create mode 100644 HySoP/src/hysop++/src/fft/fftDomainConfiguration.h create mode 100644 HySoP/src/hysop++/src/fft/fftw3.h create mode 100644 HySoP/src/hysop++/src/fft/fftwComplex.h create mode 100644 HySoP/src/hysop++/src/fft/planner.h create mode 100644 HySoP/src/hysop++/src/fft/transform.cpp create mode 100644 HySoP/src/hysop++/src/fft/transform.h create mode 100644 HySoP/src/hysop++/src/maths/polynomial.h create mode 100644 HySoP/src/hysop++/src/solver/diffSolver.h create mode 100644 HySoP/src/hysop++/src/solver/fftDiffSolver.h create mode 100644 HySoP/src/hysop++/src/solver/fftPoissonSolver.h create mode 100644 HySoP/src/hysop++/src/solver/poissonSolver.h create mode 100644 HySoP/src/hysop++/src/utils/constants.h create mode 100644 HySoP/src/hysop++/src/utils/default.h create mode 100644 HySoP/src/hysop++/src/utils/defines.h create mode 100644 HySoP/src/hysop++/src/utils/types.h create mode 100644 HySoP/src/hysop++/src/utils/utils.cpp create mode 100644 HySoP/src/hysop++/src/utils/utils.h create mode 100644 HySoP/src/hysop++/tests/CMakeLists.txt create mode 100644 HySoP/src/hysop++/tests/testDiffSolver/CMakeLists.txt create mode 100644 HySoP/src/hysop++/tests/testDiffSolver/main.cpp create mode 100644 HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp create mode 100644 HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.h create mode 100644 HySoP/src/hysop++/tests/testPlanner/CMakeLists.txt create mode 100644 HySoP/src/hysop++/tests/testPlanner/main.cpp create mode 100644 HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp create mode 100644 HySoP/src/hysop++/tests/testPlanner/testPlanner.h create mode 100644 HySoP/src/hysop++/tests/testPoissonSolver/CMakeLists.txt create mode 100644 HySoP/src/hysop++/tests/testPoissonSolver/main.cpp create mode 100644 HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp create mode 100644 HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.h create mode 100644 HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt create mode 100644 HySoP/src/hysop++/tests/testPolynoms/main.cpp create mode 100644 HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp create mode 100644 HySoP/src/hysop++/tests/testPolynoms/testPolynoms.h diff --git a/CMake/FindFFTW.cmake b/CMake/FindFFTW.cmake index bec57b85a..eb00a9802 100644 --- a/CMake/FindFFTW.cmake +++ b/CMake/FindFFTW.cmake @@ -78,6 +78,13 @@ find_library(FFTWFloat_MPI_LIBRARY NAMES fftw3f_mpi ) + +#other types for the CXX library +find_library (FFTWD_LIB "fftw3") +find_library (FFTWL_LIB "fftw3l") +find_library (FFTWQ_LIB "fftw3q") +set(FFTW_OTHER_LIBRARIES ${FFTW_LIB} ${FFTWF_LIB} ${FFTWL_LIB} ${FFTWQ_LIB}) + set(FFTW_PROCESS_INCLUDES FFTW_INCLUDE_DIR) -set(FFTW_PROCESS_LIBS FFTW_LIBRARY FFTWFloat_LIBRARY FFTW_MPI_LIBRARY FFTWFloat_MPI_LIBRARY) +set(FFTW_PROCESS_LIBS FFTW_LIBRARY FFTWFloat_LIBRARY FFTW_MPI_LIBRARY FFTWFloat_MPI_LIBRARY FFTW_OTHER_LIBRARIES) libfind_process(FFTW) diff --git a/CMake/FindPythonFull.cmake b/CMake/FindPythonFull.cmake index 604086dd8..6447a4af8 100644 --- a/CMake/FindPythonFull.cmake +++ b/CMake/FindPythonFull.cmake @@ -59,7 +59,10 @@ else() NAMES python${PYTHON_VERSION_NO_DOTS} python${PYTHON_VERSION} NO_DEFAULT_PATH - HINTS ${PYTHON_PREFIX} ${PYTHON_PREFIX}/lib/python${PYTHON_VERSION}/config + HINTS + ${PYTHON_PREFIX} + ${PYTHON_PREFIX}/lib/python${PYTHON_VERSION}/config + ${PYTHON_PREFIX}/lib/python${PYTHON_VERSION}/config-x86_64-linux-gnu PATH_SUFFIXES lib ) # find_library(PYTHON_LIBRARY diff --git a/CMake/GoogleTestHelper.cmake b/CMake/GoogleTestHelper.cmake new file mode 100644 index 000000000..44074c6dc --- /dev/null +++ b/CMake/GoogleTestHelper.cmake @@ -0,0 +1,47 @@ + +find_package(Threads REQUIRED) +include(ExternalProject) + +# download and install GoogleTest +ExternalProject_Add( + gtest + URL http://googletest.googlecode.com/files/gtest-1.7.0.zip + PREFIX ${CMAKE_CURRENT_BINARY_DIR}/gtest + INSTALL_COMMAND "" +) + +# create a libgtest target to be used as a dependency by test programs +add_library(libgtest IMPORTED STATIC GLOBAL) +add_dependencies(libgtest gtest) + +# set gtest properties +ExternalProject_Get_Property(gtest source_dir binary_dir) +set_target_properties(libgtest PROPERTIES + "IMPORTED_LOCATION" "${binary_dir}/libgtest.a" + "IMPORTED_LINK_INTERFACE_LIBRARIES" "${CMAKE_THREAD_LIBS_INIT}" +) +include_directories(SYSTEM "${source_dir}/include") + +# Download and install GoogleMock +ExternalProject_Add( + gmock + URL http://googlemock.googlecode.com/files/gmock-1.7.0.zip + PREFIX ${CMAKE_CURRENT_BINARY_DIR}/gmock + INSTALL_COMMAND "" +) + +# create a libgmock target to be used as a dependency by test programs +add_library(libgmock IMPORTED STATIC GLOBAL) +add_dependencies(libgmock gmock) + +# set gmock properties +ExternalProject_Get_Property(gmock source_dir binary_dir) +set_target_properties(libgmock PROPERTIES + "IMPORTED_LOCATION" "${binary_dir}/libgmock.a" + "IMPORTED_LINK_INTERFACE_LIBRARIES" "${CMAKE_THREAD_LIBS_INIT}" +) +include_directories(SYSTEM "${source_dir}/include") + +set(GTEST_LIBRARIES libgtest libgmock) + + diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt index 3e5ad4cd9..0491bc9e1 100644 --- a/HySoP/CMakeLists.txt +++ b/HySoP/CMakeLists.txt @@ -4,6 +4,7 @@ # It includes : # - high level python interface to HySoP routines # - HySoP fortran library (with fftw solver and scales interface) +# - HySoP C++ library (poisson solver) # # HySoP depends (optionally) on : # - scales (WITH_SCALES=ON, default) @@ -31,10 +32,12 @@ option(USE_MPI "compile and link HySoP with mpi when this mode is enable. Defaul option(WITH_TESTS "Enable testing. Default = off" ON) option(BUILD_SHARED_LIBS "Enable dynamic library build, default = ON." ON) option(WITH_LIB_FORTRAN "Generate libhysop from fortran files in src, wrapped into hysop.f2py module. Default = ON." ON) +option(WITH_LIB_CXX "Generate libhysop from fortran files in src, wrapped into hysop.cpp2hysop module. Default = ON." ON) option(WITH_SCALES "compile/create scales lib and link it with HySoP. Default = ON." ON) option(WITH_FFTW "Link with fftw library (required for some HySoP solvers), default = ON" ON) option(WITH_GPU "Use of GPU (required for some HySoP solvers), default = ON" ON) option(WITH_MAIN_FORTRAN "Create an executable (test purpose) from fortran sources in src/main, linked with libhysop, default = ON" ON) +option(WITH_MAIN_CXX "Create an executable (test purpose) from cxx sources in src/hysop++/main, linked with libhysop, default = ON" ON) option(DEBUG "Enable debug mode for HySoP (0:disabled, 1:verbose, 2:trace, 3:verbose+trace). Default = 0" 0) option(FULL_TEST "Enable all test options (pep8, mpi ...) - Default = OFF" OFF) option(PROFILE "Enable profiling mode for HySoP. 0:disabled, 1: enabled. Default = 0" 0) @@ -48,8 +51,13 @@ if(NOT WITH_LIB_FORTRAN) set(WITH_MAIN_FORTRAN "OFF") endif() +if(NOT WITH_LIB_CXX) + message(WARNING "You deactivate libhysop (cxx) generation. This will disable the Aitken-Schwarz Poisson solver.") + set(WITH_MAIN_CXX "OFF") +endif() + # We can not run scales or fftw without mpi ... -if(WITH_FFTW OR WITH_SCALES) +if(WITH_FFTW OR WITH_SCALES OR WITH_HYSOP_PLUS_PLUS) set(USE_MPI "ON") endif() @@ -138,8 +146,8 @@ include(HySoPInstallSetup) # Remark : this must be done before add_subdir below, since install process in src needs CMAKE_INSTALL_PREFIX # to be properly set. -# ====== Create non-python (fortran) libraries (fftw and scales interfaces), if required ===== -if(WITH_LIB_FORTRAN) +# ====== Create non-python (fortran, cxx) libraries (fftw and scales interfaces), if required ===== +if(WITH_LIB_FORTRAN OR WITH_LIB_CXX) add_subdirectory(src) endif() @@ -220,6 +228,7 @@ if(VERBOSE_MODE) message(STATUS " Python user include : ${PYTHON_INCLUDE_DIR}") message(STATUS " Install directory : ${CMAKE_INSTALL_PREFIX}") message(STATUS " Fortran compiler : ${CMAKE_Fortran_COMPILER}") + message(STATUS " C++ compiler : ${CMAKE_CXX_COMPILER}") message(STATUS " Sources are in : ${CMAKE_SOURCE_DIR}") message(STATUS " Project uses MPI : ${USE_MPI}") message(STATUS " Project uses Scales : ${WITH_SCALES}") diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in index 56a77e507..6c8cd130b 100644 --- a/HySoP/setup.py.in +++ b/HySoP/setup.py.in @@ -146,7 +146,7 @@ data_files = [] swig_dirs = [os.path.join('@CMAKE_SOURCE_DIR@','swig')] # C++ files and swig interface -cpp_src_dirs = ['src/fftw'] +cpp_src_dirs = ['src/fftw','src/hysop++'] for sd in cpp_src_dirs: swig_dirs.append(os.path.join('@CMAKE_SOURCE_DIR@', sd)) ext = {} diff --git a/HySoP/src/CMakeLists.txt b/HySoP/src/CMakeLists.txt index 92af38cd9..da8d5d0d1 100644 --- a/HySoP/src/CMakeLists.txt +++ b/HySoP/src/CMakeLists.txt @@ -44,7 +44,7 @@ set(EXTS_HDRS *.hpp *.h) # - ${HYSOP_LIBRARY_NAME}_BINARY_DIR : where you have run cmake, i.e. the place for compilation # - ${HYSOP_LIBRARY_NAME}_SOURCE_DIR : where sources (.f and .h and this CMakeLists.txt) are located # Note that because of OutOfSourceBuild, binary_dir and source_dir must be different. -project(${HYSOP_LIBRARY_NAME} Fortran) +project(${HYSOP_LIBRARY_NAME} Fortran CXX) # ============= Search for libraries ============= # We search for libraries HySoP depends on and @@ -63,6 +63,7 @@ set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/Modules) # Add compilation flags: #append_Fortran_FLAGS("-Wall -fPIC -ffree-line-length-none -DBLOCKING_SEND_PLUS -DBLOCKING_SEND") append_Fortran_FLAGS("-Wall -fPIC -ffree-line-length-none") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -W -Wall -Wextra -Wno-unused-variable -Wno-unused-but-set-variable -Wno-unused-parameter -std=c++11") if(USE_MPI) # Find MPI for fortran. @@ -98,6 +99,52 @@ if(WITH_FFTW) set(FFTWLIB ${dirlist} CACHE PATH "fftw libraries dir") endif() +# --- HySoP++ --- +if(WITH_LIB_CXX) + get_filename_component(CXX_DIR "${CMAKE_CURRENT_SOURCE_DIR}/hysop++" ABSOLUTE) + get_filename_component(CXX_MAIN_DIR "${CXX_DIR}/main" ABSOLUTE) + get_filename_component(CXX_TEST_DIR "${CXX_DIR}/tests" ABSOLUTE) + get_filename_component(CXX_SOURCE_DIR "${CXX_DIR}/src" ABSOLUTE) + + find_package(FFTW REQUIRED) + set(CXX_EXT_INCLUDES ${FFTW_INCLUDES}) + set(CXX_EXT_LIBS ${FFTW_LIBRARIES}) + set(CXX_EXT_LIB_DIRS "") + + file(GLOB_RECURSE cxx_header_files ${CXX_SOURCE_DIR}*.h) + file(GLOB_RECURSE cxx_source_files ${CXX_SOURCE_DIR}/*.cpp) + + list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/poissonSolver.cpp") + list(APPEND cxx_executable_sources "${CXX_MAIN_DIR}/diffSolver.cpp") + + set(cxx_source_files_no_main ${cxx_source_files}) + foreach(cxx_main_source ${cxx_executable_sources}) + list(REMOVE_ITEM cxx_source_files_no_main "${cxx_main_source}") + endforeach() + + list(APPEND cxx_local_include_dirs "${CXX_SOURCE_DIR}") + include_directories(${cxx_local_include_dirs}) + include_directories(SYSTEM ${CXX_EXT_INCLUDES}) + #link_directories(${CXX_EXT_LIB_DIRS}) + + set(HYSOP_CXX_LIBRARY "hysop++") + add_library(${HYSOP_CXX_LIBRARY} STATIC ${cxx_source_files_no_main} ${cxx_header_files}) + target_link_libraries(${HYSOP_CXX_LIBRARY} ${CXX_EXT_LIBS}) + + if(WITH_MAIN_CXX) + foreach(cxx_main_source ${cxx_executable_sources}) + get_filename_component(cxx_exec_name "${cxx_main_source}" NAME_WE) + add_executable(${cxx_exec_name} ${cxx_main_source}) + target_link_libraries(${cxx_exec_name} ${HYSOP_CXX_LIBRARY} ${CXX_EXT_LIBS}) + endforeach() + endif() + + if(WITH_TESTS) + add_subdirectory(${CXX_TEST_DIR}) + endif() + +endif() + # ============= Generates HySoPConfig.hpp ============= # The file HYSOP_LIBRARY_NAME_defines.hpp will be generated from config.hpp.cmake; if(EXISTS ${CMAKE_SOURCE_DIR}/config.hpp.cmake) diff --git a/HySoP/src/hysop++/main/diffSolver.cpp b/HySoP/src/hysop++/main/diffSolver.cpp new file mode 100644 index 000000000..0acad6456 --- /dev/null +++ b/HySoP/src/hysop++/main/diffSolver.cpp @@ -0,0 +1,216 @@ + +#include <cstdlib> + +#include "domain/domain.h" +#include "solver/fftDiffSolver.h" +#include "data/multi_array/multi_array.h" +#include "utils/constants.h" +#include "fft/extension.h" + +using namespace hysop; +using namespace hysop::domain; + +static constexpr std::size_t nExtensions = 4 ; +static constexpr std::size_t nExtensionsPair = 7 ; +static constexpr fft::Extension ext[nExtensions] = +{ fft::Extension::NONE, fft::Extension::ODD, fft::Extension::EVEN, fft::Extension::PERIODIC }; +static constexpr std::pair<fft::Extension,fft::Extension> pext[nExtensionsPair] { + std::make_pair(ext[3],ext[3]), //periodic-periodic + std::make_pair(ext[3],ext[3]), //periodic-periodic + std::make_pair(ext[2],ext[1]), //even-odd + std::make_pair(ext[1],ext[2]), //odd-even + std::make_pair(ext[2],ext[2]), //even-even + std::make_pair(ext[1],ext[1]), //odd-odd + std::make_pair(ext[0],ext[0]), //none-none +}; + +static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L }; + +template <typename T> +std::function<T(T)> func(std::size_t k) { + switch(k) { + case 0: return [=](T x) {return cos(T(freqs[0])*x);}; + case 1: return [=](T x) {return sin(T(freqs[1])*x);}; + case 2: return [=](T x) {return cos(T(freqs[2])*x);}; + case 3: return [=](T x) {return sin(T(freqs[3])*x);}; + case 4: return [=](T x) {return cos(T(freqs[4])*x);}; + case 5: return [=](T x) {return sin(T(freqs[5])*x);}; + default: return[=](T x) { return T(1); }; + } +} + +template <typename T> +std::function<T(T)> derivative(std::size_t k, int order) { + bool even = (k%2==0); + std::size_t p, offset; + T sign, coeff; + if(k>5) { + if(order != 0) + throw std::runtime_error("Non zero order !"); + return func<T>(k); + } + else if(even) { /* cos func */ + p = (order%2==0 ? k : k+1); + sign = std::pow(T(-1),(order+1)/2); + coeff = std::pow(freqs[k], order); + } + else { /* sin func */ + p = (order%2==0 ? k : k-1); + sign = std::pow(T(-1),order/2); + coeff = std::pow(freqs[k], order); + } + return [=](T x) { return sign*coeff*(func<T>(p)(x)); }; +} + +template <typename T, std::size_t Dim, bool verbose=false> +void test(std::size_t p_maxOrder, bool includePeriodicBds=false) { + Shape<Dim> shape; + typename Domain<T,Dim>::DomainSize domainSize; + Domain<T,Dim> ref, inBuffer, outBuffer; + + Domain<T,Dim>& in = inBuffer; + Domain<T,Dim>& out = outBuffer; + + std::array<int,Dim> order; + + shape.fill(8); + domainSize.fill(2*hysop::constants::pi); + + T eps = std::numeric_limits<T>::epsilon(); + const std::size_t N = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<std::size_t>()); + + ref.resize(domainSize).reshape(shape); + in = ref; + out = ref; + + Shape<Dim> maxOrder, testCases; + maxOrder.fill(p_maxOrder+1); + testCases.fill(nExtensionsPair); + Index<Dim> orderId(maxOrder); + Index<Dim> testCaseId; + std::size_t testCaseCount; + while(!(++orderId).atMaxId()) { + std::cout << " ::Order::" << orderId.ids() << (verbose ? "\n" : ""); + + std::array<T,3> meanDists; + meanDists.fill(0); + testCaseId.reset(testCases); + testCaseCount = testCaseId.maxId(); + while(!testCaseId.atMaxId()) { + std::copy(orderId.ids().begin(),orderId.ids().end(), order.begin()); + + /* generate transform configuration */ + std::array<std::pair<fft::Extension,fft::Extension>, Dim> extConfig; + for (std::size_t k=0; k<Dim; k++) { + std::size_t id = testCaseId[k]; + extConfig[k] = pext[id]; + if(pext[id].first==fft::Extension::NONE) + order[k] = 0; + } + fft::FftDomainConfiguration<Dim> domainConfig(extConfig, includePeriodicBds); + + const std::size_t orderSum = std::accumulate(order.begin(), order.end(), 0); + if(orderSum == 0) { + testCaseCount--; + ++testCaseId; + continue; + } + T orderPow = std::pow(T(10),T(orderSum)); + if(std::is_same<T,long double>::value) /* just in case long doubles are not hardware supported... */ + orderPow *= 1e3; + const auto criteria = std::make_tuple(orderPow*eps*N,orderPow*eps*sqrt(N),2*orderPow*eps); + + const auto f = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = func<T>(testCaseId[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= func<T>(testCaseId[d])(x[d]); + return val; + }; + const auto d = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = derivative<T>(testCaseId[0],order[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= derivative<T>(testCaseId[d],order[d])(x[d]); + return val; + }; + { + ref.resetDomainConfiguration(domainConfig.boundariesConfiguration()); + in = ref; + out = ref; + + in.apply(f); + ref.apply(d); + out.data().apply([](T& v){ v=T(0);}); + } + + solver::FftDiffSolver<T,Dim> solver(domainSize, domainConfig, FFTW_MEASURE, includePeriodicBds, includePeriodicBds); + solver.apply(in.data(), out.data(), order); + + std::stringstream ss; + ss << "["; + for (std::size_t k=0; k<Dim-1; k++) + ss << extConfig[k].first << "/" << extConfig[k].second << ","; + ss << extConfig[Dim-1].first << "/" << extConfig[Dim-1].second; + ss << "]"; + + const auto dist = out.distance(ref); + const bool pass = (std::get<0>(dist) < std::get<0>(criteria)) + && (std::get<1>(dist) < std::get<1>(criteria)) + && (std::get<2>(dist) < std::get<2>(criteria)); + + if((pass && verbose) || !pass) { + std::cout << (pass ? GREEN : RED); + std::cout << "\t" << std::setw(Dim*15) << ss.str() << " => " << (pass ? "OK" : "KO") + << " " << RESET << std::scientific << std::setprecision(2) << dist << std::endl; + } + if(!pass) { + //in.print("IN"); + //ref.print("REF"); + //out.print("OUT"); + std::cout << criteria << std::endl; + exit(EXIT_FAILURE); + } + + meanDists[0] += std::get<0>(dist); + meanDists[1] += std::get<1>(dist); + meanDists[2] += std::get<2>(dist); + + ++testCaseId; + } + for (std::size_t k = 0; k < 3; k++) + meanDists[k] /= T(testCaseCount); + std::cout << "=> mean distances over " << std::scientific << std::setprecision(1) << std::setw(4) + << testCaseCount << " testcases: " << meanDists; + for (std::size_t k = 0; k < 3; k++) + meanDists[k] = std::round(meanDists[k]/eps); + std::cout << " ~= " << std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; + } +} + +int main(int argc, const char *argv[]) { + + std::cout << "== TEST 1D - float ==" << std::endl; + test<float,1,false>(5); + std::cout << "== TEST 2D - float ==" << std::endl; + test<float,2,false>(3); + std::cout << "== TEST 3D - float ==" << std::endl; + test<float,3,false>(1); + std::cout << std::endl; + + std::cout << "== TEST 1D - double ==" << std::endl; + test<double,1,false>(5); + std::cout << "== TEST 2D - double ==" << std::endl; + test<double,2,false>(3); + std::cout << "== TEST 3D - double ==" << std::endl; + test<double,3,false>(1); + std::cout << std::endl; + + std::cout << "== TEST 1D - long double ==" << std::endl; + test<long double,1,false>(5); + std::cout << "== TEST 2D - long double ==" << std::endl; + test<long double,2,false>(3); + std::cout << "== TEST 3D - long double ==" << std::endl; + test<long double,3,false>(1); + std::cout << std::endl; + + return EXIT_SUCCESS; +} diff --git a/HySoP/src/hysop++/main/poissonSolver.cpp b/HySoP/src/hysop++/main/poissonSolver.cpp new file mode 100644 index 000000000..59cd79961 --- /dev/null +++ b/HySoP/src/hysop++/main/poissonSolver.cpp @@ -0,0 +1,182 @@ + +#include <cstdlib> + +#include "domain/domain.h" +#include "solver/fftPoissonSolver.h" +#include "data/multi_array/multi_array.h" +#include "utils/constants.h" +#include "domain/boundary.h" + +using namespace hysop; +using namespace hysop::domain; + +static constexpr std::size_t nBoundaries = 4; +static constexpr std::size_t nBoundaryPairs = 7; +static constexpr domain::Boundary bds[nBoundaries] = +{ domain::Boundary::NONE, domain::Boundary::HOMOGENEOUS_NEUMANN, domain::Boundary::HOMOGENEOUS_DIRICHLET, domain::Boundary::PERIODIC }; +static constexpr std::pair<domain::Boundary,domain::Boundary> pbds[nBoundaryPairs] { + std::make_pair(bds[3],bds[3]), //periodic-periodic + std::make_pair(bds[3],bds[3]), //periodic-periodic + std::make_pair(bds[2],bds[1]), //even-odd + std::make_pair(bds[1],bds[2]), //odd-even + std::make_pair(bds[2],bds[2]), //even-even + std::make_pair(bds[1],bds[1]), //odd-odd + std::make_pair(bds[0],bds[0]), //none-none +}; + +static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L }; + +template <typename T> +std::function<T(T)> func(std::size_t k) { + switch(k) { + case 0: return [=](T x) {return cos(T(freqs[0])*x);}; + case 1: return [=](T x) {return sin(T(freqs[1])*x);}; + case 2: return [=](T x) {return cos(T(freqs[2])*x);}; + case 3: return [=](T x) {return sin(T(freqs[3])*x);}; + case 4: return [=](T x) {return cos(T(freqs[4])*x);}; + case 5: return [=](T x) {return sin(T(freqs[5])*x);}; + default: return[=](T x) { return T(1); }; + } +} + +template <typename T, std::size_t Dim, bool verbose=false> +void test(bool includePeriodicBds=false) { + Shape<Dim> shape; + typename Domain<T,Dim>::DomainSize domainSize; + Domain<T,Dim> ref, inBuffer, outBuffer; + + Domain<T,Dim>& in = inBuffer; + Domain<T,Dim>& out = outBuffer; + + shape.fill(16); + domainSize.fill(2*hysop::constants::pi); + + T eps = std::numeric_limits<T>::epsilon(); + const std::size_t N = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<std::size_t>()); + + ref.resize(domainSize).reshape(shape); + in = ref; + out = ref; + + Shape<Dim> testCases; + testCases.fill(nBoundaryPairs); + Index<Dim> testCaseId(testCases); + std::array<T,3> meanDists{0}; + std::size_t testCaseCount = testCaseId.maxId()-1; + + if(verbose) + std::cout << std::endl; + + while(testCaseId() != testCaseId.maxId()-1) { + + /* generate transform configuration */ + std::size_t orderSum = 0; + std::array<std::pair<domain::Boundary,domain::Boundary>, Dim> bdsConfig; + T W2sum = T(0); + for (std::size_t k=0; k<Dim; k++) { + std::size_t id = testCaseId[k]; + bdsConfig[k] = pbds[id]; + if(bdsConfig[k].first != domain::Boundary::NONE) { + W2sum += freqs[id]*freqs[id]; + orderSum+=2; + } + } + domain::DomainConfiguration<Dim> domainConfig(bdsConfig, includePeriodicBds); + + T orderPow = std::pow(T(10),T(orderSum)); + if(std::is_same<T,long double>::value) /* just in case long doubles are not hardware supported... */ + orderPow *= 1e3; + const auto criteria = std::make_tuple(orderPow*eps*N,orderPow*eps*sqrt(N),2*orderPow*eps); + + const auto phi = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = func<T>(testCaseId[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= func<T>(testCaseId[d])(x[d]); + return val; + }; + const auto f = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + return -W2sum*phi(x); + }; + + { + ref.resetDomainConfiguration(domainConfig); + in = ref; + out = ref; + + in.apply(f); + ref.apply(phi); + out.data().apply([](T& v){ v=T(0);}); + } + + solver::FftPoissonSolver<T,Dim> solver(domainSize, domainConfig, FFTW_MEASURE, includePeriodicBds, includePeriodicBds); + solver.apply(in.data(), out.data()); + + std::stringstream ss; + ss << "["; + for (std::size_t k=0; k<Dim-1; k++) + ss << bdsConfig[k].first << "/" << bdsConfig[k].second << ","; + ss << bdsConfig[Dim-1].first << "/" << bdsConfig[Dim-1].second; + ss << "]"; + + const auto dist = out.distance(ref); + const bool pass = (std::get<0>(dist) < std::get<0>(criteria)) + && (std::get<1>(dist) < std::get<1>(criteria)) + && (std::get<2>(dist) < std::get<2>(criteria)); + + if((pass && verbose) || !pass) { + std::cout << (pass ? GREEN : RED); + std::cout << "\t" << std::setw(Dim*15) << ss.str() << " => " << (pass ? "OK" : "KO") + << " " << RESET << std::scientific << std::setprecision(2) << dist << std::endl; + } + if(!pass) { + //in.print("IN"); + //ref.print("REF"); + //out.print("OUT"); + std::cout << "\t\tTest Failed... Criteria was " << criteria << "." << std::endl; + exit(EXIT_FAILURE); + } + + meanDists[0] += std::get<0>(dist); + meanDists[1] += std::get<1>(dist); + meanDists[2] += std::get<2>(dist); + + ++testCaseId; + } + + for (std::size_t k = 0; k < 3; k++) + meanDists[k] /= T(testCaseCount); + std::cout << "\t=> mean distances over " << std::scientific << std::setprecision(1) << std::setw(4) + << testCaseCount << " testcases: " << meanDists; + for (std::size_t k = 0; k < 3; k++) + meanDists[k] = std::round(meanDists[k]/eps); + std::cout << " ~= " << std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; +} + +int main(int argc, const char *argv[]) { + + std::cout << "== TEST 1D - float =="; + test<float,1,true>(true); + std::cout << "== TEST 2D - float =="; + test<float,2,false>(); + std::cout << "== TEST 3D - float =="; + test<float,3,false>(); + std::cout << std::endl; + + std::cout << "== TEST 1D - double =="; + test<double,1,false>(); + std::cout << "== TEST 2D - double =="; + test<double,2,false>(); + std::cout << "== TEST 3D - double =="; + test<double,3,false>(); + std::cout << std::endl; + + std::cout << "== TEST 1D - long double =="; + test<long double,1,false>(); + std::cout << "== TEST 2D - long double =="; + test<long double,2,false>(); + std::cout << "== TEST 3D - long double =="; + test<long double,3,false>(); + std::cout << std::endl; + + return EXIT_SUCCESS; +} diff --git a/HySoP/src/hysop++/src/data/accumulatorIndex.h b/HySoP/src/hysop++/src/data/accumulatorIndex.h new file mode 100644 index 000000000..b794002c2 --- /dev/null +++ b/HySoP/src/hysop++/src/data/accumulatorIndex.h @@ -0,0 +1,65 @@ + +#ifndef HYSOP_ACCUMULATORINDEX_H +#define HYSOP_ACCUMULATORINDEX_H + +#include <array> +#include <cassert> +#include <vector> +#include <functional> + +#include "data/index.h" + +namespace hysop { + + template <typename T, std::size_t Dim, typename Source = std::array<std::vector<T>,Dim>, typename Functor = std::function<T(const T&, const T&)>> + struct AccumulatorIndex : public Index<Dim> { + + private: + using super = Index<Dim>; + public: + using Indices = typename super::Indices; + using Dimension = typename super::Dimension; + + public: + AccumulatorIndex(const AccumulatorIndex& idx) = default; + ~AccumulatorIndex() = default; + + template <typename DimArray=Dimension, typename AccumulatorIndexArray=Indices> + AccumulatorIndex(const DimArray& p_dim = DimArray{0}, + const Indices &p_ids = AccumulatorIndexArray{0}): + super(p_dim, p_ids), m_accumulatedData{0}, m_sourceData(nullptr), m_functor(std::plus<T>()) { + } + + AccumulatorIndex& setAccumulatorSource(const Source& p_source) { m_sourceData = &p_source; init(); return *this; } + AccumulatorIndex& setAccumulatorFunction(const Functor& p_functor) { m_functor = p_functor; init(); return *this; } + + const T& accumulatedVal() const { return m_accumulatedData[Dim-1]; } + + protected: + void init() { + if(m_sourceData != nullptr) { + m_accumulatedData[0] = (*m_sourceData)[0][this->operator[](0)]; + for (std::size_t d=0; d<Dim-1; d++) + m_accumulatedData[d+1] = m_functor(m_accumulatedData[d],(*m_sourceData)[d+1][this->operator[](d+1)]); + } + } + + virtual void onIndexChange(std::size_t pos, std::ptrdiff_t offset) final override { + if(m_sourceData != nullptr) { + assert(pos < Dim); + m_accumulatedData[pos] = (pos==0 ? (*m_sourceData)[0][this->operator[](0)] : m_functor(m_accumulatedData[pos-1],(*m_sourceData)[pos][this->operator[](pos)])); + for (std::size_t d=pos; d<Dim-1; d++) + m_accumulatedData[d+1] = m_functor(m_accumulatedData[d],(*m_sourceData)[d+1][this->operator[](d+1)]); + } + }; + + protected: + std::array<T,Dim> m_accumulatedData; + const Source* m_sourceData; + Functor m_functor; + }; + +} + +#endif /* end of include guard: HYSOP_ACCUMULATORINDEX_H */ + diff --git a/HySoP/src/hysop++/src/data/index.h b/HySoP/src/hysop++/src/data/index.h new file mode 100644 index 000000000..3166beb24 --- /dev/null +++ b/HySoP/src/hysop++/src/data/index.h @@ -0,0 +1,125 @@ + +#ifndef HYSOP_INDEX_H +#define HYSOP_INDEX_H + +#include <array> +#include <cassert> + +namespace hysop { + + template <std::size_t Dim> + struct Index { + typedef boost::array<std::ptrdiff_t,Dim> Indices; + typedef boost::array<std::size_t,Dim> Dimension; + + Index(const Index& idx) = default; + ~Index() = default; + + template <typename DimArray=Dimension, typename IndexArray=Indices> + Index(const DimArray& p_dim = DimArray{0}, + const Indices &p_ids = IndexArray{0}) : + m_dim(), m_ids(), m_id(0) { + for (std::size_t d = 0; d < Dim; d++) { + m_dim[d] = p_dim[d]; + m_ids[d] = p_ids[d]; + } + initializeId(); + } + + template <typename DimArray, typename IndexArray=Indices> + Index& reset(const DimArray& p_dim, const IndexArray &p_ids = IndexArray{0}) { + for (std::size_t d = 0; d < Dim; d++) { + m_dim[d] = p_dim[d]; + m_ids[d] = p_ids[d]; + } + initializeId(); + return *this; + } + + Index& setIndexToMinusOne() { + for (std::size_t d=0; d<Dim-1; d++) + m_ids[d] = 0; + m_ids[Dim-1] = -1; + initializeId(); + return *this; + } + + std::size_t id() const { return m_id; }; + std::size_t maxId() const { return m_maxId; } + bool atMaxId() const { return m_id == m_maxId; } + + const Indices& ids() const { return m_ids; }; + const Dimension& dim() const { return m_dim; }; + + std::ptrdiff_t operator[](std::size_t k) const { return m_ids[k]; } + std::size_t operator()() const { return m_id; } + + virtual void onIndexChange(std::size_t pos, std::ptrdiff_t offset) {}; + virtual void onIndexOverflow(std::size_t pos) {}; + virtual void onIndexUndeflow(std::size_t pos) {}; + +__attribute__((optimize("unroll-loops"))) + Index& operator++() { //prefix + for (int d = Dim-1; d >=0; d--) { + m_ids[d]++; + if(m_ids[d]==std::ptrdiff_t(m_dim[d])) { + m_ids[d]=0; + this->onIndexOverflow(d); + } + else { + this->onIndexChange(d, +1); + break; + } + } + m_id++; + return *this; + }; +__attribute__((optimize("unroll-loops"))) + Index& operator--() { //prefix + for (int d = 0; d < Dim; d++) { + if(m_ids[d]!=0) { + this->onIndexChange(d, -1); + m_ids[d]--; + break; + } + else { + this->onIndexUndeflow(d); + } + } + m_id--; + return *this; + }; + + Index operator++(int) { //postfix + Index result(*this); + ++(*this); + return result; + }; + Index operator--(int) { //postfix + Index result(*this); + --(*this); + return result; + }; + + protected: + void initializeId() { + m_id = 0; + m_maxId = 1; + for (std::size_t d = 0; d < Dim-1; d++) { + m_id = (m_id + m_ids[d]) * m_dim[d+1]; + m_maxId*=m_dim[d]; + } + m_id += m_ids[Dim-1]; + m_maxId*=m_dim[Dim-1]; + } + + private: + Dimension m_dim; + Indices m_ids; + std::size_t m_id, m_maxId; + }; + +} + +#endif /* end of include guard: HYSOP_INDEX_H */ + diff --git a/HySoP/src/hysop++/src/data/memory/fftwAllocator.h b/HySoP/src/hysop++/src/data/memory/fftwAllocator.h new file mode 100644 index 000000000..d2b8423f3 --- /dev/null +++ b/HySoP/src/hysop++/src/data/memory/fftwAllocator.h @@ -0,0 +1,53 @@ + +#ifndef HYSOP_FFTW_ALLOCATOR_H +#define HYSOP_FFTW_ALLOCATOR_H + +#include <limits> +#include <fftw3.h> + +namespace hysop { + namespace data { + namespace memory { + + /* FftwAllocator designed to correctly align data for fftw */ + template <typename T> + struct FftwAllocator { + using value_type = T; + using const_pointer = const T*; + + FftwAllocator() = default; + + template <class U> + FftwAllocator(const FftwAllocator<U>&) {} + + T* allocate(std::size_t n, const_pointer hint=nullptr) { + if (n <= std::numeric_limits<std::size_t>::max() / sizeof(T)) { + if (auto ptr = fftw_malloc(n * sizeof(T))) { + return static_cast<T*>(ptr); + } + } + throw std::bad_alloc(); + } + void deallocate(T* ptr, std::size_t n) { + fftw_free(ptr); + } + void destroy(T* ptr) { + ptr->~T(); + } + }; + + template <typename T, typename U> + inline bool operator == (const FftwAllocator<T>&, const FftwAllocator<U>&) { + return true; + } + + template <typename T, typename U> + inline bool operator != (const FftwAllocator<T>& a, const FftwAllocator<U>& b) { + return !(a == b); + } + + } /* end of namespace memory */ + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_FFTW_ALLOCATOR_H */ diff --git a/HySoP/src/hysop++/src/data/memory/minimalAllocator.h b/HySoP/src/hysop++/src/data/memory/minimalAllocator.h new file mode 100644 index 000000000..2460db0db --- /dev/null +++ b/HySoP/src/hysop++/src/data/memory/minimalAllocator.h @@ -0,0 +1,52 @@ + +#ifndef HYSOP_MINIMAL_ALLOCATOR_H +#define HYSOP_MINIMAL_ALLOCATOR_H + +#include <limits> + +namespace hysop { + namespace data { + namespace memory { + + /* Minimal MinimalAllocator for boost and std libs */ + template <typename T> + struct MinimalAllocator { + using value_type = T; + using const_pointer = const T*; + + MinimalAllocator() = default; + + template <class U> + MinimalAllocator(const MinimalAllocator<U>&) {} + + T* allocate(std::size_t n, const_pointer hint=nullptr) { + if (n <= std::numeric_limits<std::size_t>::max() / sizeof(T)) { + if (auto ptr = std::malloc(n * sizeof(T))) { + return static_cast<T*>(ptr); + } + } + throw std::bad_alloc(); + } + void deallocate(T* ptr, std::size_t n) { + std::free(ptr); + } + void destroy(T* ptr) { + ptr->~T(); + } + }; + + template <typename T, typename U> + inline bool operator == (const MinimalAllocator<T>&, const MinimalAllocator<U>&) { + return true; + } + + template <typename T, typename U> + inline bool operator != (const MinimalAllocator<T>& a, const MinimalAllocator<U>& b) { + return !(a == b); + } + + } /* end of namespace memory */ + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_MINIMAL_ALLOCATOR_H */ diff --git a/HySoP/src/hysop++/src/data/multi_array/< b/HySoP/src/hysop++/src/data/multi_array/< new file mode 100644 index 000000000..9b142c71f --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/< @@ -0,0 +1,352 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#define HYSOP_MULTI_ARRAY_H + +/**********************************/ +/*** Hysop boost multi array wrapper ***/ +/********************************/ + +#include "const_multi_array_view.h" +#include "multi_array_view.h" + +namespace hysop { + namespace data { + + /* forward declaration */ + FORWARD_DECLARE_TYPES() + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim, typename Allocator> + class multi_array : public boost_multi_array<T,Dim,Allocator> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_multi_array<T,Dim,Allocator>; + using super_ref = boost_multi_array_ref<T,Dim>; + public: + PUBLIC_CLASS_TYPES() + + public: + multi_array(const extents_gen<Dim>& extents = extents_gen<Dim>()); + multi_array(const Shape<Dim>& shape); + + multi_array(const multi_array& other); + multi_array(multi_array&& other); + + explicit multi_array(const array_ref& other); + explicit multi_array(const array_view& other); + explicit multi_array(const const_array_ref& other); + explicit multi_array(const const_array_view& other); + + explicit multi_array(const boost_multi_array<T,Dim,Allocator>& other); + explicit multi_array(const boost_multi_array_ref<T,Dim>& other); + explicit multi_array(const boost_multi_array_view<T,Dim>& other); + explicit multi_array(const boost_const_multi_array_ref<T,Dim>& other); + explicit multi_array(const boost_const_multi_array_view<T,Dim>& other); + explicit multi_array(boost_multi_array<T,Dim,Allocator>&& other); + + multi_array& operator=(const multi_array& other); + multi_array& operator=(const array_ref& ref); + multi_array& operator=(const array_view& view); + multi_array& operator=(const const_array_ref& ref); + multi_array& operator=(const const_array_view& view); + multi_array& operator=(multi_array&& other); + + operator array_ref(); + operator const_array_ref() const; + + public: + + multi_array& reshape(const Shape<Dim>& shape); + + /* print data */ + const multi_array& print(const std::string& name, std::ostream& os = std::cout, + unsigned int precision=2u, unsigned int width=6u) const; + + /* Apply function to all elements of the view */ + /* (1) apply func(T&) on all elements */ + /* (2) apply func(T&, const Index<Dim>&) on all elements */ + /* (3) apply func(T&, const Index<Dim>&, farg0, fargs...) on all elements */ + multi_array& apply(const std::function<void(T&)>& func); + multi_array& apply(const std::function<void(T&, const Index<Dim>&)>& func); + template <typename Functor, typename Arg0, typename... Args> + multi_array& apply(const Functor& func, Arg0&& farg0, Args&&... fargs); + + /* boolean array_view specific functions */ + ENABLE_IF_BOOL(bool) all() const; + ENABLE_IF_BOOL(bool) any() const; + ENABLE_IF_BOOL(bool) none() const; + + /* real and complex data accessors, usefull for FFT like transforms */ + ENABLE_IF_REAL ( T*) rdata() { return this->origin(); } + ENABLE_IF_REAL (const T*) rdata() const { return this->origin(); } + + ENABLE_IF_COMPLEX( T*) cdata() { return this->origin(); } + ENABLE_IF_COMPLEX(const T*) cdata() const { return this->origin(); } + + ENABLE_IF_COMPLEX( typename fft::fftw_complex_type<T>::std_type*) std_cdata() { + return reinterpret_cast< typename fft::fftw_complex_type<T>::std_type*>(this->origin()); + } + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::std_type*) std_cdata() const { + return reinterpret_cast<const typename fft::fftw_complex_type<T>::std_type*>(this->origin()); + } + ENABLE_IF_COMPLEX( typename fft::fftw_complex_type<T>::fftw_type*) fftw_cdata() { + return reinterpret_cast< typename fft::fftw_complex_type<T>::fftw_type*>(this->origin()); + } + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::fftw_type*) fftw_cdata() const { + return reinterpret_cast<const typename fft::fftw_complex_type<T>::fftw_type*>(this->origin()); + } + + ENABLE_IF_COMPLEX( typename fft::fftw_complex_type<T>::value_type*) asRealData() { + return reinterpret_cast< typename fft::fftw_complex_type<T>::value_type*>(this->origin()); + } + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::value_type*) asRealData() const { + return reinterpret_cast< const typename fft::fftw_complex_type<T>::value_type*>(this->origin()); + } + ENABLE_IF_REAL ( typename fft::fftw_complex_type<T>::std_type*) asStdComplexData() { + return reinterpret_cast< typename fft::fftw_complex_type<T>::std_type*>(this->origin()); + } + ENABLE_IF_REAL (const typename fft::fftw_complex_type<T>::std_type*) asStdComplexdata() const { + return reinterpret_cast<const typename fft::fftw_complex_type<T>::std_type*>(this->origin()); + } + ENABLE_IF_REAL ( typename fft::fftw_complex_type<T>::fftw_type*) asFftwComplexdata() { + return reinterpret_cast< typename fft::fftw_complex_type<T>::fftw_type*>(this->origin()); + } + ENABLE_IF_REAL (const typename fft::fftw_complex_type<T>::fftw_type*) asFftwComplexData() const { + return reinterpret_cast<const typename fft::fftw_complex_type<T>::fftw_type*>(this->origin()); + } + + protected: + static extents_gen<Dim> shapeToExtents(const Shape<Dim> &shape); + }; + + + /* Implementation */ + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const extents_gen<Dim>& extents): + boost_multi_array<T,Dim,Allocator>(extents) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const Shape<Dim>& shape): + boost_multi_array<T,Dim,Allocator>(shapeToExtents(shape)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const multi_array& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const array_view& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array_view<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const const_array_view& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array_view<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(multi_array&& other): + boost_multi_array<T,Dim,Allocator>(static_cast<boost_multi_array<T,Dim,Allocator>&&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_multi_array_view<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_const_multi_array_view<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_multi_array_ref<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_const_multi_array_ref<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(boost_multi_array<T,Dim,Allocator>&& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const array_view& other) { + this->reshape(other.shape()); + /* cast obligatory to avoid shape() function aliasing */ + super::operator=(dynamic_cast<const boost_multi_array_view<T,Dim>&>(other)); + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const const_array_view& other) { + this->reshape(other.shape()); + /* cast obligatory to avoid shape() function aliasing */ + super::operator=(dynamic_cast<const boost_multi_array_view<T,Dim>&>(other)); + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const multi_array& other) { + this->reshape(other.shape()); + /* cast obligatory to avoid shape() function aliasing */ + super::operator=(dynamic_cast<const boost_multi_array_ref<T,Dim>&>(other)); + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(multi_array&& other) { + super::operator=(other); + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::operator array_view() { + return this->operator[](hysop::utils::buildView<Dim>()); + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::operator const_array_view() const { + return this->operator[](hysop::utils::buildView<Dim>()); + } + + + template <typename T, std::size_t Dim, typename Allocator> + Shape<Dim> multi_array<T,Dim,Allocator>::shape() const { + Shape<Dim> shape; + const std::size_t* extents = this->super::shape(); + for (std::size_t d = 0; d < Dim; d++) + shape[d] = static_cast<std::size_t>(extents[d]); + return shape; + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::reshape(const Shape<Dim>& shape) { + boost::array<int,Dim> extents; + for (std::size_t d = 0; d < Dim; d++) + extents[d] = static_cast<int>(shape[d]); + this->resize(extents); + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + typename multi_array<T,Dim,Allocator>::template extents_gen<Dim> multi_array<T,Dim,Allocator>::shapeToExtents(const Shape<Dim> &shape) { + return utils::buildExtents(shape); + } + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::apply(const std::function<void(T&)>& func) { + T* data = this->origin(); + for (std::size_t k = 0; k < this->num_elements(); k++) { + func(data[k]); + } + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::apply(const std::function<void(T&, const Index<Dim>&)>& func) { + hysop::Index<Dim> idx(this->shape()); + T* data = this->origin(); + for (std::size_t k = 0; k < this->num_elements(); k++) { + func(data[k], idx); + ++idx; + } + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + template <typename Functor, typename Arg0, typename... Args> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::apply(const Functor& func, Arg0&& farg0, Args&&... fargs) { + hysop::Index<Dim> idx(this->shape()); + T* data = this->origin(); + for (std::size_t k = 0; k < this->num_elements(); k++) { + func(data[k], idx, std::forward<Arg0>(farg0), std::forward<Args>(fargs)...); + ++idx; + } + return *this; + } + + template <typename T, std::size_t Dim, typename Allocator> + const multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::print(const std::string& name, std::ostream& os, + unsigned int precision, unsigned int width) const { + auto S = this->super::shape(); + std::size_t id = 0; + + os << name << " = ["; + if(Dim==1 || Dim>3) { + for(std::size_t k=0; k<this->num_elements(); k++) { + T x = this->data()[id++]; + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; + } + } + else if(Dim==2) { + std::cout << std::endl; + for(std::size_t i=0; i<S[0]; i++) { + os << "\t["; + for(std::size_t j=0; j<S[1]; j++) { + T x = this->data()[id++]; + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; + } + os << "]" << std::endl; + } + } + else { + std::cout << std::endl; + for(std::size_t i=0; i<S[0]; i++) { + os << "\t[["; + for(std::size_t j=0; j<S[1]; j++) { + if(j>0) + os << "\t ["; + for(std::size_t k=0; k<S[2]; k++) { + T x = this->data()[id++]; + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; + } + if(j!=S[1]-1) + os << "]" << std::endl; + else + os << "]]," << std::endl; + } + } + } + os << "];" << std::endl; + return *this; + } + + /* boolean specific functions */ + template <typename T, std::size_t Dim, typename Allocator> + template <typename TT> + typename std::enable_if<std::is_same<TT,bool>::value, bool>::type multi_array<T,Dim,Allocator>::all() const { + for(std::size_t k=0; k<this->num_elements(); k++) { + if(!this->data()[k]) + return false; + } + return true; + } + + template <typename T, std::size_t Dim, typename Allocator> + template <typename TT> + typename std::enable_if<std::is_same<TT,bool>::value, bool>::type multi_array<T,Dim,Allocator>::any() const { + for(std::size_t k=0; k<this->num_elements(); k++) { + if(this->data()[k]) + return true; + } + return false; + } + + template <typename T, std::size_t Dim, typename Allocator> + template <typename TT> + typename std::enable_if<std::is_same<TT,bool>::value, bool>::type multi_array<T,Dim,Allocator>::none() const { + for(std::size_t k=0; k<this->num_elements(); k++) { + if(this->data()[k]) + return false; + } + return true; + } + + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#undef ENABLE_IF_BOOL +#undef ENABLE_IF_REAL +#undef ENABLE_IF_COMPLEX + +#include "multi_array_ext.h" + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_H */ diff --git a/HySoP/src/hysop++/src/data/multi_array/const_multi_array_ref.h b/HySoP/src/hysop++/src/data/multi_array/const_multi_array_ref.h new file mode 100644 index 000000000..cc0a2f2c4 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/const_multi_array_ref.h @@ -0,0 +1,53 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_CONST_MULTI_ARRAY_REF_H +#define HYSOP_CONST_MULTI_ARRAY_REF_H + +namespace hysop { + namespace data { + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim> + class const_multi_array_ref : public boost_const_multi_array_ref<T,Dim> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_const_multi_array_ref<T,Dim>; + public: + PUBLIC_CLASS_TYPES() + + public: + const_multi_array_ref(const const_multi_array_ref<T,Dim>& ref) = default; + const_multi_array_ref& operator=(const const_multi_array_ref<T,Dim>& other) = default; + + const_multi_array_ref(const boost_const_multi_array_ref<T,Dim>& ref); + const_multi_array_ref& operator=(const boost_const_multi_array_ref<T,Dim>& other); + + public: + PUBLIC_CONST_REF_INTERFACE(SINGLE_ARG(const_multi_array_ref<T,Dim>)) + }; + + + /* Implementation */ + template <typename T, std::size_t Dim> + const_multi_array_ref<T,Dim>::const_multi_array_ref(const boost_const_multi_array_ref<T,Dim>& ref) : + super(ref) { + } + + template <typename T, std::size_t Dim> + const_multi_array_ref<T,Dim>& const_multi_array_ref<T,Dim>::operator=(const boost_const_multi_array_ref<T,Dim>& other) { + super::operator=(other); + return *this; + } + + CONST_REF_IMPLEMENTATION(SINGLE_ARG(const_multi_array_ref<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_CONST_MULTI_ARRAY_REF_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/const_multi_array_view.h b/HySoP/src/hysop++/src/data/multi_array/const_multi_array_view.h new file mode 100644 index 000000000..335433d8a --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/const_multi_array_view.h @@ -0,0 +1,53 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_CONST_MULTI_ARRAY_VIEW_H +#define HYSOP_CONST_MULTI_ARRAY_VIEW_H + +namespace hysop { + namespace data { + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim> + class const_multi_array_view : public boost_const_multi_array_view<T,Dim> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_const_multi_array_view<T,Dim>; + public: + PUBLIC_CLASS_TYPES() + + public: + const_multi_array_view(const const_multi_array_view<T,Dim>& view) = default; + const_multi_array_view& operator=(const const_multi_array_view<T,Dim>& other) = default; + + const_multi_array_view(const boost_const_multi_array_view<T,Dim>& view); + const_multi_array_view& operator=(const boost_const_multi_array_view<T,Dim>& other); + + public: + PUBLIC_CONST_VIEW_INTERFACE(SINGLE_ARG(const_multi_array_view<T,Dim>)) + }; + + + /* Implementation */ + template <typename T, std::size_t Dim> + const_multi_array_view<T,Dim>::const_multi_array_view(const boost_const_multi_array_view<T,Dim>& view) : + super(view) { + } + + template <typename T, std::size_t Dim> + const_multi_array_view<T,Dim>& const_multi_array_view<T,Dim>::operator=(const boost_const_multi_array_view<T,Dim>& other) { + super::operator=(other); + return *this; + } + + CONST_VIEW_IMPLEMENTATION(SINGLE_ARG(const_multi_array_view<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_CONST_MULTI_ARRAY_VIEW_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array.h b/HySoP/src/hysop++/src/data/multi_array/multi_array.h new file mode 100644 index 000000000..6224ac4ca --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array.h @@ -0,0 +1,51 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#define HYSOP_MULTI_ARRAY_H + +#include <boost/multi_array.hpp> +#include "utils/default.h" + +/****************************************/ +/*** Hysop boost multi array wrapper ****/ +/****************************************/ + +namespace hysop { + namespace data { + + /* forward declaration of types */ + template <typename T, std::size_t Dim, typename Allocator = hysop::_default::allocator<T>> + class multi_array; + template <typename T, std::size_t Dim> + class multi_array_ref; + template <typename T, std::size_t Dim> + class const_multi_array_ref; + template <typename T, std::size_t Dim> + class multi_array_view; + template <typename T, std::size_t Dim> + class const_multi_array_view; + + template <typename T, std::size_t Dim, typename Allocator = hysop::_default::allocator<T>> + using boost_multi_array = boost::multi_array<T,Dim,Allocator>; + template <typename T, std::size_t Dim> + using boost_multi_array_ref = boost::multi_array_ref<T,Dim>; + template <typename T, std::size_t Dim> + using boost_const_multi_array_ref = boost::const_multi_array_ref<T,Dim>; + template <typename T, std::size_t Dim> + using boost_multi_array_view = boost::detail::multi_array::multi_array_view<T,Dim>; + template <typename T, std::size_t Dim> + using boost_const_multi_array_view = boost::detail::multi_array::const_multi_array_view<T,Dim>; + } +} + +#include "data/multi_array/multi_array_defines.h" + +#include "data/multi_array/const_multi_array_view.h" +#include "data/multi_array/multi_array_view.h" +#include "data/multi_array/const_multi_array_ref.h" +#include "data/multi_array/multi_array_ref.h" +#include "data/multi_array/multi_array_impl.h" +#include "data/multi_array/multi_array_ext.h" + +#include "data/multi_array/multi_array_clean.h" + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_H */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_clean.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_clean.h new file mode 100644 index 000000000..0f8c988cc --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_clean.h @@ -0,0 +1,47 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_CLEAN_H +#define HYSOP_MULTI_ARRAY_CLEAN_H + +/* clean all macros used to generate hysop multiarray */ + +#undef SINGLE_ARG +#undef PUBLIC_CLASS_TYPES + +#undef ENABLE_IF_BOOL +#undef ENABLE_IF_REAL +#undef ENABLE_IF_COMPLEX + +#undef PUBLIC_COMMON_CONST_INTEFACE +#undef PUBLIC_CONST_REF_INTERFACE +#undef PUBLIC_CONST_VIEW_INTERFACE +#undef PUBLIC_COMMON_NON_CONST_INTERFACE +#undef PUBLIC_NON_CONST_REF_INTERFACE +#undef PUBLIC_NON_CONST_VIEW_INTERFACE + +#undef LOOP_VARNAME +#undef LOOP_OVER_ALL_REF_ELEMENTS +#undef REF_DATA_ACCESS +#undef LOOP_OVER_ALL_VIEW_ELEMENTS +#undef VIEW_DATA_ACCESS + +#undef COMMON_CONST_IMPLEMENTATION +#undef LOOP_DEPENDENT_CONST_IMPLEMENTATION +#undef CONST_REF_IMPL +#undef CONST_VIEW_IMPL +#undef CONST_REF_IMPLEMENTATION +#undef CONST_VIEW_IMPLEMENTATION + +#undef COMMON_NON_CONST_IMPLEMENTATION +#undef LOOP_DEPENDENT_NON_CONST_IMPLEMENTATION +#undef NON_CONST_REF_IMPL +#undef NON_CONST_VIEW_IMPL +#undef NON_CONST_REF_IMPLEMENTATION +#undef NON_CONST_VIEW_IMPLEMENTATION + +#endif /* end of include guard: MULTI_ARRAY_CLEAN_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h new file mode 100644 index 000000000..c28bd8ec0 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h @@ -0,0 +1,374 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_DEFINES_H +#define HYSOP_MULTI_ARRAY_DEFINES_H + +#include <algorithm> +#include <iostream> +#include <iomanip> +#include <cmath> + +#include "utils/utils.h" +#include "utils/types.h" +#include "data/index.h" +#include "fft/fftwComplex.h" + + +/* helper macros to generate all hysop multi_array wrappers */ +#define SINGLE_ARG(...) __VA_ARGS__ + +/* declaration of all public class types */ +#define PUBLIC_CLASS_TYPES() \ + template <typename Alloc> \ + using array = multi_array<T,Dim,Alloc>; \ + \ + using array_ref = multi_array_ref<T,Dim>; \ + using array_view = multi_array_view<T,Dim>; \ + using const_array_ref = const_multi_array_ref<T,Dim>; \ + using const_array_view = const_multi_array_view<T,Dim>; \ + \ + template <std::size_t NExtents> \ + using extents_gen = boost::detail::multi_array::extent_gen<NExtents>; \ + \ + template <std::size_t NumRanges, std::size_t NumDim> \ + using index_gen = boost::detail::multi_array::index_gen<NumRanges, NumDim>; \ + \ + using index_range = boost::multi_array_types::index_range; + + +/* various enable if macros */ +#define ENABLE_IF_BOOL(TYPE,DEFAULT) \ + template <typename TT DEFAULT> \ + typename std::enable_if< \ + std::is_same<TT,bool>::value \ + , TYPE>::type +#define ENABLE_IF_REAL(TYPE,DEFAULT) \ + template <typename TT DEFAULT> \ + typename std::enable_if< \ + std::is_same<TT,float> ::value \ + || std::is_same<TT,double> ::value \ + || std::is_same<TT,long double>::value \ + , TYPE>::type +#define ENABLE_IF_COMPLEX(TYPE,DEFAULT) \ + template <typename TT DEFAULT> \ + typename std::enable_if< \ + std::is_same<TT,std::complex<float>> ::value \ + || std::is_same<TT,std::complex<double>> ::value \ + || std::is_same<TT,std::complex<long double>>::value \ + || std::is_same<TT,fftwf_complex> ::value \ + || std::is_same<TT, fftw_complex> ::value \ + || std::is_same<TT,fftwl_complex> ::value \ + || std::is_same<TT,fftwq_complex> ::value \ + , TYPE>::type + +/* CLASS INTERFACES */ +/* public const common interfaces for views and referencess */ +#define PUBLIC_COMMON_CONST_INTERFACE(TYPENAME) \ + /* shape related */ \ + Shape<Dim> shape() const; \ + \ + /* print data */ \ + const TYPENAME& print(const std::string& name, std::ostream& os = std::cout, \ + unsigned int precision=2u, unsigned int width=6u) const; \ + \ + /* boolean array_view specific functions */ \ + ENABLE_IF_BOOL(bool,=T) all() const; \ + ENABLE_IF_BOOL(bool,=T) any() const; \ + ENABLE_IF_BOOL(bool,=T) none() const; + + +/* public const reference interface */ +#define PUBLIC_CONST_REF_INTERFACE(TYPENAME) \ + /* common const interface */ \ + PUBLIC_COMMON_CONST_INTERFACE(SINGLE_ARG(TYPENAME)) \ + /* real and complex data accessors, usefull for FFT like transforms */ \ + ENABLE_IF_REAL (const T*,=T) rdata() const; \ + ENABLE_IF_REAL (const typename fft::fftw_complex_type<T>::std_type*,=T) asStdComplexData() const; \ + ENABLE_IF_REAL (const typename fft::fftw_complex_type<T>::fftw_type*,=T) asFftwComplexData() const; \ + ENABLE_IF_COMPLEX(const T*,=T) cdata() const; \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::std_type*,=T) std_cdata() const; \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::fftw_type*,=T) fftw_cdata() const; \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::value_type*,=T) asRealData() const; + + +/* public const view interface */ +#define PUBLIC_CONST_VIEW_INTERFACE(TYPENAME) \ + /* common const interface */ \ + PUBLIC_COMMON_CONST_INTERFACE(SINGLE_ARG(TYPENAME)) + + +/* non const interfaces */ +#define PUBLIC_COMMON_NON_CONST_INTERFACE(TYPENAME) \ + /* Apply function to all elements */ \ + TYPENAME& apply(const std::function<void(T&)>& func); \ + TYPENAME& apply(const std::function<void(T&, const Index<Dim>&)>& func); \ + template <typename Functor, typename Arg0, typename... Args> \ + TYPENAME& apply(const Functor& func, Arg0&& farg0, Args&&... fargs); + + +#define PUBLIC_NON_CONST_REF_INTERFACE(TYPENAME) \ + /* common non const interface */ \ + PUBLIC_COMMON_NON_CONST_INTERFACE(SINGLE_ARG(TYPENAME)) \ + /* real and complex data accessors, usefull for FFT like transforms */ \ + ENABLE_IF_REAL(T*,=T) rdata(); \ + ENABLE_IF_REAL(typename fft::fftw_complex_type<T>::std_type*,=T) asStdComplexData(); \ + ENABLE_IF_REAL(typename fft::fftw_complex_type<T>::fftw_type*,=T) asFftwComplexData(); \ + ENABLE_IF_COMPLEX(T*,=T) cdata(); \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::std_type*,=T) std_cdata(); \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::fftw_type*,=T) fftw_cdata(); \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::value_type*,=T) asRealData(); + + +#define PUBLIC_NON_CONST_VIEW_INTERFACE(TYPENAME) \ + /* common non const interface */ \ + PUBLIC_COMMON_NON_CONST_INTERFACE(SINGLE_ARG(TYPENAME)) + + + +/* CLASS IMPLEMENTATIONS */ + +/* Loop dependant implementation macros (references contain contiguous data but not the views) */ +#define LOOP_VARNAME multi_array_index +#define NO_DEFAULT_TEMPLATES + +/* A reference has a contiguous memory layout and can be accessed by data offsets */ +#define LOOP_OVER_ALL_REF_ELEMENTS(ARRAY) \ + for (std::size_t LOOP_VARNAME=0; LOOP_VARNAME<((ARRAY).num_elements()); LOOP_VARNAME++) + +#define REF_DATA_ACCESS(ARRAY) (ARRAY).data()[LOOP_VARNAME] + +/* A view is non a contiguous memory access and can be accessed only by index list */ +#define LOOP_OVER_ALL_VIEW_ELEMENTS(ARRAY) \ + Index<Dim> LOOP_VARNAME((ARRAY).shape()); \ + LOOP_VARNAME.setIndexToMinusOne(); \ + while((++LOOP_VARNAME)() != LOOP_VARNAME.maxId()) + +#define VIEW_DATA_ACCESS(ARRAY) (ARRAY).operator()(LOOP_VARNAME.ids()) + + +/* CONST IMPLEMENTATIONS */ +/* Common const implementation */ +#define COMMON_CONST_IMPLEMENTATION(TYPENAME,TEMPLATES) \ + /* shape related */ \ + TEMPLATES \ + Shape<Dim> TYPENAME::shape() const { \ + Shape<Dim> shape; \ + const std::size_t* extents = this->super::shape(); \ + for (std::size_t d = 0; d < Dim; d++) \ + shape[d] = static_cast<std::size_t>(extents[d]); \ + return shape; \ + } + +/* Loop dependant const implementation */ +#define LOOP_DEPENDENT_CONST_IMPLEMENTATION(TYPENAME,TEMPLATES,LOOP_OVER_ALL_ELEMENTS,DATA_ACCESS) \ + /* print data */ \ + TEMPLATES \ + const TYPENAME& TYPENAME::print(const std::string& name, std::ostream& os, \ + unsigned int precision, unsigned int width) const { \ + auto S = this->super::shape(); \ + std::size_t id = 0; \ + \ + os << name << " = ["; \ + if(Dim==1) { \ + for(std::size_t k=0; k<this->num_elements(); k++) { \ + T x = this->operator()(boost::array<std::size_t,1>{k}); \ + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; \ + } \ + } \ + else if(Dim==2) { \ + std::cout << std::endl; \ + for(std::size_t i=0; i<S[0]; i++) { \ + os << "\t["; \ + for(std::size_t j=0; j<S[1]; j++) { \ + T x = this->operator()(boost::array<std::size_t,2>{i,j}); \ + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; \ + } \ + os << "]" << std::endl; \ + } \ + } \ + else if(Dim==3) { \ + std::cout << std::endl; \ + for(std::size_t i=0; i<S[0]; i++) { \ + os << "\t[["; \ + for(std::size_t j=0; j<S[1]; j++) { \ + if(j>0) \ + os << "\t ["; \ + for(std::size_t k=0; k<S[2]; k++) { \ + T x = this->operator()(boost::array<std::size_t,3>{i,j,k}); \ + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; \ + } \ + if(j!=S[1]-1) \ + os << "]" << std::endl; \ + else \ + os << "]]," << std::endl; \ + } \ + } \ + } \ + else { \ + LOOP_OVER_ALL_ELEMENTS(*this) { \ + T x = DATA_ACCESS(*this); \ + os << std::fixed << std::showpos << std::setprecision(precision) << std::setw(width) << x << " "; \ + } \ + } \ + os << "];" << std::endl; \ + return *this; \ + } \ + \ + /* boolean array_view specific functions */ \ + TEMPLATES ENABLE_IF_BOOL(bool,NO_DEFAULT_TEMPLATES) TYPENAME::all() const { \ + LOOP_OVER_ALL_ELEMENTS(*this) { \ + const bool val = DATA_ACCESS(*this); \ + if(!val) \ + return false; \ + } \ + return true; \ + } \ + TEMPLATES ENABLE_IF_BOOL(bool,NO_DEFAULT_TEMPLATES) TYPENAME::any() const { \ + LOOP_OVER_ALL_ELEMENTS(*this) { \ + const bool val = DATA_ACCESS(*this); \ + if(val) \ + return true; \ + } \ + return false; \ + } \ + TEMPLATES ENABLE_IF_BOOL(bool,NO_DEFAULT_TEMPLATES) TYPENAME::none() const { \ + LOOP_OVER_ALL_ELEMENTS(*this) { \ + const bool val = DATA_ACCESS(*this); \ + if(val) \ + return false; \ + } \ + return true; \ + } + +/* Reference specific const implementation */ +#define CONST_REF_IMPL(TYPENAME,TEMPLATES) \ + TEMPLATES \ + ENABLE_IF_REAL(const T*,NO_DEFAULT_TEMPLATES) TYPENAME::rdata() const { \ + return this->data(); \ + } \ + TEMPLATES \ + ENABLE_IF_REAL(const typename fft::fftw_complex_type<T>::std_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asStdComplexData() const { \ + return reinterpret_cast<const typename fft::fftw_complex_type<T>::std_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_REAL(const typename fft::fftw_complex_type<T>::fftw_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asFftwComplexData() const { \ + return reinterpret_cast<const typename fft::fftw_complex_type<T>::fftw_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(const T*,NO_DEFAULT_TEMPLATES) TYPENAME::cdata() const { \ + return this->data(); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::std_type*,NO_DEFAULT_TEMPLATES) TYPENAME::std_cdata() const { \ + return reinterpret_cast<const typename fft::fftw_complex_type<T>::std_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::fftw_type*,NO_DEFAULT_TEMPLATES) TYPENAME::fftw_cdata() const { \ + return reinterpret_cast<const typename fft::fftw_complex_type<T>::fftw_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(const typename fft::fftw_complex_type<T>::value_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asRealData() const { \ + return reinterpret_cast<const typename fft::fftw_complex_type<T>::value_type*>(this->data()); \ + } + +/* View specific const implementation */ +#define CONST_VIEW_IMPL(TYPENAME,TEMPLATES) + +/* All reference const implementation */ +#define CONST_REF_IMPLEMENTATION(TYPENAME,TEMPLATES) \ + COMMON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) \ + LOOP_DEPENDENT_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES),LOOP_OVER_ALL_REF_ELEMENTS,REF_DATA_ACCESS) \ + CONST_REF_IMPL(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) + +/* All view const implementation */ +#define CONST_VIEW_IMPLEMENTATION(TYPENAME,TEMPLATES) \ + COMMON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) \ + LOOP_DEPENDENT_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS) \ + CONST_VIEW_IMPL(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) + + + + +/* NON CONST IMPLEMENTATIONS */ +/* Common non const implementation */ +#define COMMON_NON_CONST_IMPLEMENTATION(TYPENAME,TEMPLATES) + +/* Loop dependant non const implementation */ +#define LOOP_DEPENDENT_NON_CONST_IMPLEMENTATION(TYPENAME,TEMPLATES,LOOP_OVER_ALL_ELEMENTS,DATA_ACCESS) \ + /* Apply function to all elements */ \ + TEMPLATES \ + TYPENAME& TYPENAME::apply(const std::function<void(T&)>& func) { \ + LOOP_OVER_ALL_ELEMENTS(*this) { \ + func( DATA_ACCESS(*this) ); \ + } \ + return *this; \ + } \ + TEMPLATES \ + TYPENAME& TYPENAME::apply(const std::function<void(T&, const Index<Dim>&)>& func) { \ + LOOP_OVER_ALL_VIEW_ELEMENTS(*this) { \ + func( VIEW_DATA_ACCESS(*this), LOOP_VARNAME ); \ + } \ + return *this; \ + } \ + TEMPLATES \ + template <typename Functor, typename Arg0, typename... Args> \ + TYPENAME& TYPENAME::apply(const Functor& func, Arg0&& farg0, Args&&... fargs) { \ + LOOP_OVER_ALL_VIEW_ELEMENTS(*this) { \ + func( VIEW_DATA_ACCESS(*this), LOOP_VARNAME, std::forward<Arg0>(farg0), std::forward<Args>(fargs)... ); \ + } \ + return *this; \ + } + +/* Reference specific non const implementation */ +#define NON_CONST_REF_IMPL(TYPENAME,TEMPLATES) \ + TEMPLATES \ + ENABLE_IF_REAL(T*,NO_DEFAULT_TEMPLATES) TYPENAME::rdata() { \ + return this->data(); \ + } \ + TEMPLATES \ + ENABLE_IF_REAL(typename fft::fftw_complex_type<T>::std_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asStdComplexData() { \ + return reinterpret_cast<typename fft::fftw_complex_type<T>::std_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_REAL(typename fft::fftw_complex_type<T>::fftw_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asFftwComplexData() { \ + return reinterpret_cast<typename fft::fftw_complex_type<T>::fftw_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(T*,NO_DEFAULT_TEMPLATES) TYPENAME::cdata() { \ + return this->data(); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::std_type*,NO_DEFAULT_TEMPLATES) TYPENAME::std_cdata() { \ + return reinterpret_cast<typename fft::fftw_complex_type<T>::std_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::fftw_type*,NO_DEFAULT_TEMPLATES) TYPENAME::fftw_cdata() { \ + return reinterpret_cast<typename fft::fftw_complex_type<T>::fftw_type*>(this->data()); \ + } \ + TEMPLATES \ + ENABLE_IF_COMPLEX(typename fft::fftw_complex_type<T>::value_type*,NO_DEFAULT_TEMPLATES) TYPENAME::asRealData() { \ + return reinterpret_cast<typename fft::fftw_complex_type<T>::value_type*>(this->data()); \ + } + +/* View specific non const implementation */ +#define NON_CONST_VIEW_IMPL(TYPENAME,TEMPLATES) + +/* Reference specific non const implementation */ +#define NON_CONST_REF_IMPLEMENTATION(TYPENAME,TEMPLATES) \ + COMMON_NON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) \ + LOOP_DEPENDENT_NON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES),LOOP_OVER_ALL_REF_ELEMENTS,REF_DATA_ACCESS) \ + NON_CONST_REF_IMPL(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) + +/* View specific non const implementation */ +#define NON_CONST_VIEW_IMPLEMENTATION(TYPENAME,TEMPLATES) \ + COMMON_NON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) \ + LOOP_DEPENDENT_NON_CONST_IMPLEMENTATION(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS) \ + NON_CONST_VIEW_IMPL(SINGLE_ARG(TYPENAME),SINGLE_ARG(TEMPLATES)) + + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_DEFINES_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_ext.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_ext.h new file mode 100644 index 000000000..da2f8f696 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_ext.h @@ -0,0 +1,247 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_EXT_H +#define HYSOP_MULTI_ARRAY_EXT_H + +namespace hysop { + namespace data { + + /* distances */ + template <typename T, std::size_t Dim> + T distance_L1(const const_multi_array_ref<T,Dim> &lhs, const const_multi_array_ref<T,Dim> &rhs); + template <typename T, std::size_t Dim> + T distance_L2(const const_multi_array_ref<T,Dim> &lhs, const const_multi_array_ref<T,Dim> &rhs); + template <typename T, std::size_t Dim> + T distance_Linf(const const_multi_array_ref<T,Dim> &lhs, const const_multi_array_ref<T,Dim> &rhs); + + /* unary operators */ + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator+(const const_multi_array_ref<T,Dim>& arr); + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator-(const const_multi_array_ref<T,Dim>& arr); + template <std::size_t Dim> + multi_array<bool,Dim> operator~(const const_multi_array_ref<bool,Dim>& arr); + + /* elementwise arithmetic operations */ + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator-(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator+(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator*(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator/(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<T,Dim> operator%(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + + /* elementwise boolean operations */ + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator& (const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator| (const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator^ (const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + + + /* element wise ordering test - near equality test for floating point types, see utils/utils.h::areEqual<T> */ + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator==(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator!=(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator>=(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator<=(const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator> (const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + template <typename T, std::size_t Dim> + multi_array<bool,Dim> operator< (const const_multi_array_ref<T,Dim>& lhs, const const_multi_array_ref<T,Dim>& rhs); + + + + /* And all their view and rvalue references variants */ + /* Ref - Ref */ + /* Ref - View */ + /* View - Ref */ + /* View - View */ + /* */ + /* Rvalue - View */ + /* Rvalue - Ref */ + /* View - Rvalue */ + /* Ref - Rvalue */ + /* ... */ + + + /* Implementation */ + + #define CHECK_SHAPES(LHS,RHS) assert(LHS.shape() == RHS.shape()) + + #define BINARY_OP(TEMPLATES,T,R,RES,OPNAME,BINOP,LHS,RHS,LOOP_OVER_ALL_ELEMENTS,DATA_ACCESS,CREATE_BUFFER,FROM_BUFFER,RET_OP) \ + TEMPLATES \ + RES OPNAME(LHS lhs, RHS rhs) { \ + CHECK_SHAPES(lhs,rhs); \ + CREATE_BUFFER(FROM_BUFFER,R); \ + LOOP_OVER_ALL_ELEMENTS(lhs) { \ + const T& lhsVal = DATA_ACCESS(lhs); \ + const T& rhsVal = DATA_ACCESS(rhs); \ + BINOP; \ + } \ + return RET_OP(BUFFER_NAME); \ + } + + #define BINARY_OP_VIEW_VIEW(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + VIEW(T),VIEW(T),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + #define BINARY_OP_VIEW_REF(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + VIEW(T),REF(T),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + #define BINARY_OP_REF_VIEW(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + REF(T),VIEW(T),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + #define BINARY_OP_REF_REF(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + REF(T),REF(T),LOOP_OVER_ALL_REF_ELEMENTS,REF_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + #define BINARY_OP_REF_RVALUE(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + REF(T),RVALUE(T),LOOP_OVER_ALL_REF_ELEMENTS,REF_DATA_ACCESS,CREATE_BUFFER,rhs,RET_OP) + #define BINARY_OP_RVALUE_REF(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + RVALUE(T),REF(T),LOOP_OVER_ALL_REF_ELEMENTS,REF_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + #define BINARY_OP_VIEW_RVALUE(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + VIEW(T),RVALUE(T),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS,CREATE_BUFFER,rhs,RET_OP) + #define BINARY_OP_RVALUE_VIEW(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),\ + RVALUE(T),VIEW(T),LOOP_OVER_ALL_VIEW_ELEMENTS,VIEW_DATA_ACCESS,CREATE_BUFFER,lhs,RET_OP) + + #define LVALUE_BINARY_OPS(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP_VIEW_VIEW(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_VIEW_REF(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_REF_VIEW(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_REF_REF(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) + + #define RVALUE_BINARY_OPS(TEMPLATES,T,R,RET,OPNAME,BINOP,CREATE_BUFFER,RET_OP) \ + BINARY_OP_REF_RVALUE(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_RVALUE_REF(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_VIEW_RVALUE(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) \ + BINARY_OP_RVALUE_VIEW(SINGLE_ARG(TEMPLATES),SINGLE_ARG(T),SINGLE_ARG(R),SINGLE_ARG(RET),OPNAME,SINGLE_ARG(BINOP),CREATE_BUFFER,RET_OP) + + + /* Code generator */ + #define BUFFER_NAME macro_generated_local_buffer + #define CREATE_MA_BUFFER(FROM_MA,R) multi_array<R,Dim> BUFFER_NAME(FROM_MA.shape()) + #define CREATE_R_BUFFER(FROM_MA,R) R BUFFER_NAME = R(0) + #define NO_MA_BUFFER(FROM_MA,R) multi_array<R,Dim>& BUFFER_NAME = FROM_MA + + #define SIMPLE_BINOP(OP) DATA_ACCESS(BUFFER_NAME) = ((lhsVal) OP (rhsVal)) + #define SIMPLE_FBINOP(FOP) DATA_ACCESS(BUFFER_NAME) = FOP((lhsVal),(rhsVal)) + #define SUP_OR_EQUAL() DATA_ACCESS(BUFFER_NAME) = (((lhsVal) > (rhsVal)) || hysop::utils::areEqual<T>((rhsVal),(lhsVal))) + #define INF_OR_EQUAL() DATA_ACCESS(BUFFER_NAME) = (((lhsVal) < (rhsVal)) || hysop::utils::areEqual<T>((rhsVal),(lhsVal))) + #define IDENTITY(X) (X) + + #define T1 SINGLE_ARG(template<std::size_t Dim>) + #define T1bis SINGLE_ARG(template<std::size_t Dim, typename Allocator>) + #define T2 SINGLE_ARG(template<typename T, std::size_t Dim>) + #define T3 SINGLE_ARG(template<typename T, std::size_t Dim, typename Allocator>) + #define VIEW(T) SINGLE_ARG(const const_multi_array_view<T,Dim>&) + #define REF(T) SINGLE_ARG(const const_multi_array_ref<T,Dim>&) + #define RVALUE(T) SINGLE_ARG(multi_array<T,Dim,Allocator>&&) + #define LVALUE(T) SINGLE_ARG(multi_array<T,Dim,Allocator>) + #define DLVALUE(T) SINGLE_ARG(multi_array<T,Dim>) + + + /* distances */ + LVALUE_BINARY_OPS(T2,T,T,T,distance_L1,\ + const T val = std::abs<T>(rhsVal-lhsVal); BUFFER_NAME += val,\ + CREATE_R_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,T,T,distance_L2,\ + const T val = std::abs<T>(rhsVal-lhsVal); BUFFER_NAME += val*val,\ + CREATE_R_BUFFER,std::sqrt) + LVALUE_BINARY_OPS(T2,T,T,T,distance_Linf,\ + SINGLE_ARG(const T val = std::abs<T>(rhsVal-lhsVal); BUFFER_NAME = std::max<T>(BUFFER_NAME,val)),\ + CREATE_R_BUFFER,IDENTITY) + + + /* elementwise arithmetic operations */ + LVALUE_BINARY_OPS(T2,T,T,DLVALUE(T),operator+,SIMPLE_BINOP(+),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,T,DLVALUE(T),operator-,SIMPLE_BINOP(-),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,T,DLVALUE(T),operator*,SIMPLE_BINOP(*),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,T,DLVALUE(T),operator/,SIMPLE_BINOP(/),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,T,DLVALUE(T),operator%,SIMPLE_BINOP(%),CREATE_MA_BUFFER,IDENTITY) + + RVALUE_BINARY_OPS(T3,T,T,LVALUE(T),operator+,SIMPLE_BINOP(+),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T3,T,T,LVALUE(T),operator-,SIMPLE_BINOP(-),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T3,T,T,LVALUE(T),operator*,SIMPLE_BINOP(*),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T3,T,T,LVALUE(T),operator/,SIMPLE_BINOP(/),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T3,T,T,LVALUE(T),operator%,SIMPLE_BINOP(%),NO_MA_BUFFER,IDENTITY) + + /* elementwise boolean like operations */ + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator&,SIMPLE_BINOP(&),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator|,SIMPLE_BINOP(|),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator^,SIMPLE_BINOP(^),CREATE_MA_BUFFER,IDENTITY) + + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator&,SIMPLE_BINOP(&),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator|,SIMPLE_BINOP(|),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator^,SIMPLE_BINOP(^),NO_MA_BUFFER,IDENTITY) + + /* element wise ordering test - near equality test for floating point types, see utils/utils.h::areEqual<T> */ + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator==,SIMPLE_FBINOP(hysop::utils::areEqual<T>),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator!=,SIMPLE_FBINOP(hysop::utils::areNotEqual<T>),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator<,SIMPLE_BINOP(<),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator>,SIMPLE_BINOP(>),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator>=,SUP_OR_EQUAL(),CREATE_MA_BUFFER,IDENTITY) + LVALUE_BINARY_OPS(T2,T,bool,DLVALUE(bool),operator<=,INF_OR_EQUAL(),CREATE_MA_BUFFER,IDENTITY) + + /* Comparisson for booleans... */ + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator==,SIMPLE_BINOP(==),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator!=,SIMPLE_BINOP(!=),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator<,SIMPLE_BINOP(<),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator>,SIMPLE_BINOP(>),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator>=,SIMPLE_BINOP(>=),NO_MA_BUFFER,IDENTITY) + RVALUE_BINARY_OPS(T1bis,bool,bool,LVALUE(bool),operator<=,SIMPLE_BINOP(<=),NO_MA_BUFFER,IDENTITY) + + /* clean macros */ + #undef LVALUE_BINARY_OPS + #undef RVALUE_BINARY_OPS + + #undef BINARY_OP_RVALUE_VIEW + #undef BINARY_OP_RVALUE_REF + #undef BINARY_OP_VIEW_RVALUE + #undef BINARY_OP_REF_RVALUE + + #undef BINARY_OP_REF_REF + #undef BINARY_OP_REF_VIEW + #undef BINARY_OP_VIEW_REF + #undef BINARY_OP_VIEW_VIEW + + #undef IDENTITY + #undef SIMPLE_BINOP + #undef SIMPLE_FBINOP + #undef BINARY_OP + + #undef DLVALUE + #undef LVALUE + #undef RVALUE + #undef REF + #undef VIEW + #undef T3 + #undef T2 + #undef T1bis + #undef T1 + + #undef NO_MA_BUFFER + #undef CREATE_R_BUFFER + #undef CREATE_MA_BUFFER + #undef CHECK_SHAPES + #undef BUFFER_NAME + + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_EXT_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h new file mode 100644 index 000000000..870595368 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h @@ -0,0 +1,192 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_IMPL_H +#define HYSOP_MULTI_ARRAY_IMPL_H + +namespace hysop { + namespace data { + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim, typename Allocator> + class multi_array : public boost_multi_array<T,Dim,Allocator> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_multi_array<T,Dim,Allocator>; + public: + PUBLIC_CLASS_TYPES() + + public: + multi_array(const extents_gen<Dim>& extents = extents_gen<Dim>()); + multi_array(const Shape<Dim>& shape); + + multi_array(const multi_array& other); + multi_array(multi_array&& other); + + explicit multi_array(const array_ref& other); + explicit multi_array(const array_view& other); + explicit multi_array(const const_array_ref& other); + explicit multi_array(const const_array_view& other); + + explicit multi_array(const boost_multi_array<T,Dim,Allocator>& other); + explicit multi_array(const boost_multi_array_ref<T,Dim>& other); + explicit multi_array(const boost_multi_array_view<T,Dim>& other); + explicit multi_array(const boost_const_multi_array_ref<T,Dim>& other); + explicit multi_array(const boost_const_multi_array_view<T,Dim>& other); + explicit multi_array(boost_multi_array<T,Dim,Allocator>&& other); + + multi_array& operator=(const multi_array& other); + multi_array& operator=(const array_ref& ref); + multi_array& operator=(const array_view& view); + multi_array& operator=(const const_array_ref& ref); + multi_array& operator=(const const_array_view& view); + multi_array& operator=(multi_array&& other); + + operator array_ref(); + operator const_array_ref() const; + + public: + PUBLIC_CONST_REF_INTERFACE(SINGLE_ARG(multi_array<T,Dim,Allocator>)) + PUBLIC_NON_CONST_REF_INTERFACE(SINGLE_ARG(multi_array<T,Dim,Allocator>)) + + multi_array& reshape(const Shape<Dim>& shape); + + protected: + static extents_gen<Dim> shapeToExtents(const Shape<Dim> &shape); + }; + + + /* Implementation */ + + /* constructors */ + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const extents_gen<Dim>& extents): + boost_multi_array<T,Dim,Allocator>(extents) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const Shape<Dim>& shape): + boost_multi_array<T,Dim,Allocator>(shapeToExtents(shape)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const multi_array& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array<T,Dim,Allocator>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const array_view& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array_view<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const const_array_view& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_const_multi_array_view<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const array_ref& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_multi_array_ref<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const const_array_ref& other): + boost_multi_array<T,Dim,Allocator>(static_cast<const boost_const_multi_array_ref<T,Dim>&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(multi_array&& other): + boost_multi_array<T,Dim,Allocator>(static_cast<boost_multi_array<T,Dim,Allocator>&&>(other)) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_multi_array_view<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_const_multi_array_view<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_multi_array_ref<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(const boost_const_multi_array_ref<T,Dim>& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::multi_array(boost_multi_array<T,Dim,Allocator>&& other): + boost_multi_array<T,Dim,Allocator>(other) {} + + + /* operator = */ + /* cast obligatory to avoid shape() function aliasing */ + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const multi_array<T,Dim,Allocator>& other) { + this->reshape(other.shape()); + super::operator=(dynamic_cast<const boost_multi_array<T,Dim,Allocator>&>(other)); + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const array_view& other) { + this->reshape(other.shape()); + super::operator=(dynamic_cast<const boost_multi_array_view<T,Dim>&>(other)); + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const const_array_view& other) { + this->reshape(other.shape()); + super::operator=(dynamic_cast<const boost_const_multi_array_view<T,Dim>&>(other)); + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const array_ref& other) { + this->reshape(other.shape()); + super::operator=(dynamic_cast<const boost_multi_array_ref<T,Dim>&>(other)); + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(const const_array_ref& other) { + this->reshape(other.shape()); + super::operator=(dynamic_cast<const boost_const_multi_array_ref<T,Dim>&>(other)); + return *this; + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::operator=(multi_array&& other) { + super::operator=(other); + return *this; + } + + /* casting operators */ + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::operator multi_array_ref<T,Dim>() { + return static_cast<boost_multi_array_ref<T,Dim>>(*this); + } + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>::operator const_multi_array_ref<T,Dim>() const { + return static_cast<boost_const_multi_array_ref<T,Dim>>(*this); + } + + /* static members */ + template <typename T, std::size_t Dim, typename Allocator> + typename multi_array<T,Dim,Allocator>::template extents_gen<Dim> multi_array<T,Dim,Allocator>::shapeToExtents(const Shape<Dim> &shape) { + return utils::buildExtents(shape); + } + + /* multiarray const & non const reference implementation */ + CONST_REF_IMPLEMENTATION(SINGLE_ARG(multi_array<T,Dim,Allocator>), SINGLE_ARG(template <typename T, std::size_t Dim, typename Allocator>)) + NON_CONST_REF_IMPLEMENTATION(SINGLE_ARG(multi_array<T,Dim,Allocator>), SINGLE_ARG(template <typename T, std::size_t Dim, typename Allocator>)) + + /* multi array specific */ + template <typename T, std::size_t Dim, typename Allocator> + multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::reshape(const Shape<Dim>& shape) { + boost::array<int,Dim> extents; + for (std::size_t d = 0; d < Dim; d++) + extents[d] = static_cast<int>(shape[d]); + this->resize(extents); + return *this; + } + + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_IMPL_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_ref.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_ref.h new file mode 100644 index 000000000..804185136 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_ref.h @@ -0,0 +1,70 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_REF_H +#define HYSOP_MULTI_ARRAY_REF_H + +namespace hysop { + namespace data { + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim> + class multi_array_ref : public boost_multi_array_ref<T,Dim> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_multi_array_ref<T,Dim>; + public: + PUBLIC_CLASS_TYPES() + + public: + multi_array_ref(T* base=nullptr, const extents_gen<Dim>& ranges = hysop::utils::buildExtents(std::array<std::size_t,Dim>{0})); + multi_array_ref(const multi_array_ref<T,Dim>& ref) = default; + multi_array_ref& operator=(const multi_array_ref<T,Dim>& other) = default; + + multi_array_ref(const boost_multi_array_ref<T,Dim>& ref); + multi_array_ref& operator=(const boost_multi_array_ref<T,Dim>& other); + + + operator const_array_ref() const; + + public: + PUBLIC_CONST_REF_INTERFACE(SINGLE_ARG(multi_array_ref<T,Dim>)) + PUBLIC_NON_CONST_REF_INTERFACE(SINGLE_ARG(multi_array_ref<T,Dim>)) + }; + + + /* Implementation */ + template <typename T, std::size_t Dim> + multi_array_ref<T,Dim>::multi_array_ref(T* base, const extents_gen<Dim>& ranges): + super(base,ranges) { + } + + template <typename T, std::size_t Dim> + multi_array_ref<T,Dim>::multi_array_ref(const boost_multi_array_ref<T,Dim>& ref) : + super(ref) { + } + + template <typename T, std::size_t Dim> + multi_array_ref<T,Dim>& multi_array_ref<T,Dim>::operator=(const boost_multi_array_ref<T,Dim>& other) { + super::operator=(other); + return *this; + } + + template <typename T, std::size_t Dim> + multi_array_ref<T,Dim>::operator const_array_ref() const { + return static_cast<boost_const_multi_array_ref<T,Dim>>(*this); + } + + CONST_REF_IMPLEMENTATION(SINGLE_ARG(multi_array_ref<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + NON_CONST_REF_IMPLEMENTATION(SINGLE_ARG(multi_array_ref<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + + } /* end of namespace data */ +} /* end of namespace hysop */ + + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_REF_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/multi_array_view.h b/HySoP/src/hysop++/src/data/multi_array/multi_array_view.h new file mode 100644 index 000000000..75d4fc89b --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_view.h @@ -0,0 +1,66 @@ + +#ifndef HYSOP_MULTI_ARRAY_H +#include "data/multi_array/multi_array.h" +#else + +#ifndef HYSOP_MULTI_ARRAY_VIEW_H +#define HYSOP_MULTI_ARRAY_VIEW_H + +namespace hysop { + namespace data { + + /* class hysop::data::multi_array */ + template <typename T, std::size_t Dim> + class multi_array_view : public boost_multi_array_view<T,Dim> { + static_assert(Dim>0, "Dim cannot be zero !"); + + private: + using super = boost_multi_array_view<T,Dim>; + public: + PUBLIC_CLASS_TYPES() + + public: + multi_array_view(const multi_array_view<T,Dim>& view) = default; + multi_array_view& operator=(const multi_array_view<T,Dim>& other) = default; + + multi_array_view(const boost_multi_array_view<T,Dim>& view); + multi_array_view& operator=(const boost_multi_array_view<T,Dim>& other); + + operator const_array_view() const; + + public: + PUBLIC_CONST_VIEW_INTERFACE(SINGLE_ARG(multi_array_view<T,Dim>)) + PUBLIC_NON_CONST_VIEW_INTERFACE(SINGLE_ARG(multi_array_view<T,Dim>)) + }; + + + /* Implementation */ + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wmaybe-uninitialized" + template <typename T, std::size_t Dim> + multi_array_view<T,Dim>::multi_array_view(const boost_multi_array_view<T,Dim>& view) : + super(view) { + } +#pragma GCC diagnostic pop + + template <typename T, std::size_t Dim> + multi_array_view<T,Dim>& multi_array_view<T,Dim>::operator=(const boost_multi_array_view<T,Dim>& other) { + super::operator=(other); + return *this; + } + + template <typename T, std::size_t Dim> + multi_array_view<T,Dim>::operator const_array_view() const { + return static_cast<boost_const_multi_array_view<T,Dim>>(*this); + } + + CONST_VIEW_IMPLEMENTATION(SINGLE_ARG(multi_array_view<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + NON_CONST_VIEW_IMPLEMENTATION(SINGLE_ARG(multi_array_view<T,Dim>), SINGLE_ARG(template <typename T, std::size_t Dim>)) + + } /* end of namespace data */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_MULTI_ARRAY_VIEW_H */ + +#endif /* end of MULTI_ARRAY include guard */ diff --git a/HySoP/src/hysop++/src/data/multi_array/tags b/HySoP/src/hysop++/src/data/multi_array/tags new file mode 100644 index 000000000..33ed9e868 --- /dev/null +++ b/HySoP/src/hysop++/src/data/multi_array/tags @@ -0,0 +1,6 @@ +!_TAG_FILE_FORMAT 2 /extended format; --format=1 will not append ;" to lines/ +!_TAG_FILE_SORTED 1 /0=unsorted, 1=sorted, 2=foldcase/ +!_TAG_PROGRAM_AUTHOR Darren Hiebert /dhiebert@users.sourceforge.net/ +!_TAG_PROGRAM_NAME Exuberant Ctags // +!_TAG_PROGRAM_URL http://ctags.sourceforge.net /official site/ +!_TAG_PROGRAM_VERSION 5.9~svn20110310 // diff --git a/HySoP/src/hysop++/src/detail/index_seq.h b/HySoP/src/hysop++/src/detail/index_seq.h new file mode 100644 index 000000000..2db6ab5c5 --- /dev/null +++ b/HySoP/src/hysop++/src/detail/index_seq.h @@ -0,0 +1,43 @@ + +#ifndef HYSOP_INDEX_SEQ_H +#define HYSOP_INDEX_SEQ_H + +namespace hysop { + namespace detail { + + template <int...> + struct index_seq {}; + + template <int k, std::size_t d, int... I> + struct constant_seq_impl { + typedef typename constant_seq_impl<k,d-1,k,I...>::type type; + }; + + template <int k, int... I> + struct constant_seq_impl<k,0,I...> { + typedef index_seq<I...> type; + }; + + template <std::size_t count, int step, int current, int... I> + struct index_seq_impl { + typedef typename index_seq_impl<count-1,step,current+step,I...,current>::type type; + }; + template <int step, int current, int... I> + struct index_seq_impl<0,step,current,I...> { + typedef index_seq<I...> type; + }; + + + + template <std::size_t count, int i0=0, int step=1> + using index_seq_gen = typename index_seq_impl<count,step,i0>::type; + + template <int constant, std::size_t count> + using constant_seq_gen = typename constant_seq_impl<constant,count>::type; + } +} + + + +#endif /* end of include guard: HYSOP_INDEX_SEQ_H */ + diff --git a/HySoP/src/hysop++/src/domain/boundary.cpp b/HySoP/src/hysop++/src/domain/boundary.cpp new file mode 100644 index 000000000..4cb3f491b --- /dev/null +++ b/HySoP/src/hysop++/src/domain/boundary.cpp @@ -0,0 +1,26 @@ + +#include "domain/boundary.h" + +namespace hysop { + namespace domain { + + static const char* boundary_strings[6] = { + "NONE", + "HOMOGENEOUS_DIRICHLET", + "HOMOGENEOUS_NEUMANN", + "DIRICHLET", + "NEUMANN", + "PERIODIC" + }; + + const char* toStringBoundary(Boundary bd) { + return boundary_strings[bd+1]; + } + + std::ostream& operator<<(std::ostream& os, const Boundary& bd) { + os << toStringBoundary(bd); + return os; + } + + } +} diff --git a/HySoP/src/hysop++/src/domain/boundary.h b/HySoP/src/hysop++/src/domain/boundary.h new file mode 100644 index 000000000..d0aa55438 --- /dev/null +++ b/HySoP/src/hysop++/src/domain/boundary.h @@ -0,0 +1,23 @@ +#ifndef BOUNDARY_H +#define BOUNDARY_H + +#include <iostream> + +namespace hysop { + namespace domain { + + enum Boundary : int { + NONE = -1, + HOMOGENEOUS_DIRICHLET = 0, + HOMOGENEOUS_NEUMANN = 1, + DIRICHLET = 2, + NEUMANN = 3, + PERIODIC = 4 + }; + + const char* toStringBoundary(Boundary bd); + std::ostream& operator<<(std::ostream& os, const Boundary& bd); + } +} + +#endif /* end of include guard: BOUNDARY_H */ diff --git a/HySoP/src/hysop++/src/domain/domain.h b/HySoP/src/hysop++/src/domain/domain.h new file mode 100644 index 000000000..ff228c3d1 --- /dev/null +++ b/HySoP/src/hysop++/src/domain/domain.h @@ -0,0 +1,154 @@ + +#ifndef HYSOP_DOMAIN_H +#define HYSOP_DOMAIN_H + +#include "data/multi_array/multi_array.h" +#include "utils/types.h" +#include "utils/utils.h" +#include "fft/fftDomainConfiguration.h" + +namespace hysop { + namespace domain { + + template <typename T, std::size_t Dim> + class Domain; + + template <typename T, std::size_t Dim> + std::ostream& operator<< (std::ostream& os, const Domain<T,Dim>& domain); + + template <typename T, std::size_t Dim> + class Domain { + public: + using DomainSize = std::array<T, Dim>; + using SpaceStep = std::array<T, Dim>; + using SpaceVariable = std::array<T, Dim>; + + public: + Domain() : + m_shape{0}, m_dataShape{0}, m_leftDataOffset{0}, + m_domainConfig(), m_domainSize{0.0L}, m_spaceStep{0.0L}, m_data() {} + + Domain(const Shape<Dim>& p_shape, const DomainConfiguration<Dim>& p_domainConfig, const DomainSize& p_domainSize) : + m_shape{0}, m_dataShape{0}, m_leftDataOffset{0}, + m_domainConfig(p_domainConfig), m_domainSize(p_domainSize), + m_spaceStep{0.0L}, m_data() { + this->reshape(p_shape); + } + + Domain(Domain<T,Dim>&& other) = default; + Domain<T,Dim>& operator=(const Domain<T,Dim>& other) = default; + Domain<T,Dim>& operator=(Domain<T,Dim>&& other) = default; + + Domain<T,Dim>& reshape(const Shape<Dim>& p_domainShape) { + m_shape = p_domainShape; + m_domainConfig.getShapeConfiguration(m_shape, m_dataShape, m_leftDataOffset); + m_data.reshape(m_dataShape); + this->computeSpaceStep(); + return *this; + } + + Domain<T,Dim>& resize(const DomainSize& p_size) { + m_domainSize = p_size; + this->computeSpaceStep(); + return *this; + } + + Domain<T,Dim>& resetDomainConfiguration(const DomainConfiguration<Dim>& p_domainConfig) { + m_domainConfig = p_domainConfig; + this->reshape(m_shape); + return *this; + } + + Domain<T,Dim>& print(const std::string &p_name) { + m_data.print(p_name); + return *this; + } + + fft::FftDomainConfiguration<Dim> fftDomainConfiguration() const { + return fft::FftDomainConfiguration<Dim>(m_domainConfig); + } + + const Shape<Dim>& shape() const { return m_shape; } + const Shape<Dim>& dataShape() const { return m_dataShape; } + const SpaceStep& spaceStep() const { return m_spaceStep; } + const DomainSize& domainSize() const { return m_domainSize; } + const Offset<Dim>& leftDataOffset() const { return m_leftDataOffset; } + const DomainConfiguration<Dim> boundaryConfig() const { return m_domainConfig; } + + const hysop::multi_array<T,Dim>& data() const { return m_data; } + hysop::multi_array_ref<T,Dim> data() { return m_data; } + + + /* Apply f(X, fargs...) on the whole domain where X = [x_0, x_1, ..., x_{Dim-1}] is the space variable */ + /* The result of the Functor f should be convertible to domain real data type T */ + template <typename Functor, typename... Args> + Domain<T,Dim>& apply(const Functor& f, Args&&... fargs); + + T distance_L1(const Domain<T,Dim> &other) { + assert(this->dataShape() == other.dataShape()); + return hysop::data::distance_L1<T,Dim>(this->data(), other.data()); + } + T distance_L2(const Domain<T,Dim> &other) { + assert(this->dataShape() == other.dataShape()); + return hysop::data::distance_L2<T,Dim>(this->data(), other.data()); + } + T distance_Linf(const Domain<T,Dim> &other) { + assert(this->dataShape() == other.dataShape()); + return hysop::data::distance_Linf<T,Dim>(this->data(), other.data()); + } + std::tuple<T,T,T> distance(const Domain<T,Dim>& other) { + return std::tuple<T,T,T>(distance_L1(other), distance_L2(other), distance_Linf(other)); + } + + protected: + void computeSpaceStep() { + for (std::size_t d = 0; d < Dim; d++) { + //std::size_t N = ((m_domainConfig[d].first==domain::Boundary::PERIODIC && !m_domainConfig.includePeriodicBoundaries()) ? m_shape[d]-1 : m_shape[d]-1); + m_spaceStep[d] = m_domainSize[d] / (m_shape[d]-1); + } + } + + protected: + Shape<Dim> m_shape, m_dataShape; + Offset<Dim> m_leftDataOffset; + + DomainConfiguration<Dim> m_domainConfig; + DomainSize m_domainSize; + SpaceStep m_spaceStep; + + hysop::multi_array<T, Dim> m_data; + }; + + /* Apply f(X, args...) on the whole domain where X = [x_0, x_1, ..., x_{Dim-1}] is the space variable */ + template <typename T, std::size_t Dim> + template <typename Functor, typename... Args> + Domain<T,Dim>& Domain<T,Dim>::apply(const Functor& f, Args&&... fargs) { + hysop::Index<Dim> idx(m_dataShape); + std::array<T,Dim> X{0.0L}; + T* data = m_data.origin(); + for (std::size_t k = 0; k < m_data.num_elements(); k++) { + data[k] = static_cast<T>(f(X, std::forward<Args>(fargs)...)); + ++idx; + for (std::size_t d = 0; d < Dim; d++) + X[d] = (idx[d]+m_leftDataOffset[d])*m_spaceStep[d]; + } + return *this; + } + + template <typename T, std::size_t Dim> + std::ostream& operator<< (std::ostream& os, const Domain<T,Dim>& domain) { + os << "== Domain ==" << std::endl; + os << "\tShape : " << domain.shape() << std::endl; + os << "\tSize : " << domain.domainSize() << std::endl; + os << "\tSpaceStep : " << domain.spaceStep() << std::endl; + os << domain.boundaryConfig(); + os << "\tLeftDataOffset: " << domain.leftDataOffset() << std::endl; + os << "\tDataShape : " << domain.dataShape() << std::endl; + return os; + } + + } +} + +#endif /* end of include guard: HYSOP_DOMAIN_H */ + diff --git a/HySoP/src/hysop++/src/domain/domainConfiguration.h b/HySoP/src/hysop++/src/domain/domainConfiguration.h new file mode 100644 index 000000000..42af207fa --- /dev/null +++ b/HySoP/src/hysop++/src/domain/domainConfiguration.h @@ -0,0 +1,98 @@ + +#ifndef HYSOP_DOMAIN_CONFIGURATION_H +#define HYSOP_DOMAIN_CONFIGURATION_H + +#include <array> +#include <stdexcept> + +#include "utils/types.h" +#include "domain/boundary.h" +#include "detail/index_seq.h" + +namespace hysop { + namespace domain { + + template <std::size_t Dim> + class DomainConfiguration; + + template <std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const DomainConfiguration<Dim>& config); + + template <std::size_t Dim> + class DomainConfiguration { + + public: + using BoundaryPair = std::pair<Boundary,Boundary>; + using BoundaryArray = std::array<BoundaryPair,Dim>; + public: + DomainConfiguration(const BoundaryArray& p_boundaries = defaultDomainBoundaries(), + //bool p_includeDirichletBoundaries = true, + bool p_includePeriodicBoundaries = true) : + m_boundaries(p_boundaries), + //m_includeDirichletBoundaries(p_includeDirichletBoundaries), + m_includePeriodicBoundaries(p_includePeriodicBoundaries) { + checkBoundaries(); + } + + const BoundaryArray& boundaries() const { return m_boundaries; } + //bool includeDirichletBoundaries() const { return m_includeDirichletBoundaries; } + bool includePeriodicBoundaries() const { return m_includePeriodicBoundaries; } + + const BoundaryPair& operator[](std::size_t k) const { return m_boundaries[k]; } + + void getShapeConfiguration(const Shape<Dim> &p_fullShape, Shape<Dim> &p_realShape, Offset<Dim> &p_leftOffset) const { + for (std::size_t d = 0; d < Dim; d++) { + const BoundaryPair& pair = m_boundaries[d]; + //bool hasDirichletLeftOffset = (pair.first == DIRICHLET || pair.first == HOMOGENEOUS_DIRICHLET); + //bool hasDirichletRightOffset = (pair.second == DIRICHLET || pair.second == HOMOGENEOUS_DIRICHLET); + //std::size_t dirichletLeftOffset = hasDirichletLeftOffset && !this->includeDirichletBoundaries(); + //std::size_t dirichletRightOffset = hasDirichletRightOffset && !this->includeDirichletBoundaries(); + bool hasPeriodicRightOffset = (pair.second == PERIODIC); + std::size_t periodicRightOffset = hasPeriodicRightOffset && !this->includePeriodicBoundaries(); + std::size_t leftOffset = 0; + std::size_t rightOffset = periodicRightOffset; + if(p_fullShape[d] <= (leftOffset + rightOffset)) + throw std::runtime_error("Domain shape is to small on axe " + std::to_string(d) + " for prescribed boundaries !"); + p_leftOffset[d] = std::ptrdiff_t(leftOffset); + p_realShape[d] = p_fullShape[d] - leftOffset - rightOffset; + } + } + + protected: + void checkBoundaries() const { + for (std::size_t d = 0; d < Dim; d++) { + const BoundaryPair& pair = m_boundaries[d]; + if((pair.first == PERIODIC) ^ (pair.second == PERIODIC)) + throw std::runtime_error("Bad boundaries configuration on axe " + std::to_string(d) + " !"); + } + } + + template <int... I> + static const std::array<BoundaryPair, Dim> defaultDomainBoundariesImpl(hysop::detail::index_seq<I...>) { + const BoundaryPair defaultVal[1] = { std::make_pair(Boundary::NONE,Boundary::NONE) }; + return { defaultVal[I]..., }; + } + static const std::array<BoundaryPair, Dim> defaultDomainBoundaries() { + return defaultDomainBoundariesImpl(hysop::detail::constant_seq_gen<0,Dim>()); + } + + protected: + BoundaryArray m_boundaries; + bool m_includePeriodicBoundaries; + //bool m_includeDirichletBoundaries; + }; + + template <std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const DomainConfiguration<Dim>& config) { + os << "== DomainConfiguration<Dim=" << std::to_string(Dim) << "> ==" << std::endl; + for (std::size_t d = 0; d < Dim; d++) + os << "\taxe[" << d << "]: " << config[d].first << "/" << config[d].second << std::endl; + //os << "\tDirichlet boundaries included ? " << std::boolalpha << config.includeDirichletBoundaries() << std::endl; + os << "\tPeriodic boundaries included ? " << std::boolalpha << config.includePeriodicBoundaries() << std::endl; + return os; + } + + } +} + +#endif /* end of include guard: HYSOP_DOMAIN_CONFIGURATION_H */ diff --git a/HySoP/src/hysop++/src/fft/extension.cpp b/HySoP/src/hysop++/src/fft/extension.cpp new file mode 100644 index 000000000..eeacb339b --- /dev/null +++ b/HySoP/src/hysop++/src/fft/extension.cpp @@ -0,0 +1,24 @@ + +#include "fft/extension.h" + +namespace hysop { + namespace fft { + + static const char* extension_strings[4] { + "NONE", + "EVEN", + "ODD", + "PERIODIC" + }; + + const char* toStringExtension(Extension ext) { + return extension_strings[ext+1]; + } + + std::ostream& operator<<(std::ostream& os, const Extension& ext) { + os << toStringExtension(ext); + return os; + } + + } +} diff --git a/HySoP/src/hysop++/src/fft/extension.h b/HySoP/src/hysop++/src/fft/extension.h new file mode 100644 index 000000000..815f25a65 --- /dev/null +++ b/HySoP/src/hysop++/src/fft/extension.h @@ -0,0 +1,23 @@ + +#ifndef FFT_EXTENSION_H +#define FFT_EXTENSION_H + +#include <ostream> + +namespace hysop { + namespace fft { + + enum Extension : int { + NONE=-1, + EVEN=0, + ODD=1, + PERIODIC=2 + }; + + const char* toStringExtension(Extension ext); + std::ostream& operator<<(std::ostream& os, const Extension& ext); + } +} + + +#endif /* end of include guard: FFT_EXTENSION_H */ diff --git a/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h b/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h new file mode 100644 index 000000000..f3ccea896 --- /dev/null +++ b/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h @@ -0,0 +1,114 @@ + +#ifndef HYSOP_FFTDOMAINCONFIGURATION_H +#define HYSOP_FFTDOMAINCONFIGURATION_H + +#include "utils/defines.h" +#include "domain/domainConfiguration.h" +#include "fft/extension.h" + +namespace hysop { + namespace fft { + + template <std::size_t Dim> + class FftDomainConfiguration; + + template <std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const FftDomainConfiguration<Dim>& config); + + template <std::size_t Dim> + class FftDomainConfiguration { + + public: + using BoundaryPair = typename domain::DomainConfiguration<Dim>::BoundaryPair; + using BoundaryArray = typename domain::DomainConfiguration<Dim>::BoundaryArray; + using ExtensionPair = std::pair<fft::Extension,fft::Extension>; + using ExtensionArray = std::array<ExtensionPair,Dim>; + public: + FftDomainConfiguration(const ExtensionArray& p_extensions, bool p_includePeriodicBoundaries) : + m_extensions(p_extensions), m_includePeriodicBoundaries(p_includePeriodicBoundaries) { + } + + FftDomainConfiguration(const domain::DomainConfiguration<Dim>& p_domain): + m_extensions(boundariesToExtensions(p_domain.boundaries())), + m_includePeriodicBoundaries(p_domain.includePeriodicBoundaries()) { + } + + const ExtensionArray& extensions() const { return m_extensions; } + bool includePeriodicBoundaries() const { return m_includePeriodicBoundaries; } + + ExtensionPair operator[](std::size_t k) const { + return m_extensions[k]; + } + + domain::DomainConfiguration<Dim> boundariesConfiguration() const { + return domain::DomainConfiguration<Dim>(fftExtensionsToBoundaries(m_extensions), m_includePeriodicBoundaries); + } + + static BoundaryArray fftExtensionsToBoundaries(const ExtensionArray& extArray) { + BoundaryArray bdArray; + for(std::size_t d=0; d<Dim; d++) { + const ExtensionPair& extPair = extArray[d]; + bdArray[d] = std::make_pair(fftExtensionToBoundary(extPair.first), fftExtensionToBoundary(extPair.second)); + } + return bdArray; + } + + static ExtensionArray boundariesToExtensions(const BoundaryArray& bdArray) { + ExtensionArray extArray; + for(std::size_t d=0; d<Dim; d++) { + const BoundaryPair& bdPair = bdArray[d]; + extArray[d] = std::make_pair(boundaryToExtension(bdPair.first), boundaryToExtension(bdPair.second)); + } + return extArray; + } + + static domain::Boundary fftExtensionToBoundary(fft::Extension ext) { + switch(ext) { + case(fft::Extension::PERIODIC): + return domain::Boundary::PERIODIC; + case(fft::Extension::EVEN): + return domain::Boundary::HOMOGENEOUS_DIRICHLET; + case(fft::Extension::ODD): + return domain::Boundary::HOMOGENEOUS_NEUMANN; + case(fft::Extension::NONE): + return domain::Boundary::NONE; + default: + NOT_IMPLEMENTED_YET; + } + } + + static fft::Extension boundaryToExtension(domain::Boundary bd) { + switch(bd) { + case(domain::Boundary::PERIODIC): + return fft::Extension::PERIODIC; + case(domain::Boundary::HOMOGENEOUS_DIRICHLET): + return fft::Extension::EVEN; + case(domain::Boundary::HOMOGENEOUS_NEUMANN): + return fft::Extension::ODD; + case(domain::Boundary::NONE): + return fft::Extension::NONE; + //throw std::runtime_error("Cannot build a FftDomainConfiguration based on a boundary of type 'NONE' !"); + default: + NOT_IMPLEMENTED_YET; + } + } + + protected: + ExtensionArray m_extensions; + bool m_includePeriodicBoundaries; + }; + + template <std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const FftDomainConfiguration<Dim>& config) { + os << "== FftDomainConfiguration<Dim=" << std::to_string(Dim) << "> ==" << std::endl; + for (std::size_t d = 0; d < Dim; d++) + os << "\taxe[" << d << "]: " << config[d].first << "/" << config[d].second; + os << "\tPeriodic boundaries included ?" << std::boolalpha << config.includePeriodicBoundaries() << std::endl; + return os; + } + + } +} + +#endif /* end of include guard: HYSOP_FFTDOMAINCONFIGURATION_H */ + diff --git a/HySoP/src/hysop++/src/fft/fftw3.h b/HySoP/src/hysop++/src/fft/fftw3.h new file mode 100644 index 000000000..6d214c5a4 --- /dev/null +++ b/HySoP/src/hysop++/src/fft/fftw3.h @@ -0,0 +1,440 @@ + +#ifndef HYSOP_FFTW3_H +#define HYSOP_FFTW3_H + +#include <complex> +#include <fftw3.h> + +#ifdef HYSOP_ENABLE_QUAD_MATH +#include <quadmath.h> +#endif + +/* */ +/* fftw3 c++ wrapper based on original macros in <fftw3.h> */ +/* */ + + +/* macros normally already defined in <fftw3.h> */ +#ifndef FFTW_CONCAT +#define FFTW_CONCAT(prefix, name) prefix ## name +#endif +#ifndef FFTW_MANGLE_FLOAT +#define FFTW_MANGLE_FLOAT(name) FFTW_CONCAT(fftwf_, name) +#endif +#ifndef FFTW_MANGLE_LONG_DOUBLE +#define FFTW_MANGLE_LONG_DOUBLE(name) FFTW_CONCAT(fftwl_, name) +#endif +#ifndef FFTW_MANGLE_QUAD +#define FFTW_MANGLE_QUAD(name) FFTW_CONCAT(fftwq_, name) +#endif +/***********************************************/ + +#undef FFTW_MANGLE_DOUBLE +#define FFTW_MANGLE_DOUBLE(name) FFTW_CONCAT(::fftw_, name) + +/* prefix for function wrappers inside the class */ +#define FFTW_MANGLE_CLASS(name) FFTW_CONCAT(fftw_, name) + +/* macro used to generate a full template specialisation of class Fftw3 for each type */ +#define FFTW_DEFINE_CXX_API(X, Y, REAL) \ +template <> \ + struct Fftw3<REAL> { \ + using R = REAL; \ + using C = Y(complex); \ + using plan = Y(plan); \ + using iodim = Y(iodim); \ + using iodim64 = Y(iodim64); \ + using r2r_kind = Y(r2r_kind); \ + using read_char_func = Y(read_char_func); \ + using write_char_func = Y(write_char_func); \ + using stdC = std::complex<REAL>; \ + \ + void X(execute)(const plan p) const { \ + Y(execute)(p); \ + } \ + \ + plan X(plan_dft)(int rank, const int *n, \ + C *in, C *out, int sign, unsigned int flags) const { \ + return Y(plan_dft)(rank, n, in, out, sign, flags); \ + } \ + \ + plan X(plan_dft_1d)(int n, C *in, C *out, int sign, \ + unsigned int flags) const { \ + return Y(plan_dft_1d)(n, in, out, sign, flags); \ + } \ + plan X(plan_dft_2d)(int n0, int n1, \ + C *in, C *out, int sign, unsigned int flags) const { \ + return Y(plan_dft_2d)(n0, n1, in, out, sign, flags); \ + } \ + plan X(plan_dft_3d)(int n0, int n1, int n2, \ + C *in, C *out, int sign, unsigned int flags) const { \ + return Y(plan_dft_3d)(n0, n1, n2, in, out, sign, flags); \ + } \ + \ + plan X(plan_many_dft)(int rank, const int *n, \ + int howmany, \ + C *in, const int *inembed, \ + int istride, int idist, \ + C *out, const int *onembed, \ + int ostride, int odist, \ + int sign, unsigned int flags) const { \ + return Y(plan_many_dft)(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, sign, flags); \ + } \ + \ + plan X(plan_guru_dft)(int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + C *in, C *out, \ + int sign, unsigned int flags) const { \ + return Y(plan_guru_dft)(rank, dims, howmany_rank, howmany_dims, in, out, sign, flags); \ + } \ + plan X(plan_guru_split_dft)(int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + R *ri, R *ii, R *ro, R *io, \ + unsigned int flags) const { \ + return Y(plan_guru_split_dft)(rank, dims, howmany_rank, howmany_dims, ri, ii, ro, io, flags); \ + } \ + \ + plan X(plan_guru64_dft)(int rank, \ + const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + C *in, C *out, \ + int sign, unsigned int flags) const { \ + return Y(plan_guru64_dft)(rank, dims, howmany_rank, howmany_dims, in, out, sign, flags); \ + } \ + plan X(plan_guru64_split_dft)(int rank, \ + const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + R *ri, R *ii, R *ro, R *io, \ + unsigned int flags) const { \ + return Y(plan_guru64_split_dft)(rank, dims, howmany_rank, howmany_dims, ri, ii, ro, io, flags); \ + } \ + \ + void X(execute_dft)(const plan p, C *in, C *out) const { \ + Y(execute_dft)(p, in, out); \ + } \ + void X(execute_split_dft)(const plan p, R *ri, R *ii, R *ro, R *io) const { \ + Y(execute_split_dft)(p, ri, ii, ro, io); \ + } \ + \ + plan X(plan_many_dft_r2c)(int rank, const int *n, \ + int howmany, \ + R *in, const int *inembed, \ + int istride, int idist, \ + C *out, const int *onembed, \ + int ostride, int odist, \ + unsigned int flags) const { \ + return Y(plan_many_dft_r2c)(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, flags); \ + } \ + \ + plan X(plan_dft_r2c)(int rank, const int *n, \ + R *in, C *out, unsigned int flags) const { \ + return Y(plan_dft_r2c)(rank, n, in, out, flags); \ + } \ + \ + plan X(plan_dft_r2c_1d)(int n,R *in,C *out,unsigned int flags) const { \ + return Y(plan_dft_r2c_1d)(n,in,out,flags); \ + } \ + plan X(plan_dft_r2c_2d)(int n0, int n1, \ + R *in, C *out, unsigned int flags) const { \ + return Y(plan_dft_r2c_2d)(n0, n1, in, out, flags); \ + } \ + plan X(plan_dft_r2c_3d)(int n0, int n1, \ + int n2, \ + R *in, C *out, unsigned int flags) const { \ + return Y(plan_dft_r2c_3d)(n0, n1, n2, in, out, flags); \ + } \ + \ + \ + plan X(plan_many_dft_c2r)(int rank, const int *n, \ + int howmany, \ + C *in, const int *inembed, \ + int istride, int idist, \ + R *out, const int *onembed, \ + int ostride, int odist, \ + unsigned int flags) const { \ + return Y(plan_many_dft_c2r)(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, flags); \ + } \ + \ + plan X(plan_dft_c2r)(int rank, const int *n, \ + C *in, R *out, unsigned int flags) const { \ + return Y(plan_dft_c2r)(rank, n, in, out, flags); \ + } \ + \ + plan X(plan_dft_c2r_1d)(int n,C *in,R *out,unsigned int flags) const { \ + return Y(plan_dft_c2r_1d)(n,in,out,flags); \ + } \ + plan X(plan_dft_c2r_2d)(int n0, int n1, \ + C *in, R *out, unsigned int flags) const { \ + return Y(plan_dft_c2r_2d)(n0, n1, in, out, flags); \ + } \ + plan X(plan_dft_c2r_3d)(int n0, int n1, \ + int n2, \ + C *in, R *out, unsigned int flags) const { \ + return Y(plan_dft_c2r_3d)(n0, n1, n2, in, out, flags); \ + } \ + \ + plan X(plan_guru_dft_r2c)(int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + R *in, C *out, \ + unsigned int flags) const { \ + return Y(plan_guru_dft_r2c)(rank, dims, howmany_rank, howmany_dims, in, out, flags); \ + } \ + plan X(plan_guru_dft_c2r)(int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + C *in, R *out, \ + unsigned int flags) const { \ + return Y(plan_guru_dft_c2r)(rank, dims, howmany_rank, howmany_dims, in, out, flags); \ + } \ + \ + plan X(plan_guru_split_dft_r2c)( \ + int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + R *in, R *ro, R *io, \ + unsigned int flags) const { \ + return Y(plan_guru_split_dft_r2c)( rank, dims, howmany_rank, howmany_dims, in, ro, io, flags); \ + } \ + plan X(plan_guru_split_dft_c2r)( \ + int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + R *ri, R *ii, R *out, \ + unsigned int flags) const { \ + return Y(plan_guru_split_dft_c2r)( rank, dims, howmany_rank, howmany_dims, ri, ii, out, flags); \ + } \ + \ + plan X(plan_guru64_dft_r2c)(int rank, \ + const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + R *in, C *out, \ + unsigned int flags) const { \ + return Y(plan_guru64_dft_r2c)(rank, dims, howmany_rank, howmany_dims, in, out, flags); \ + } \ + plan X(plan_guru64_dft_c2r)(int rank, \ + const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + C *in, R *out, \ + unsigned int flags) const { \ + return Y(plan_guru64_dft_c2r)(rank, dims, howmany_rank, howmany_dims, in, out, flags); \ + } \ + \ + plan X(plan_guru64_split_dft_r2c)( \ + int rank, const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + R *in, R *ro, R *io, \ + unsigned int flags) const { \ + return Y(plan_guru64_split_dft_r2c)( rank, dims, howmany_rank, howmany_dims, in, ro, io, flags); \ + } \ + plan X(plan_guru64_split_dft_c2r)( \ + int rank, const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + R *ri, R *ii, R *out, \ + unsigned int flags) const { \ + return Y(plan_guru64_split_dft_c2r)( rank, dims, howmany_rank, howmany_dims, ri, ii, out, flags); \ + } \ + \ + void X(execute_dft_r2c)(const plan p, R *in, C *out) const { \ + Y(execute_dft_r2c)(p, in, out); \ + } \ + void X(execute_dft_c2r)(const plan p, C *in, R *out) const { \ + Y(execute_dft_c2r)(p, in, out); \ + } \ + \ + void X(execute_split_dft_r2c)(const plan p, \ + R *in, R *ro, R *io) const { \ + return Y(execute_split_dft_r2c)(p, in, ro, io); \ + } \ + void X(execute_split_dft_c2r)(const plan p, \ + R *ri, R *ii, R *out) const { \ + return Y(execute_split_dft_c2r)(p, ri, ii, out); \ + } \ + \ + plan X(plan_many_r2r)(int rank, const int *n, \ + int howmany, \ + R *in, const int *inembed, \ + int istride, int idist, \ + R *out, const int *onembed, \ + int ostride, int odist, \ + const r2r_kind *kind, unsigned int flags) const { \ + return Y(plan_many_r2r)(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, kind, flags); \ + } \ + \ + plan X(plan_r2r)(int rank, const int *n, R *in, R *out, \ + const r2r_kind *kind, unsigned int flags) const { \ + return Y(plan_r2r)(rank, n, in, out, kind, flags); \ + } \ + \ + plan X(plan_r2r_1d)(int n, R *in, R *out, \ + r2r_kind kind, unsigned int flags) const { \ + return Y(plan_r2r_1d)(n, in, out, kind, flags); \ + } \ + plan X(plan_r2r_2d)(int n0, int n1, R *in, R *out, \ + r2r_kind kind0, r2r_kind kind1, \ + unsigned int flags) const { \ + return Y(plan_r2r_2d)(n0, n1, in, out, kind0, kind1, flags); \ + } \ + plan X(plan_r2r_3d)(int n0, int n1, int n2, \ + R *in, R *out, r2r_kind kind0, \ + r2r_kind kind1, r2r_kind kind2, \ + unsigned int flags) const { \ + return Y(plan_r2r_3d)(n0, n1, n2, in, out, kind0, kind1, kind2, flags); \ + } \ + \ + plan X(plan_guru_r2r)(int rank, const iodim *dims, \ + int howmany_rank, \ + const iodim *howmany_dims, \ + R *in, R *out, \ + const r2r_kind *kind, unsigned int flags) const { \ + return Y(plan_guru_r2r)(rank, dims, howmany_rank, howmany_dims, in, out, kind, flags); \ + } \ + \ + plan X(plan_guru64_r2r)(int rank, const iodim64 *dims, \ + int howmany_rank, \ + const iodim64 *howmany_dims, \ + R *in, R *out, \ + const r2r_kind *kind, unsigned int flags) const { \ + return Y(plan_guru64_r2r)(rank, dims, howmany_rank, howmany_dims, in, out, kind, flags); \ + } \ + \ + void X(execute_r2r)(const plan p, R *in, R *out) const { \ + Y(execute_r2r)(p, in, out); \ + } \ + \ + void X(destroy_plan)(plan p) const { \ + Y(destroy_plan)(p); \ + } \ + void X(forget_wisdom)() const { \ + Y(forget_wisdom)(); \ + } \ + void X(cleanup)() const { \ + Y(cleanup)(); \ + } \ + \ + void X(set_timelimit)(double t) const { \ + Y(set_timelimit)(t); \ + } \ + \ + void X(plan_with_nthreads)(int nthreads) const { \ + Y(plan_with_nthreads)(nthreads); \ + } \ + int X(init_threads)() const { \ + return Y(init_threads)(); \ + } \ + void X(cleanup_threads)() const { \ + Y(cleanup_threads)(); \ + } \ + \ + int X(export_wisdom_to_filename)(const char *filename) const { \ + return Y(export_wisdom_to_filename)(filename); \ + } \ + void X(export_wisdom_to_file)(FILE *output_file) const { \ + Y(export_wisdom_to_file)(output_file); \ + } \ + char *X(export_wisdom_to_string)() const { \ + return Y(export_wisdom_to_string)(); \ + } \ + void X(export_wisdom)(write_char_func write_char, \ + void *data) const { \ + return Y(export_wisdom)(write_char, data); \ + } \ + int X(import_system_wisdom)() const { \ + return Y(import_system_wisdom)(); \ + } \ + int X(import_wisdom_from_filename)(const char *filename) const { \ + return Y(import_wisdom_from_filename)(filename); \ + } \ + int X(import_wisdom_from_file)(FILE *input_file) const { \ + return Y(import_wisdom_from_file)(input_file); \ + } \ + int X(import_wisdom_from_string)(const char *input_string) const { \ + return Y(import_wisdom_from_string)(input_string); \ + } \ + int X(import_wisdom)(read_char_func read_char, void *data) const { \ + return Y(import_wisdom)(read_char, data); \ + } \ + \ + void X(fprint_plan)(const plan p, FILE *output_file) const { \ + Y(fprint_plan)(p, output_file); \ + } \ + void X(print_plan)(const plan p) const { \ + Y(print_plan)(p); \ + } \ + \ + void *X(malloc)(size_t n) const { \ + return Y(malloc)(n); \ + } \ + R *X(alloc_real)(size_t n) const { \ + return Y(alloc_real)(n); \ + } \ + C *X(alloc_complex)(size_t n) const { \ + return Y(alloc_complex)(n); \ + } \ + void X(free)(void *p) const { \ + Y(free)(p); \ + } \ + \ + void X(flops)(const plan p, \ + double *add, double *mul, double *fmas) const { \ + return Y(flops)(p, add, mul, fmas); \ + } \ + double X(estimate_cost)(const plan p) const { \ + return Y(estimate_cost)(p); \ + } \ + double X(cost)(const plan p) const { \ + return Y(cost)(p); \ + } \ + }; + + +namespace hysop { + namespace fft { + + template <typename T> struct is_fftw_supported_type { static constexpr bool value = false; }; + template <> struct is_fftw_supported_type<float> { static constexpr bool value = true; }; + template <> struct is_fftw_supported_type<double> { static constexpr bool value = true; }; + template <> struct is_fftw_supported_type<long double> { static constexpr bool value = true; }; +#ifdef HYSOP_ENABLE_QUAD_MATH + template <> struct is_fftw_supported_type<__float128> { static constexpr bool value = true; }; +#endif + + template <typename T> + struct Fftw3 { + Fftw3() { + throw std::runtime_error( + "Can only use Fftw3 wrapper with types {float, double, long double, __float128} ! " + "Special __float128 type is enabled via the macro HYSOP_ENABLE_QUAD_MATH." + ); + } + }; + + /* Generate Ffftw<> template specialisations */ + /* (1) link with lib fftw3f (-lfftw3f) */ + /* (2) link with lib fftw3 (-lfftw3) */ + /* (3) link with lib fftw3l (-lfftw3l) */ + /* (4) link with lib fftw3q and quadmath (-lfftw3q -lquadmath) */ + + FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_FLOAT, float) + FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_DOUBLE, double) + FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_LONG_DOUBLE, long double) +#ifdef HYSOP_ENABLE_QUAD_MATH + FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_QUAD, __float128) +#endif + + } +} + +#undef FFTW_DEFINE_CXX_API +#undef FFTW_MANGLE_CLASS + +#endif /* end of include guard: HYSOP_FFTW3_H */ + diff --git a/HySoP/src/hysop++/src/fft/fftwComplex.h b/HySoP/src/hysop++/src/fft/fftwComplex.h new file mode 100644 index 000000000..711be112d --- /dev/null +++ b/HySoP/src/hysop++/src/fft/fftwComplex.h @@ -0,0 +1,69 @@ + +#ifndef FFTWCOMPLEX_H +#define FFTWCOMPLEX_H + +#include <type_traits> +#include <complex> +#include <fftw3.h> + +namespace hysop { + namespace fft { + + template <typename T> struct fftw_complex_type { + typedef void value_type; + typedef void fftw_type; + typedef void std_type; + }; + + template <> struct fftw_complex_type<long double> { + typedef long double value_type; + typedef fftwl_complex fftw_type; + typedef std::complex<long double> std_type; + }; + template <> struct fftw_complex_type<double> { + typedef double value_type; + typedef fftw_complex fftw_type; + typedef std::complex<double> std_type; + }; + template <> struct fftw_complex_type<float> { + typedef float value_type; + typedef fftwf_complex fftw_type; + typedef std::complex<float> std_type; + }; + + template <> struct fftw_complex_type<std::complex<long double>> { + typedef long double value_type; + typedef fftwl_complex fftw_type; + typedef std::complex<long double> std_type; + }; + template <> struct fftw_complex_type<std::complex<double>> { + typedef double value_type; + typedef fftw_complex fftw_type; + typedef std::complex<double> std_type; + }; + template <> struct fftw_complex_type<std::complex<float>> { + typedef float value_type; + typedef fftwf_complex fftw_type; + typedef std::complex<float> std_type; + }; + + template <> struct fftw_complex_type<fftwl_complex> { + typedef long double value_type; + typedef fftwl_complex fftw_type; + typedef std::complex<long double> std_type; + }; + template <> struct fftw_complex_type<fftw_complex> { + typedef double value_type; + typedef fftw_complex fftw_type; + typedef std::complex<double> std_type; + }; + template <> struct fftw_complex_type<fftwf_complex> { + typedef float value_type; + typedef fftwf_complex fftw_type; + typedef std::complex<float> std_type; + }; + + } +} + +#endif /* end of include guard: FFTWCOMPLEX_H */ diff --git a/HySoP/src/hysop++/src/fft/planner.h b/HySoP/src/hysop++/src/fft/planner.h new file mode 100644 index 000000000..c80197cae --- /dev/null +++ b/HySoP/src/hysop++/src/fft/planner.h @@ -0,0 +1,745 @@ + +#ifndef HYSOP_PLANNER_H +#define HYSOP_PLANNER_H + +#include <list> +#include <cassert> +#include <functional> + +#include "fft/fftw3.h" +#include "fft/transform.h" +#include "fft/fftDomainConfiguration.h" +#include "domain/domain.h" +#include "utils/default.h" + +namespace hysop { + namespace fft { + + enum FftTransformType : int { + FFT_NONE=-1, + FFT_R2R, + FFT_R2C + }; + + template <typename T, std::size_t Dim> + class Planner; + + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const Planner<T,Dim>& planner); + + template <typename T, std::size_t Dim> + class Planner : private Fftw3<T> { + static_assert(hysop::fft::is_fftw_supported_type<T>::value, + "Planner data type is not currently supported by fftw !"); + + public: + using real = typename Fftw3<T>::R; + using fftw_complex = typename Fftw3<T>::C; + using std_complex = typename Fftw3<T>::stdC; + using fftw_plan = typename Fftw3<T>::plan; + using fftw_iodim = typename Fftw3<T>::iodim; + + protected: + using RealArray = hysop::multi_array<T,Dim>; + using TransformArray = std::array<fft::Transform,Dim>; + + public: + Planner(std::complex<T> p_fixedAxeWaveNumbers = std::complex<T>(1,0), + std::complex<T> p_fixedAxeWaveNumberPows = std::complex<T>(1,0)); + Planner(const Planner& other) = delete; /* cannot copy fftw plans to prevent inter instance plan destroying */ + Planner(Planner&& other) = default; + virtual ~Planner(); + + Planner& operator=(const Planner& planner) = delete; + Planner& operator=(Planner&& planner) = default; + + /* Plans a forward and a backward DCT/DST/FFT on each axe of the domain. */ + /* An axe with a NONE/NONE fft fomain configuration or order=0 is not transformed. */ + /* Input and output can be the same but in case of inplace transform, a temporary buffer is created */ + /* during planning. */ + /* If transforms include a complex transform (FFT) or if the transform is inplace, additional data may */ + /* be allocated inside the class. */ + /* fftw_flags are used to pass flags to FFTW. */ + /* Return true if planning was successfull AND if there is any transform, else return false. */ + /* Note: To get max performance use hysop::_default::fft_allocator<T> with hysop::multi_array */ + bool plan(hysop::const_multi_array_ref<T,Dim> input_rdata, + hysop::multi_array_ref<T,Dim> output_rdata, + const fft::FftDomainConfiguration<Dim>& inputFftDomainConfig, + const std::array<int,Dim>& order, + const std::array<T ,Dim>& domainSize, + unsigned int fftw_flags, + bool includeOutputPeriodicBoundaries=false, bool mirrorOutputPeriodicBoundaries=false); + + void executeForwardTransform(); + void executeBackwardTransform(); + + FftTransformType transformType() const; + + T normalisationFactor() const; + const std::array<T,Dim> signs() const; + const std::array<std::vector<std::complex<T>>,Dim>& waveNumbers() const; + const std::array<std::vector<std::complex<T>>,Dim>& waveNumbersPows() const; + + hysop::multi_array_view<T,Dim> transformedRealData(); /* view because possibly non contiguous */ + hysop::multi_array_ref<std::complex<T>,Dim> transformedComplexData(); /* ref because contiguous (allocated inside class) */ + + /* Get planned transform description */ + const std::string& toString() const; + + /* Modify untransformed axe generated wave numbers to simplify algorithm implementations based on the planner */ + Planner<T,Dim>& setFixedAxeWaveNumbers(std::complex<T> p_fixedAxeWaveNumber); + Planner<T,Dim>& setFixedAxeWaveNumberPows(std::complex<T> p_fixedAxeWaveNumberPow); + + protected: + void checkDomainCompatibility(const fft::FftDomainConfiguration<Dim>& domainConfig) const; + void checkExtensionCompatibility(const fft::Extension &lhs, const fft::Extension &rhs) const; + + fft::Transform findTransform(const std::pair<Extension,Extension>& ed) const; + + void reset(); + + template <typename Transfo> + void toStreamTransform(std::ostream& os, const Transfo& tr, + int rank, int howmany_rank, + const fftw_iodim* dims, const fftw_iodim *howmany_dims, + int input_data_offset=0, int output_data_offset=0); + + protected: + std::complex<T> m_fixedAxeWaveNumbers, m_fixedAxeWaveNumberPows; + + hysop::multi_array_ref<T,Dim> m_realBuffer; + hysop::multi_array<std_complex,Dim,hysop::_default::fft_allocator<std_complex>> m_complexBuffer; + typename hysop::multi_array<T,Dim>::template index_gen<Dim,Dim> m_transformedRealBufferView; + + std::list<fftw_plan> m_forward_R2R_plans, m_R2C_plans; + std::list<fftw_plan> m_backward_R2R_plans, m_C2R_plans; + + FftTransformType m_transformType; + std::vector<std::size_t> m_periodicAxes; + + std::array<T,Dim> m_signs; + std::array<std::vector<std_complex>,Dim> m_waveNumbers; + std::array<std::vector<std_complex>,Dim> m_waveNumbersPows; + T m_normalisationFactor; + + bool m_mirrorOutputPeriodicBoundaries; + std::string m_plannedTransformStr; + }; + + template <typename T, std::size_t Dim> + Planner<T,Dim>::Planner(std::complex<T> p_fixedAxeWaveNumbers, std::complex<T> p_fixedAxeWaveNumberPows): + m_fixedAxeWaveNumbers(p_fixedAxeWaveNumbers), + m_fixedAxeWaveNumberPows(p_fixedAxeWaveNumberPows), + m_realBuffer(), m_complexBuffer(), + m_forward_R2R_plans(), m_R2C_plans(), + m_backward_R2R_plans(), m_C2R_plans(), + m_transformType(), + m_signs(), m_waveNumbers(), m_waveNumbersPows(), m_normalisationFactor(), + m_mirrorOutputPeriodicBoundaries(), + m_plannedTransformStr() { + reset(); + } + + template <typename T, std::size_t Dim> + Planner<T,Dim>::~Planner() { + reset(); + } + + template <typename T, std::size_t Dim> + void Planner<T,Dim>::reset() { + new (&m_realBuffer) multi_array_ref<T,Dim>(); + m_complexBuffer.reshape(Shape<Dim>{0}); + for(fftw_plan plan : m_forward_R2R_plans) + this->fftw_destroy_plan(plan); + for(fftw_plan plan : m_backward_R2R_plans) + this->fftw_destroy_plan(plan); + for(fftw_plan plan : m_C2R_plans) + this->fftw_destroy_plan(plan); + for(fftw_plan plan : m_R2C_plans) + this->fftw_destroy_plan(plan); + m_forward_R2R_plans.clear(); + m_R2C_plans.clear(); + m_backward_R2R_plans.clear(); + m_C2R_plans.clear(); + m_periodicAxes.clear(); + m_transformType = FftTransformType::FFT_NONE; + for (std::size_t d = 0; d < Dim; d++) { + m_waveNumbers[d].clear(); + m_waveNumbersPows[d].clear(); + } + m_signs.fill(0); + m_normalisationFactor = T(0); + m_mirrorOutputPeriodicBoundaries = false; + m_plannedTransformStr = "** No planned transforms **"; + } + + template <typename T, std::size_t Dim> + void Planner<T,Dim>::checkDomainCompatibility(const fft::FftDomainConfiguration<Dim>& domainConfig) const { + for (std::size_t d = 0; d < Dim; d++) + checkExtensionCompatibility(domainConfig[d].first, domainConfig[d].second); + } + + template <typename T, std::size_t Dim> + void Planner<T,Dim>::checkExtensionCompatibility(const fft::Extension &lhs, const fft::Extension &rhs) const { + if ((lhs==fft::Extension::PERIODIC) ^ (rhs==fft::Extension::PERIODIC)) + throw std::runtime_error("Planner error: Periodic domain extensions are not compatible !"); + if ((lhs==fft::Extension::NONE) ^ (rhs==fft::Extension::NONE)) + throw std::runtime_error("Planner error: None domain extensions are not compatible !"); + } + + template <typename T, std::size_t Dim> + bool Planner<T,Dim>::plan(hysop::const_multi_array_ref<T,Dim> input_rdata, + hysop::multi_array_ref<T,Dim> output_rdata, + const fft::FftDomainConfiguration<Dim>& inputFftDomainConfig, + const std::array<int,Dim>& order, + const std::array<T ,Dim>& domainSize, + unsigned int fftw_flags, + const bool includeOutputPeriodicBoundaries, + const bool mirrorOutputPeriodicBoundaries) { + + this->checkDomainCompatibility(inputFftDomainConfig); + this->reset(); + + bool inplace = (input_rdata.data() == output_rdata.data()); + hysop::multi_array<T,Dim> input_rdata_buffer; + if(inplace) + input_rdata_buffer = input_rdata; + + m_mirrorOutputPeriodicBoundaries = includeOutputPeriodicBoundaries && mirrorOutputPeriodicBoundaries; + + TransformArray forwardTransforms, backwardTransforms; + Shape<Dim> realBufferShape, complexBufferShape; + std::array<int,Dim> forward_transform_size, backward_transform_size, complex_transform_size; + std::array<int,Dim> forward_input_offset, forward_output_offset; + std::array<int,Dim> backward_input_offset, backward_output_offset; + std::array<bool,Dim> oddOrder, axeTransformed; + + bool hasRealTransforms=false, hasComplexTransforms=false; + int lastComplexTransformAxe=-1; + std::stringstream ss; + + ss << "=== Planner ===" << std::endl; + ss << "** input **" << std::endl; + ss << "\tInputDataDim: " << input_rdata.shape() << std::endl; + ss << "\tDomainConfig: " << inputFftDomainConfig; + ss << "\tOrder : " << order << std::endl; + ss << "\tDomainSize : " << domainSize << std::endl; + ss << "\tipb/mpb : " << std::boolalpha << includeOutputPeriodicBoundaries << "/" + << std::boolalpha << mirrorOutputPeriodicBoundaries << std::endl; + + for (std::size_t d = 0; d < Dim; d++) { + const FftDomainConfiguration<Dim> &inputConfig = inputFftDomainConfig; + const auto &inputExtensions = inputConfig.extensions()[d]; + const bool axeIsPeriodic = (inputExtensions.first==Extension::PERIODIC); + const bool axeIsOddOrder = (order[d]%2!=0); + const bool axeIsTransformed = !(order[d]==0 || inputExtensions.first == fft::Extension::NONE); + const std::size_t inputSize = input_rdata.shape()[d] + - std::size_t(axeIsPeriodic && inputConfig.includePeriodicBoundaries()); + + const fft::Transform tr = this->findTransform(inputExtensions); + + forwardTransforms[d] = tr; + if(axeIsOddOrder) + backwardTransforms[d] = tr.conjugateInverseTransform(); + else + backwardTransforms[d] = tr.inverseTransform(); + + realBufferShape[d] = inputSize + + std::size_t((inputExtensions.second==Extension::PERIODIC) && includeOutputPeriodicBoundaries); + complexBufferShape[d] = inputSize; + forward_transform_size[d] = inputSize; + backward_transform_size[d] = inputSize; + complex_transform_size[d] = inputSize; + + forward_input_offset[d] = 0; + forward_output_offset[d] = 0; + backward_input_offset[d] = 0; + backward_output_offset[d] = 0; + + oddOrder[d] = axeIsOddOrder; + axeTransformed[d] = axeIsTransformed; + + if(axeIsPeriodic) + m_periodicAxes.push_back(d); + + if(!axeIsTransformed) { + continue; + } + else { + if(tr.isR2R()) { + hasRealTransforms = true; + const std::size_t firstExtOdd = (inputExtensions.first == ODD); + const std::size_t secondExtOdd = (inputExtensions.second == ODD); + const std::size_t firstExtEven = (inputExtensions.first == EVEN); + const std::size_t secondExtEven = (inputExtensions.second == EVEN); + + forward_transform_size[d] -= (firstExtOdd + secondExtOdd); + forward_input_offset[d] = std::ptrdiff_t(firstExtOdd); + forward_output_offset[d] = std::ptrdiff_t(firstExtOdd); + + complexBufferShape[d] = forward_transform_size[d]; + complex_transform_size[d] = forward_transform_size[d]; + + if(axeIsOddOrder) { + backward_transform_size[d] -= (firstExtEven + secondExtEven); + backward_input_offset[d] = std::ptrdiff_t( (tr.kind == FFTW_RODFT01) + || ((tr.kind != FFTW_REDFT01) && firstExtEven) ); + backward_output_offset[d] = std::ptrdiff_t(firstExtEven); + } + else { + backward_transform_size[d] = forward_transform_size[d]; + backward_input_offset[d] = forward_output_offset[d]; + backward_output_offset[d] = forward_input_offset[d]; + } + } + else { + hasComplexTransforms = true; + lastComplexTransformAxe = d; + } + } + } + + + /* Allocate complex buffer only if necessary, return if no transforms */ + if(hasComplexTransforms) { + m_transformType = FFT_R2C; + + assert(lastComplexTransformAxe >= 0); + complexBufferShape[lastComplexTransformAxe] = complexBufferShape[lastComplexTransformAxe]/2 + 1; + m_complexBuffer.reshape(complexBufferShape); + } + else if(hasRealTransforms) { + m_transformType = FFT_R2R; + } + else { + m_transformType = FFT_NONE; + if(!inplace) + std::copy(input_rdata.begin(), input_rdata.end(), output_rdata.begin()); + this->reset(); + return false; + } + + /* Check if output buffer shape match */ + { + if(output_rdata.shape() != realBufferShape) + throw std::runtime_error("Output buffer shape should match the planned one !"); + new (&m_realBuffer) multi_array_ref<T,Dim>(output_rdata); + } + + /* Compute normalisation factor and wave numbers */ + m_normalisationFactor = T(1); + for(std::size_t d=0; d<Dim; d++) { + const fft::Transform tr = forwardTransforms[d]; + const std::size_t N = (tr.isR2R() ? realBufferShape[d] : complexBufferShape[d]); + const std::size_t K = (tr.isR2R() ? forward_transform_size[d] : complexBufferShape[d]); + const std::size_t S = (tr.isR2R() ? forward_transform_size[d] : complex_transform_size[d]); + const real L = domainSize[d]; + T sign = T(1); + if(!axeTransformed[d]) { + m_waveNumbers[d].resize(N); + m_waveNumbersPows[d].resize(N); + m_normalisationFactor*=T(1); + for (std::size_t k = 0; k < N; k++) { + m_waveNumbers[d][k] = m_fixedAxeWaveNumbers; + m_waveNumbersPows[d][k] = m_fixedAxeWaveNumberPows; + } + } + else { + m_waveNumbers[d].resize(K); + m_waveNumbersPows[d].resize(K); + m_normalisationFactor*= tr.normalisation<T>(S); + if(tr.isR2R()) { + const std::size_t sign_offset = (tr.basefunc()==fft::BaseFunc::COS ? 1 : 0); + sign *= std::pow(T(-1),(order[d]+sign_offset)/2); + } + for (std::size_t k = 0; k < K; k++) { + m_waveNumbers[d][k] = tr.omega<T>(k,K,L,d==std::size_t(lastComplexTransformAxe)); + m_waveNumbersPows[d][k] = sign*std::pow(m_waveNumbers[d][k], order[d]); + } + } + m_signs[d] = sign; + } + + /* Build planned transforms description string */ + ss << "** configuration **" << std::endl; + ss << "\thasRealTransforms : " << std::boolalpha << hasRealTransforms << std::endl; + ss << "\thasComplexTransforms : " << std::boolalpha << hasComplexTransforms << std::endl; + ss << "\tTransformed axes : " << std::noboolalpha << axeTransformed << std::endl; + ss << "\tOddOrder : " << std::noboolalpha << oddOrder << std::endl; + ss << "\tForward transforms : " << forwardTransforms << std::endl; + ss << "\tBackward transforms : " << backwardTransforms << std::endl; + ss << "\tLast cplx trans. axe : " << lastComplexTransformAxe << std::endl; + + ss << "** buffer shape **" << std::endl; + ss << "\tReal buffer shape : " << realBufferShape << std::endl; + ss << "\tComplex buffer shape : " << complexBufferShape << std::endl; + + ss << "** transform size and offsets **" << std::endl; + ss << "\tForward transf. size : " << forward_transform_size << std::endl; + ss << "\tForward input offset: " << forward_input_offset << std::endl; + ss << "\tForward output offset: " << forward_output_offset << std::endl; + ss << "\tComplex transf. size : " << complex_transform_size << std::endl; + ss << "\tBackward transf. size : " << backward_transform_size << std::endl; + ss << "\tBackward input offset: " << backward_input_offset << std::endl; + ss << "\tBackward output offset: " << backward_output_offset << std::endl; + + ss << "** normalisation & wave numbers **" << std::endl; + ss << "\tNormalisation: " << m_normalisationFactor << std::endl; + ss << "\tSigns: " << m_signs << std::endl; + ss << "\t--Wave numbers--"<< std::endl; + for (std::size_t k = 0; k < Dim; k++) + ss << "\t\taxe" << k << ": " << m_waveNumbers[k] << std::endl; + ss << "\t--wave numbers powers--"<< std::endl; + for (std::size_t k = 0; k < Dim; k++) + ss << "\t\taxe" << k << ": " << m_waveNumbersPows[k] << std::endl; + + + /* Compute complex plans */ + if(hasComplexTransforms) { + ss << "** complex transforms detail **" << std::endl; + fftw_plan R2C_plan=NULL, C2R_plan=NULL; + const int rank = lastComplexTransformAxe+1; + const int howmany_rank = Dim; + fftw_iodim dims[Dim]; + fftw_iodim howmany_dims[Dim]; + int forwardInputOffset=0, backwardOutputOffset=0; + { + int local_io_stride = 1; + for(int d=Dim-1; d>=0; d--) { + forwardInputOffset += forward_output_offset[d]*local_io_stride; + backwardOutputOffset += forward_output_offset[d]*local_io_stride; + local_io_stride *= realBufferShape[d]; + } + } + { + int input_stride=1, output_stride=1; + for(int d=Dim-1, dd=rank; d>=0; d--) { + const bool isR2C = forwardTransforms[d].isR2C(); + const bool isAxeTransformed = axeTransformed[d]; + if(isR2C && isAxeTransformed) { + dims[d] = fftw_iodim{complex_transform_size[d], input_stride, output_stride}; + howmany_dims[d] = fftw_iodim{ 1, input_stride, output_stride}; + } + else { + dims[d] = fftw_iodim{ 1, input_stride, output_stride}; + howmany_dims[d] = fftw_iodim{forward_transform_size[d], input_stride, output_stride}; + } + input_stride *= realBufferShape[d]; + output_stride *= complexBufferShape[d]; + } + R2C_plan = this->fftw_plan_guru_dft_r2c( + rank, dims, + howmany_rank, howmany_dims, + m_realBuffer.rdata()+forwardInputOffset, m_complexBuffer.fftw_cdata(), + fftw_flags); + this->toStreamTransform(ss, "FFTW_FORWARD", rank, howmany_rank, dims, howmany_dims, + forwardInputOffset, 0); + if(!R2C_plan) { + ss << "=> R2C plan creation FAILED !" << std::endl; + m_plannedTransformStr = ss.str(); + return false; + } + } + { + int input_stride=1, output_stride=1; + for(int d=Dim-1, dd=rank; d>=0; d--) { + const bool isR2C = forwardTransforms[d].isR2C(); + const bool isAxeTransformed = axeTransformed[d]; + if(isR2C && isAxeTransformed) { + dims[d] = fftw_iodim{complex_transform_size[d], input_stride, output_stride}; + howmany_dims[d] = fftw_iodim{ 1, input_stride, output_stride}; + } + else { + dims[d] = fftw_iodim{ 1, input_stride, output_stride}; + howmany_dims[d] = fftw_iodim{ forward_transform_size[d], input_stride, output_stride}; + } + input_stride *= complexBufferShape[d]; + output_stride *= realBufferShape[d]; + } + C2R_plan = this->fftw_plan_guru_dft_c2r( + rank, dims, + howmany_rank, howmany_dims, + m_complexBuffer.fftw_cdata(), m_realBuffer.rdata()+backwardOutputOffset, + fftw_flags); + this->toStreamTransform(ss, "FFTW_BACKWARD", rank, howmany_rank, dims, howmany_dims, + 0, backwardOutputOffset); + if(!C2R_plan) { + ss << "=> C2R plan creation FAILED !" << std::endl; + m_plannedTransformStr = ss.str(); + return false; + } + } + m_R2C_plans.push_back(R2C_plan); + m_C2R_plans.push_back(C2R_plan); + } + + /* Compute real plans */ + if(hasRealTransforms) { + ss << "** real transforms detail **" << std::endl; + const int rank = 1; + const int howmany_rank = Dim; + fftw_r2r_kind kind[rank]; + fftw_iodim dims[rank]; + fftw_iodim howmany_dims[howmany_rank]; + int io_stride = 1; + for(int k=Dim-1; k>=0; k--) { + const fft::Transform& ftr = forwardTransforms[k]; + const fft::Transform& btr = backwardTransforms[k]; + const bool isR2R = ftr.isR2R() && btr.isR2R(); + const bool isAxeTransformed = axeTransformed[k]; + //int forwardInputOffset=0, forwardOutputOffset=0, backwardInputOffset=0, backwardOutputOffset=0; + if(isR2R && isAxeTransformed) { + ss << "\tTRANSFORM (" << std::to_string(k) << "):" << std::endl; + fftw_plan forward_plan=NULL, backward_plan=NULL; + int local_io_stride = 1; + for(int d=Dim-1; d>=0; d--) { + //howmany_dims[d] = fftw_iodim{ (d==k ? 1 : forward_transform_size[d]), local_io_stride, local_io_stride }; + //forwardInputOffset += forward_input_offset[d] *local_io_stride; + //forwardOutputOffset += forward_output_offset[d] *local_io_stride; + //backwardInputOffset += backward_input_offset[d] *local_io_stride; + //backwardOutputOffset += backward_output_offset[d]*local_io_stride; + howmany_dims[d] = fftw_iodim{ (d==k ? 1 : int(realBufferShape[d])), local_io_stride, local_io_stride }; + local_io_stride *= realBufferShape[d]; + } + { + kind[0] = static_cast<fftw_r2r_kind>(ftr.kind); + dims[0] = { fftw_iodim{forward_transform_size[k], io_stride, io_stride} }; + forward_plan = this->fftw_plan_guru_r2r( + rank, dims, + howmany_rank, howmany_dims, + //m_realBuffer.rdata()+forwardInputOffset, + //m_realBuffer.rdata()+forwardOutputOffset, + m_realBuffer.rdata()+forward_input_offset[k]*io_stride, + m_realBuffer.rdata()+forward_output_offset[k]*io_stride, + kind, fftw_flags); + //this->toStreamTransform(ss, ftr, rank, howmany_rank, dims, howmany_dims, + //forwardInputOffset, forwardOutputOffset); + this->toStreamTransform(ss, ftr, rank, howmany_rank, dims, howmany_dims, + forward_input_offset[k]*io_stride, forward_output_offset[k]*io_stride); + if(!forward_plan) { + ss << "=> Forward R2R plan creation FAILED !" << std::endl; + m_plannedTransformStr = ss.str(); + return false; + } + } + { + kind[0] = static_cast<fftw_r2r_kind>(btr.kind); + dims[0] = { fftw_iodim{backward_transform_size[k], io_stride, io_stride} }; + backward_plan = this->fftw_plan_guru_r2r( + rank, dims, + howmany_rank, howmany_dims, + //m_realBuffer.rdata()+backwardInputOffset, + //m_realBuffer.rdata()+backwardOutputOffset, + m_realBuffer.rdata()+backward_input_offset[k]*io_stride, + m_realBuffer.rdata()+backward_output_offset[k]*io_stride, + kind, fftw_flags); + //this->toStreamTransform(ss, btr, rank, howmany_rank, dims, howmany_dims, + //backwardInputOffset, backwardOutputOffset); + this->toStreamTransform(ss, btr, rank, howmany_rank, dims, howmany_dims, + backward_input_offset[k]*io_stride, backward_output_offset[k]*io_stride); + if(!backward_plan) { + ss << "=> Backward R2R plan creation FAILED !" << std::endl; + m_plannedTransformStr = ss.str(); + return false; + } + } + m_forward_R2R_plans.push_back(forward_plan); + m_backward_R2R_plans.push_back(backward_plan); + } + io_stride *= realBufferShape[k]; + } + } + + /* Copy input data into the buffer */ + if(inputFftDomainConfig.includePeriodicBoundaries() ^ includeOutputPeriodicBoundaries) { + NOT_IMPLEMENTED_YET; + } + else { + if(inplace) + std::copy(input_rdata_buffer.begin(), input_rdata_buffer.end(), m_realBuffer.begin()); + else + std::copy(input_rdata.begin() , input_rdata.end() , m_realBuffer.begin()); + } + + /* Create real buffer subview */ + ss << "** real buffer view **" << std::endl; + { + ss << "\tview = index_gen"; + std::array<boost::multi_array_types::index_range,Dim> ranges; + for (std::size_t d = 0; d < Dim; d++) { + const int offset = forward_output_offset[d]; + ranges[d] = boost::multi_array_types::index_range( + offset, + offset+forward_transform_size[d]); + + ss << "[range(" << ranges[d].start() << "," << ranges[d].finish() << ")]"; + } + ss << std::endl; + m_transformedRealBufferView = hysop::utils::buildIndices<Dim>(ranges); + } + + ss << "===============" << std::endl; + m_plannedTransformStr = ss.str(); + + return true; + } + + template <typename T, std::size_t Dim> + void Planner<T,Dim>::executeForwardTransform() { + for(auto& plan : m_forward_R2R_plans) + this->fftw_execute(plan); + for(auto& plan : m_R2C_plans) + this->fftw_execute(plan); + } + + template <typename T, std::size_t Dim> + void Planner<T,Dim>::executeBackwardTransform() { + for(auto& plan : m_C2R_plans) + this->fftw_execute(plan); + for(auto plan = m_backward_R2R_plans.rbegin(); plan!=m_backward_R2R_plans.rend(); ++plan) + this->fftw_execute(*plan); + + if(m_mirrorOutputPeriodicBoundaries) { + const Shape<Dim> rshape = m_realBuffer.shape(); + T *data = m_realBuffer.rdata(); + for(const int axe : m_periodicAxes) { + const int N = rshape[axe]; + const int num_elem = m_realBuffer.num_elements()/N; + std::array<int,Dim> ids{0}; + int id = 0; + int offset=1; + { + for(int d=Dim-1; d>axe; d--) offset*= rshape[d]; + offset *= (N-1); + } + for(int i=0; i<num_elem; i++) { + data[id+offset]=data[id]; + for (int d=Dim-1; d>=0; d--) { + if(d==axe) + continue; + ids[d]++; + if(ids[d]==int(rshape[d])) + ids[d]=0; + else + break; + } + id = ids[0]; + for (std::size_t d=1; d < Dim; d++) + id = id*rshape[d] + ids[d]; + } + } + } + } + + template <typename T, std::size_t Dim> + FftTransformType Planner<T,Dim>::transformType() const { + return m_transformType; + } + + template <typename T, std::size_t Dim> + T Planner<T,Dim>::normalisationFactor() const { + return m_normalisationFactor; + } + + template <typename T, std::size_t Dim> + const std::array<T,Dim> Planner<T,Dim>::signs() const { + return m_signs; + } + + template <typename T, std::size_t Dim> + const std::array<std::vector<std::complex<T>>,Dim>& Planner<T,Dim>::waveNumbers() const { + return m_waveNumbers; + } + + template <typename T, std::size_t Dim> + const std::array<std::vector<std::complex<T>>,Dim>& Planner<T,Dim>::waveNumbersPows() const { + return m_waveNumbersPows; + } + + template <typename T, std::size_t Dim> + hysop::multi_array_view<T,Dim> Planner<T,Dim>::transformedRealData() { + if(m_transformType==FFT_R2C) + throw std::runtime_error("Requesting planner real data but planned transform is real to complex !"); + else if(m_transformType==FFT_NONE) + throw std::runtime_error("Requesting planner real data but there was no successfull planned transforms !"); + return m_realBuffer[m_transformedRealBufferView]; + } + + template <typename T, std::size_t Dim> + hysop::multi_array_ref<std::complex<T>,Dim> Planner<T,Dim>::transformedComplexData() { + if(m_transformType==FFT_R2R) + throw std::runtime_error("Requesting planner complex data but planned transform is real to real !"); + else if(m_transformType==FFT_NONE) + throw std::runtime_error("Requesting planner real data but there was no successfull planned transforms !"); + return m_complexBuffer; + } + + template <typename T, std::size_t Dim> + fft::Transform Planner<T,Dim>::findTransform(const std::pair<Extension,Extension>& ed) const { + if(ed.first == EVEN) { + if(ed.second == EVEN) + return fft::Transform(FFTW_REDFT00); + else + return fft::Transform(FFTW_REDFT01); + } + else if(ed.first == ODD) { + if(ed.second == EVEN) + return fft::Transform(FFTW_RODFT01); + else + return fft::Transform(FFTW_RODFT00); + } + else { + return fft::Transform(FFTW_FORWARD); + } + } + + template <typename T, std::size_t Dim> + template <typename Transfo> + void Planner<T,Dim>::toStreamTransform(std::ostream& os, const Transfo& tr, + int rank, int howmany_rank, + const fftw_iodim* dims, const fftw_iodim* howmany_dims, + const int input_data_offset, const int output_data_offset) { + os << "\t --" << tr << "--" << std::endl; + os << "\t\tdims[" << rank << "] = {" << std::endl; + for (int i = 0; i < rank-1; i++) + os << "\t\t " << dims[i] << "," << std::endl; + os << "\t\t " << dims[rank-1] << std::endl; + os << "\t\t};" << std::endl; + os << "\t\thowmany[" << howmany_rank << "] = {" << std::endl; + for (int i = 0; i < howmany_rank-1; i++) + os << "\t\t " << howmany_dims[i] << "," << std::endl; + os << "\t\t " << howmany_dims[howmany_rank-1] << std::endl; + os << "\t\t};" << std::endl; + os << "\t\tinput data offset: " << input_data_offset << std::endl; + os << "\t\toutput data offset: " << output_data_offset << std::endl; + } + + /* Get planned transform description */ + template <typename T, std::size_t Dim> + const std::string& Planner<T,Dim>::toString() const { + return m_plannedTransformStr; + } + + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const Planner<T,Dim>& planner) { + os << planner.toString(); + return os; + } + + template <typename T, std::size_t Dim> + Planner<T,Dim>& Planner<T,Dim>::setFixedAxeWaveNumbers(std::complex<T> p_fixedAxeWaveNumber) { + m_fixedAxeWaveNumbers = p_fixedAxeWaveNumber; + return *this; + } + + template <typename T, std::size_t Dim> + Planner<T,Dim>& Planner<T,Dim>::setFixedAxeWaveNumberPows(std::complex<T> p_fixedAxeWaveNumberPow) { + m_fixedAxeWaveNumberPows = p_fixedAxeWaveNumberPow; + return *this; + } + + } /* end of namespace fft */ +} /* end of namespace hysop */ + + +#endif /* end of include guard: HYSOP_PLANNER_H */ diff --git a/HySoP/src/hysop++/src/fft/transform.cpp b/HySoP/src/hysop++/src/fft/transform.cpp new file mode 100644 index 000000000..53531952f --- /dev/null +++ b/HySoP/src/hysop++/src/fft/transform.cpp @@ -0,0 +1,14 @@ + +#include "fft/transform.h" + +namespace hysop { + namespace fft { + + std::ostream& operator<<(std::ostream& os, const Transform &tr) { + os << tr.toString(); + return os; + } + + } +} + diff --git a/HySoP/src/hysop++/src/fft/transform.h b/HySoP/src/hysop++/src/fft/transform.h new file mode 100644 index 000000000..06cc24333 --- /dev/null +++ b/HySoP/src/hysop++/src/fft/transform.h @@ -0,0 +1,188 @@ + +#ifndef FFTTRANSFORM_H +#define FFTTRANSFORM_H + +#include <complex> +#include <iostream> +#include <fftw3.h> + +#include "utils/constants.h" + +namespace hysop { + namespace fft { + + struct Transform; + std::ostream& operator<<(std::ostream& os, const Transform &tr); + + enum BaseFunc { + CEXP, + SIN, + COS + }; + + struct Transform { + public: + int kind; + + Transform(int p_kind = 0): + kind(p_kind) {} + + bool isR2C() const { + return kind == FFTW_FORWARD; + } + bool isR2R() const { + return !this->isR2C(); + } + + BaseFunc basefunc() const { + switch(kind) { + case(FFTW_REDFT00): + case(FFTW_REDFT01): + case(FFTW_REDFT10): + case(FFTW_REDFT11): + return COS; + case(FFTW_RODFT00): + case(FFTW_RODFT01): + case(FFTW_RODFT10): + case(FFTW_RODFT11): + return SIN; + case(FFTW_R2HC): + case(FFTW_HC2R): + case(FFTW_FORWARD): + return CEXP; + default: + throw std::runtime_error("Unknown transform !"); + } + } + + template <typename T> + std::complex<T> omega(std::size_t k, std::size_t N, T L = T(1.0L), bool lastDim=false) const { + using namespace hysop::constants; + switch(kind) { + case(FFTW_FORWARD): + if(lastDim) { + return std::complex<T>(T(0),T(2)*pi*T(k)/L); + } + else { + T kk; + if(k <= N/2 -1) + kk = T(k); + else if(k==N/2) + kk = T(0); + else + kk = double(k)-double(N); + return std::complex<T>(T(0),T(2)*pi*kk/L); + } + case(FFTW_REDFT00): + return std::complex<T>(pi*T(k)/L, T(0)); + case(FFTW_RODFT00): /* offset +1 */ + return std::complex<T>(pi*T(k+1)/L, T(0)); + case(FFTW_REDFT01): + return std::complex<T>(pi*(T(k)+T(0.5))/L, T(0)); + case(FFTW_RODFT01): /* -0.5 + 1 offset = +0.5 */ + return std::complex<T>(pi*(T(k)+T(0.5))/L, T(0)); + default: + throw std::runtime_error("Not implemented yet !"); + } + } + + template <typename T> + T normalisation(std::size_t n) const { + switch(kind) { + case(FFTW_FORWARD): + return T(n); + case(FFTW_RODFT00): + return T(2*(n+1)); + case(FFTW_REDFT00): + return T(2*(n-1)); + case(FFTW_REDFT01): + case(FFTW_REDFT10): + case(FFTW_REDFT11): + case(FFTW_RODFT01): + case(FFTW_RODFT10): + case(FFTW_RODFT11): + return T(2*n); + default: + return T(n); + } + } + + int inverseTransform() const { + switch(kind) { + case(FFTW_REDFT00): + case(FFTW_RODFT00): + return kind; + case(FFTW_REDFT01): + return FFTW_REDFT10; + case(FFTW_REDFT10): + return FFTW_REDFT01; + case(FFTW_RODFT01): + return FFTW_RODFT10; + case(FFTW_RODFT10): + return FFTW_RODFT01; + case(FFTW_R2HC): + return FFTW_HC2R; + case(FFTW_HC2R): + return FFTW_R2HC; + case(FFTW_FORWARD): + return FFTW_BACKWARD; + default: + throw std::runtime_error("Unknown transform !"); + } + } + + int conjugateInverseTransform() const { + switch(kind) { + case(FFTW_REDFT00): + return FFTW_RODFT00; + case(FFTW_REDFT01): + return FFTW_RODFT10; + case(FFTW_REDFT10): + return FFTW_RODFT01; + case(FFTW_RODFT00): + return FFTW_REDFT00; + case(FFTW_RODFT01): + return FFTW_REDFT10; + case(FFTW_RODFT10): + return FFTW_REDFT01; + default: + return this->inverseTransform(); + } + } + + std::string toString() const { + switch(kind) { + case(FFTW_REDFT00): + return "FFTW_REDFT00"; + case(FFTW_RODFT00): + return "FFTW_RODFT00"; + case(FFTW_REDFT01): + return "FFTW_REDFT01"; + case(FFTW_REDFT11): + return "FFTW_REDFT11"; + case(FFTW_REDFT10): + return "FFTW_REDFT10"; + case(FFTW_RODFT01): + return "FFTW_RODFT01"; + case(FFTW_RODFT10): + return "FFTW_RODFT10"; + case(FFTW_RODFT11): + return "FFTW_RODFT11"; + case(FFTW_R2HC): + return "FFTW_R2HC"; + //case(FFTW_HC2R): + //return "FFTW_HC2R"; + case(FFTW_BACKWARD): + return "FFTW_BACKWARD"; + case(FFTW_FORWARD): + return "FFTW_FORWARD"; + default: + return "FFTW_TRANSFORM_KIND_STRING_NOT_FOUND"; + } + } + }; + + } +} + +#endif /* end of include guard: FFTTRANSFORM_H */ diff --git a/HySoP/src/hysop++/src/maths/polynomial.h b/HySoP/src/hysop++/src/maths/polynomial.h new file mode 100644 index 000000000..d0bfaa977 --- /dev/null +++ b/HySoP/src/hysop++/src/maths/polynomial.h @@ -0,0 +1,427 @@ + +#ifndef HYSOP_POLYNOMIAL_H +#define HYSOP_POLYNOMIAL_H + +#include <array> +#include "data/multi_array/multi_array.h" + +namespace hysop { + namespace maths { + + /* Polynomials in dimension Dim with coefficients of type T */ + /* Basic polynomial operations are provided */ + /* TODO: Implement fast Nlog(N) multiplication by FFT */ + /* TODO: Implement polynomial division */ + template <typename T, std::size_t Dim> + class Polynomial; + + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const Polynomial<T,Dim>& poly); + + template <typename T, std::size_t Dim> + class Polynomial { + public: + /* constructors, destructors & operator= */ + Polynomial() = default; + Polynomial(const Polynomial& p_other) = default; + Polynomial( Polynomial&& p_other) = default; + Polynomial& operator=(const Polynomial& p_other) = default; + Polynomial& operator=( Polynomial&& p_other) = default; + ~Polynomial() = default; + + Polynomial(const std::array<std::size_t,Dim>& p_order); + explicit Polynomial(const hysop::multi_array<T,Dim>& p_coeffs); + explicit Polynomial( hysop::multi_array<T,Dim>&& p_coeffs); + + template <typename U> + explicit Polynomial(const Polynomial<U,Dim>& p_other); + + /* accessors */ + const hysop::multi_array<T,Dim>& coefficients() const; + hysop::multi_array<T,Dim>& coefficients(); + const std::array<std::size_t,Dim>& order() const; + std::array<std::size_t,Dim> shape() const; /* == order + 1 */ + + /* mutators */ + Polynomial& reshape(const std::array<std::size_t,Dim>& p_shape); + Polynomial& setOrder(const std::array<std::size_t,Dim>& p_order); + + Polynomial& applyToCoefficients(const std::function<void(T&)>& func); + Polynomial& applyToCoefficients(const std::function<void(T&, const Index<Dim>&)>& func); + + /* apply func(T&, const Index<Dim>&, farg0, fargs...) on all coefficients */ + template <typename Functor, typename Arg0, typename... Args> + Polynomial& applyToCoefficients(const Functor& func, Arg0&& farg0, Args&&... fargs); + + /* elementwise access to coefficients */ + const T& operator[](const Index<Dim> &p_id) const; + T& operator[](const Index<Dim> &p_id); + const T& operator[](std::size_t k) const; + T& operator[](std::size_t k); + + /* polynomial function evaluation with arbitrary type */ + template <typename U1, typename... U, typename R = typename std::common_type<T,U1,U...>::type> + R operator()(const U1& x1, const U&... xs) const; + template <typename U, typename R=typename std::common_type<T,U>::type> + R operator()(const std::array<U,Dim> &x) const; + template <typename U, typename R=typename std::common_type<T,U>::type> + R operator()(const U* x) const; + + /* basic elementwise operations */ + Polynomial& operator+=(const Polynomial& p_other); + Polynomial& operator-=(const Polynomial& p_other); + + Polynomial& operator*=(const T& p_val); + Polynomial& operator/=(const T& p_val); + Polynomial& operator%=(const T& p_val); + + /* polynomial multiplication and division */ + Polynomial& operator*=(const Polynomial& p_other); + Polynomial& operator/=(const Polynomial& p_other); + + /* integral and derivatives */ + Polynomial& integrate (std::size_t dim, int order); + Polynomial& differentiate(std::size_t dim, int order); + + template <typename I> typename std::enable_if<std::is_integral<I>::value, Polynomial&>::type + differentiate(const std::array<I,Dim>& order); + template <typename I> typename std::enable_if<std::is_integral<I>::value, Polynomial&>::type + integrate(const std::array<I,Dim>& order); + + template <typename I> typename std::enable_if<std::is_integral<I>::value, Polynomial&>::type + operator >>=(const std::array<I,Dim>& order); + template <typename I> typename std::enable_if<std::is_integral<I>::value, Polynomial&>::type + operator <<=(const std::array<I,Dim>& order); + + /* comparisson operators - uses near equality if T is a floating point type */ + bool operator==(const Polynomial& other); + bool operator!=(const Polynomial& other); + + /* misc */ + std::string toString(unsigned int p_precision=2, unsigned int p_width=6) const; + + protected: + /* misc */ + template <std::size_t D> + std::string toStringImpl(const T& p_coeff, unsigned int p_precision, unsigned int p_width, bool p_begin, bool p_end) const { + std::stringstream ss; + ss << (p_begin ? "" : " ") << std::fixed << std::showpos << std::setprecision(p_precision) << p_coeff; + return ss.str(); + } + template <std::size_t D, typename ArrayView, std::size_t K=Dim-D> + std::string toStringImpl(const ArrayView& p_view, unsigned int p_precision, unsigned int p_width, + bool=false, bool=false) const { + static const char varNames[3] = { 'z','y','x' }; + static const int offset = (Dim==1 ? 2 : (Dim==2 ? 1 : (Dim==3 ? 0 : -1))); + static const char delimiters[3][2] = { {'[',']'}, + {'{','}'}, + {'(',')'} }; + + + std::string str; + for (std::ptrdiff_t k=m_coeffs.shape()[K]-1; k>=0; k--) { + std::string localStr = toStringImpl<D-1>( + p_view[k], p_precision, p_width, k==std::ptrdiff_t(m_coeffs.shape()[K]-1), k==0); + if(localStr!="") { + if(D>1) + str += delimiters[D%3][0]; + str += localStr; + if(D>1) + str += delimiters[D%3][1]; + + std::string varName; + if(Dim<=3) + varName = varNames[K+offset]; + else + varName = "x_"+std::to_string(D); + if(k==0) + ; + else if(k==1) + str += varName; + else + str += varName + "^" + std::to_string(k); + if(k>0 && D>1) + str += " + "; + } + } + return str; + } + + /* static members */ + static std::array<std::size_t,Dim> orderFromShape(const std::array<std::size_t,Dim>& p_shape); + static std::array<std::size_t,Dim> shapeFromOrder(const std::array<std::size_t,Dim>& p_order); + + public: + template <typename X> + struct PolynomIndex : public Index<Dim> { + + public: + using typename Index<Dim>::Dimension; + using typename Index<Dim>::Indices; + public: + template <typename DimArray=typename Index<Dim>::Dimension, typename IndexArray=typename Index<Dim>::Indices> + PolynomIndex(std::array<X,Dim> p_spaceVar, + const DimArray& p_dim = Dimension{0}, + const IndexArray& p_ids = Indices{0}); + const std::array<X,Dim>& spaceVariable() const; /* returns {X[0],...,X[Dim-1]} */ + X value() const; /* returns X[0]^id[0] * X1^id[1] * ... * X[Dim-1]^id[Dim-1]*/ + protected: + void initialize(); + virtual void onIndexChange (std::size_t p_pos, std::ptrdiff_t p_offset) final override; + virtual void onIndexOverflow (std::size_t p_pos) final override; + + protected: + const std::array<X,Dim> m_spaceVar; + std::array<X,Dim> m_powers; + X m_value; + }; + + protected: + hysop::multi_array<T,Dim> m_coeffs; + std::array<std::size_t,Dim> m_order; + }; + + /* unary operations */ + template <typename T, std::size_t Dim> + Polynomial<T,Dim> operator+(const Polynomial<T,Dim>& poly); + template <typename T, std::size_t Dim> + Polynomial<T,Dim> operator-(const Polynomial<T,Dim>& poly); + + /* basic operations */ + template <typename T1, typename T2, std::size_t Dim, typename T = typename std::common_type<T1,T2>::type> + Polynomial<T,Dim> operator+(const Polynomial<T1,Dim>& lhs, const Polynomial<T2,Dim>& rhs); + template <typename T1, typename T2, std::size_t Dim, typename T = typename std::common_type<T1,T2>::type> + Polynomial<T,Dim> operator-(const Polynomial<T1,Dim>& lhs, const Polynomial<T2,Dim>& rhs); + template <typename T1, typename T2, std::size_t Dim, typename T = typename std::common_type<T1,T2>::type> + Polynomial<T,Dim> operator*(const Polynomial<T1,Dim>& lhs, const Polynomial<T2,Dim>& rhs); + template <typename T1, typename T2, std::size_t Dim, typename T = typename std::common_type<T1,T2>::type> + Polynomial<T,Dim> operator/(const Polynomial<T1,Dim>& lhs, const Polynomial<T2,Dim>& rhs); + + + /* tensor product of polynomials */ + template <typename T1, typename T2, std::size_t Dim1, std::size_t Dim2, + typename T=typename std::common_type<T1,T2>::type, std::size_t Dim=Dim1+Dim2> + Polynomial<T,Dim> operator|(const Polynomial<T1,Dim1>& lhs, const Polynomial<T2,Dim2>& rhs); + + + /* integral and derivatives */ + template <typename T, std::size_t Dim, typename I> + typename std::enable_if<std::is_integral<I>::value, Polynomial<T,Dim>>::type + operator<<(const Polynomial<T,Dim>& lhs, const std::array<I,Dim>& k); + template <typename T, std::size_t Dim, typename I> + typename std::enable_if<std::is_integral<I>::value, Polynomial<T,Dim>>::type + operator>>(const Polynomial<T,Dim>& lhs, const std::array<I,Dim>& k); + + template <typename T, std::size_t Dim, typename I> + typename std::enable_if<std::is_integral<I>::value, Polynomial<T,Dim>>::type + operator<<(Polynomial<T,Dim>&& lhs, const std::array<I,Dim>& k); + template <typename T, std::size_t Dim, typename I> + typename std::enable_if<std::is_integral<I>::value, Polynomial<T,Dim>>::type + operator>>(Polynomial<T,Dim>&& lhs, const std::array<I,Dim>& k); + + + + /********************/ + /** IMPLEMENTATION **/ + /********************/ + + /* static members */ + template <typename T, std::size_t Dim> + std::array<std::size_t,Dim> Polynomial<T,Dim>::orderFromShape(const std::array<std::size_t,Dim>& p_shape) { + std::array<std::size_t,Dim> order; + for (std::size_t d = 0; d < Dim; d++) + order[d] = p_shape[d]-1; + return order; + } + + template <typename T, std::size_t Dim> + std::array<std::size_t,Dim> Polynomial<T,Dim>::shapeFromOrder(const std::array<std::size_t,Dim>& p_order) { + std::array<std::size_t,Dim> shape; + for (std::size_t d = 0; d < Dim; d++) + shape[d] = p_order[d]+1; + return shape; + } + + /* constructors, destructors & operator= */ + template <typename T, std::size_t Dim> + Polynomial<T,Dim>::Polynomial(const std::array<std::size_t,Dim>& p_shape) : + m_coeffs(), m_order() { + this->reshape(p_shape); + } + + template <typename T, std::size_t Dim> + Polynomial<T,Dim>::Polynomial(const hysop::multi_array<T,Dim>& p_coeffs) : + m_coeffs(p_coeffs), m_order(orderFromShape(m_coeffs.shape())) { + } + + template <typename T, std::size_t Dim> + Polynomial<T,Dim>::Polynomial(hysop::multi_array<T,Dim>&& p_coeffs) : + m_coeffs(std::move(p_coeffs)), m_order(orderFromShape(m_coeffs.shape())) { + } + + template <typename T, std::size_t Dim> + template <typename U> + Polynomial<T,Dim>::Polynomial(const Polynomial<U,Dim>& p_other) { + this->reshape(p_other.shape()); + for (std::size_t k=0; k < m_coeffs.num_elements(); k++) + m_coeffs.data()[k] = static_cast<T>(p_other.data()[k]); + } + + /* accessors */ + template <typename T, std::size_t Dim> + const hysop::multi_array<T,Dim>& Polynomial<T,Dim>::coefficients() const { + return m_coeffs; + } + template <typename T, std::size_t Dim> + hysop::multi_array<T,Dim>& Polynomial<T,Dim>::coefficients() { + return m_coeffs; + } + template <typename T, std::size_t Dim> + const std::array<std::size_t,Dim>& Polynomial<T,Dim>::order() const { + return m_order; + } + template <typename T, std::size_t Dim> + std::array<std::size_t,Dim> Polynomial<T,Dim>::shape() const { + return m_coeffs.shape(); + } + + /* mutators */ + template <typename T, std::size_t Dim> + Polynomial<T,Dim>& Polynomial<T,Dim>::reshape(const std::array<std::size_t,Dim>& p_shape) { + m_order = orderFromShape(p_shape); + m_coeffs.reshape(p_shape); + return *this; + } + + template <typename T, std::size_t Dim> + Polynomial<T,Dim>& Polynomial<T,Dim>::setOrder(const std::array<std::size_t,Dim>& p_order) { + m_order = p_order; + m_coeffs.reshape(shapeFromOrder(p_order)); + return *this; + } + + + template <typename T, std::size_t Dim> + Polynomial<T,Dim>& Polynomial<T,Dim>::applyToCoefficients(const std::function<void(T&)>& func) { + m_coeffs.apply(func); + return *this; + } + template <typename T, std::size_t Dim> + Polynomial<T,Dim>& Polynomial<T,Dim>::applyToCoefficients(const std::function<void(T&, const Index<Dim>&)>& func) { + m_coeffs.apply(func); + return *this; + } + + /* apply func(T&, const Index<Dim>&, farg0, fargs...) on all coefficients */ + template <typename T, std::size_t Dim> + template <typename Functor, typename Arg0, typename... Args> + Polynomial<T,Dim>& Polynomial<T,Dim>::applyToCoefficients(const Functor& func, Arg0&& farg0, Args&&... fargs) { + m_coeffs.apply(func, farg0, fargs...); + return *this; + } + + /* access to coefficients */ + template <typename T, std::size_t Dim> + const T& Polynomial<T,Dim>::operator[](std::size_t k) const { + return m_coeffs.data()[k]; + } + template <typename T, std::size_t Dim> + T& Polynomial<T,Dim>::operator[](std::size_t k) { + return m_coeffs.data()[k]; + } + template <typename T, std::size_t Dim> + const T& Polynomial<T,Dim>::operator[](const Index<Dim> &p_id) const { + return m_coeffs.data()[p_id.id()]; + } + + template <typename T, std::size_t Dim> + T& Polynomial<T,Dim>::operator[](const Index<Dim> &p_id) { + return m_coeffs.data()[p_id.id()]; + } + + /* polynomial evaluation */ + template <typename T, std::size_t Dim> + template <typename U1, typename... U, typename R> + R Polynomial<T,Dim>::operator()(const U1& x1, const U&... xs) const { + return this->operator()(std::array<R,Dim>{x1,xs...}); + } + template <typename T, std::size_t Dim> + template <typename U, typename R> + R Polynomial<T,Dim>::operator()(const U* p_x) const { + return this->operator()(std::array<U,Dim>(p_x)); + } + template <typename T, std::size_t Dim> + template <typename U, typename R> + R Polynomial<T,Dim>::operator()(const std::array<U,Dim> &p_x) const { + /* compute result */ + R res = R(0); + const T* coeffs = m_coeffs.data(); + PolynomIndex<U> idx(p_x, this->shape()); + while(!idx.atMaxId()) { + res += coeffs[idx()]*idx.value(); + ++idx; + } + return res; + } + + template <typename T, std::size_t Dim> + std::string Polynomial<T,Dim>::toString(unsigned int p_precision, unsigned int p_width) const { + return toStringImpl<Dim>(m_coeffs, p_precision, p_width); + } + + /* struct PolynomIndex */ + template <typename T, std::size_t Dim> + template <typename X> + template <typename DimArray, typename IndexArray> + Polynomial<T,Dim>::PolynomIndex<X>:: + PolynomIndex(std::array<X,Dim> p_spaceVar, const DimArray& p_dim, const IndexArray& p_ids): + Index<Dim>(p_dim, p_ids), m_spaceVar(p_spaceVar), m_powers{0}, m_value(0) { + this->initialize(); + } + template <typename T, std::size_t Dim> + template <typename X> + void Polynomial<T,Dim>::PolynomIndex<X>::initialize() { + m_value = X(1); + for (std::size_t d=0; d<Dim; d++) { + X power = std::pow(m_spaceVar[d],this->operator[](d)); + m_powers[d] = power; + m_value *= power; + } + } + template <typename T, std::size_t Dim> + template <typename X> + const std::array<X,Dim>& Polynomial<T,Dim>::PolynomIndex<X>::spaceVariable() const { + return m_spaceVar; + } + template <typename T, std::size_t Dim> + template <typename X> + X Polynomial<T,Dim>::PolynomIndex<X>::value() const { + return m_value; + } + template <typename T, std::size_t Dim> + template <typename X> + void Polynomial<T,Dim>::PolynomIndex<X>::onIndexChange(std::size_t p_pos, std::ptrdiff_t p_offset) { + assert(p_offset == 1); + m_powers[p_pos] = m_powers[p_pos]*m_spaceVar[p_pos]; + m_value *= m_spaceVar[p_pos]; + } + template <typename T, std::size_t Dim> + template <typename X> + void Polynomial<T,Dim>::PolynomIndex<X>::onIndexOverflow (std::size_t p_pos) { + m_powers[p_pos] = X(1); + m_value = X(1); + for (std::size_t d=0; d < Dim; d++) + m_value *= m_powers[d]; + } + + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const Polynomial<T,Dim>& poly) { + os << poly.toString(); + return os; + } + + } /* end of namespace maths */ +} /* end of namesapce hysop */ + + +#endif /* end of include guard: HYSOP_POLYNOMIAL_H */ diff --git a/HySoP/src/hysop++/src/solver/diffSolver.h b/HySoP/src/hysop++/src/solver/diffSolver.h new file mode 100644 index 000000000..7b6908803 --- /dev/null +++ b/HySoP/src/hysop++/src/solver/diffSolver.h @@ -0,0 +1,29 @@ + +#ifndef HYSOP_DIFFSOLVER_H +#define HYSOP_DIFFSOLVER_H + +#include "data/multi_array/multi_array.h" + +namespace hysop { + namespace solver { + template <typename T, std::size_t Dim> + class DiffSolver { + + public: + virtual void apply(hysop::const_multi_array_ref<T,Dim> input, + hysop::multi_array_ref<T,Dim> output, + const std::array<int,Dim> &order) const = 0; + + void operator()(hysop::const_multi_array_ref<T,Dim> input, + hysop::multi_array_ref<T,Dim> output, + const std::array<int,Dim> &order) { + this->apply(input,output,order); + } + + }; + + } /* end of namespace solver */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_DIFFSOLVER_H */ + diff --git a/HySoP/src/hysop++/src/solver/fftDiffSolver.h b/HySoP/src/hysop++/src/solver/fftDiffSolver.h new file mode 100644 index 000000000..d9951e7a7 --- /dev/null +++ b/HySoP/src/hysop++/src/solver/fftDiffSolver.h @@ -0,0 +1,191 @@ + +#ifndef HYSOP_FFTDIFFSOLVER_H +#define HYSOP_FFTDIFFSOLVER_H + +#include "solver/diffSolver.h" +#include "fft/planner.h" +#include "data/accumulatorIndex.h" + +namespace hysop { + namespace solver { + + template <typename T, std::size_t Dim> + class FftDiffSolver : public DiffSolver<T,Dim> { + + private: + using super = DiffSolver<T,Dim>; + public: + FftDiffSolver() = default; + FftDiffSolver(const FftDiffSolver& p_other) = default; + FftDiffSolver( FftDiffSolver&& p_other) = default; + ~FftDiffSolver() = default; + + FftDiffSolver& operator=(const FftDiffSolver& p_other) = default; + FftDiffSolver& operator=( FftDiffSolver&& p_other) = default; + + FftDiffSolver(const std::array<T,Dim> &p_domainSize, const fft::FftDomainConfiguration<Dim> &p_inputFftConfig, + unsigned int p_fftFlags = FFTW_MEASURE, + bool p_includeOutputPeriodicBds=false, bool p_mirrorOutputPeriodicBds=false); + + /* Mutators */ + FftDiffSolver& setDomainSize(const std::array<T,Dim>& p_domainSize); + FftDiffSolver& setFftDomainConfiguration(const fft::FftDomainConfiguration<Dim>& p_fftConfig); + + FftDiffSolver& setFftFlags(unsigned int p_flags); + FftDiffSolver& appendFftFlags(unsigned int p_flags); + + FftDiffSolver& includeOutputPeriodicBoundaries(bool p_val = true); + FftDiffSolver& excludeOutputPeriodicBoundaries(); + + FftDiffSolver& enableOutputPeriodicBoundariesMirroring(bool p_val = true); + FftDiffSolver& disableOutputPeriodicBoundariesMirroring(); + + /* Accessors */ + std::array<T,Dim> domainSize() const; + fft::FftDomainConfiguration<Dim> inputFftConfig() const; + + unsigned int fftFlags() const; + bool includeOutputPeriodicBds() const; + bool mirrorOutputPeriodicBds() const; + + /* Apply operator */ + virtual void apply(hysop::const_multi_array_ref<T,Dim> p_input, + hysop::multi_array_ref<T,Dim> p_output, + const std::array<int,Dim> &p_order) const final override; + + protected: + std::array<T,Dim> m_domainSize; + fft::FftDomainConfiguration<Dim> m_inputFftConfig; + unsigned int m_fftFlags; + bool m_includeOutputPeriodicBds, m_mirrorOutputPeriodicBds; + }; + + + + /* Implementation */ + + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>::FftDiffSolver(const std::array<T,Dim> &p_domainSize, const fft::FftDomainConfiguration<Dim> &p_inputFftConfig, + unsigned int p_fftFlags, + bool p_includeOutputPeriodicBds, bool p_mirrorOutputPeriodicBds): + m_domainSize(p_domainSize), m_inputFftConfig(p_inputFftConfig), + m_fftFlags(p_fftFlags), + m_includeOutputPeriodicBds(p_includeOutputPeriodicBds), m_mirrorOutputPeriodicBds(p_mirrorOutputPeriodicBds) { + } + + + /* Mutators */ + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::setDomainSize(const std::array<T,Dim>& p_domainSize) { + m_domainSize = p_domainSize; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::setFftDomainConfiguration(const fft::FftDomainConfiguration<Dim>& p_fftConfig) { + m_inputFftConfig = p_fftConfig; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::includeOutputPeriodicBoundaries(bool p_val) { + m_includeOutputPeriodicBds = p_val; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::excludeOutputPeriodicBoundaries() { + return this->includeOutputPeriodicBoundaries(false); + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::setFftFlags(unsigned int p_flags) { + m_fftFlags = p_flags; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::appendFftFlags(unsigned int p_flags) { + m_fftFlags |= p_flags; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::enableOutputPeriodicBoundariesMirroring(bool p_val) { + m_mirrorOutputPeriodicBds = p_val; + return *this; + } + template <typename T, std::size_t Dim> + FftDiffSolver<T,Dim>& FftDiffSolver<T,Dim>::disableOutputPeriodicBoundariesMirroring() { + return this->enableOutputPeriodicBoundariesMirroring(false); + } + + /* Accessors */ + template <typename T, std::size_t Dim> + std::array<T,Dim> FftDiffSolver<T,Dim>::domainSize() const { + return m_domainSize; + } + template <typename T, std::size_t Dim> + fft::FftDomainConfiguration<Dim> FftDiffSolver<T,Dim>::inputFftConfig() const { + return m_inputFftConfig; + } + template <typename T, std::size_t Dim> + unsigned int FftDiffSolver<T,Dim>::fftFlags() const { + return m_fftFlags; + } + template <typename T, std::size_t Dim> + bool FftDiffSolver<T,Dim>::includeOutputPeriodicBds() const { + return m_includeOutputPeriodicBds; + } + template <typename T, std::size_t Dim> + bool FftDiffSolver<T,Dim>::mirrorOutputPeriodicBds() const { + return m_mirrorOutputPeriodicBds; + } + + /* Apply operator */ + template <typename T, std::size_t Dim> + void FftDiffSolver<T,Dim>::apply( + hysop::const_multi_array_ref<T,Dim> p_input, + hysop::multi_array_ref<T,Dim> p_output, + const std::array<int,Dim> &p_order) const { + + fft::Planner<T,Dim> planner; + bool success = planner.plan(p_input, p_output, m_inputFftConfig, p_order, m_domainSize, m_fftFlags, + m_includeOutputPeriodicBds, m_mirrorOutputPeriodicBds); + if(!success) + throw std::runtime_error("Failed to plan transforms !"); + + planner.executeForwardTransform(); + { + AccumulatorIndex<std::complex<T>,Dim> idx; + idx.setAccumulatorSource(planner.waveNumbersPows()).setAccumulatorFunction(std::multiplies<std::complex<T>>()); + + if(planner.transformType() == fft::FftTransformType::FFT_R2R) { + multi_array_view<T,Dim> view = planner.transformedRealData(); + //view.print("PRE-RDATA"); + + idx.reset(view.shape()); + while(!idx.atMaxId()) { + view(idx.ids()) *= idx.accumulatedVal().real()/planner.normalisationFactor(); + ++idx; + } + + //view.print("POST-RDATA"); + } + else if(planner.transformType() == fft::FftTransformType::FFT_R2C) { + multi_array_ref<std::complex<T>,Dim> ref = planner.transformedComplexData(); + std::complex<T> *data = ref.data(); + //ref.print("PRE-CDATA"); + + idx.reset(ref.shape()); + while(!idx.atMaxId()) { + data[idx()] *= idx.accumulatedVal()/planner.normalisationFactor(); + ++idx; + } + + //ref.print("POST-CDATA"); + } + } + planner.executeBackwardTransform(); + } + + } /* end of namespace solver */ +} /* end of namespace hysop */ + + + +#endif /* end of include guard: HYSOP_FFTDIFFSOLVER_H */ diff --git a/HySoP/src/hysop++/src/solver/fftPoissonSolver.h b/HySoP/src/hysop++/src/solver/fftPoissonSolver.h new file mode 100644 index 000000000..4662efb83 --- /dev/null +++ b/HySoP/src/hysop++/src/solver/fftPoissonSolver.h @@ -0,0 +1,205 @@ + +#ifndef HYSOP_FFTPOISSONSOLVER_H +#define HYSOP_FFTPOISSONSOLVER_H + +#include <cmath> + +#include "solver/poissonSolver.h" +#include "fft/planner.h" +#include "data/accumulatorIndex.h" + +namespace hysop { + namespace solver { + + template <typename T, std::size_t Dim> + class FftPoissonSolver : public PoissonSolver<T,Dim> { + + private: + using super = PoissonSolver<T,Dim>; + public: + FftPoissonSolver() = default; + FftPoissonSolver(const FftPoissonSolver& p_other) = default; + FftPoissonSolver( FftPoissonSolver&& p_other) = default; + ~FftPoissonSolver() = default; + + FftPoissonSolver& operator=(const FftPoissonSolver& p_other) = default; + FftPoissonSolver& operator=( FftPoissonSolver&& p_other) = default; + + FftPoissonSolver(const std::array<T,Dim> &p_domainSize, const domain::DomainConfiguration<Dim> &p_domainConfig, + unsigned int p_fftFlags = FFTW_MEASURE, + bool p_includeOutputPeriodicBds=false, bool p_mirrorOutputPeriodicBds=false); + + /* Mutators */ + FftPoissonSolver& setDomainSize(const std::array<T,Dim>& p_domainSize); + FftPoissonSolver& setInputDomainConfiguration(const domain::DomainConfiguration<Dim>& p_domainConfig); + + FftPoissonSolver& setFftFlags(unsigned int p_flags); + FftPoissonSolver& appendFftFlags(unsigned int p_flags); + + FftPoissonSolver& includeOutputPeriodicBoundaries(bool p_val = true); + FftPoissonSolver& excludeOutputPeriodicBoundaries(); + + FftPoissonSolver& enableOutputPeriodicBoundariesMirroring(bool p_val = true); + FftPoissonSolver& disableOutputPeriodicBoundariesMirroring(); + + /* Accessors */ + std::array<T,Dim> domainSize() const; + domain::DomainConfiguration<Dim> inputDomainConfig() const; + + unsigned int fftFlags() const; + bool includeOutputPeriodicBds() const; + bool mirrorOutputPeriodicBds() const; + + /* Apply operator */ + virtual void apply(hysop::const_multi_array_ref<T,Dim> p_input, + hysop::multi_array_ref<T,Dim> p_output) const final override; + + protected: + std::array<T,Dim> m_domainSize; + domain::DomainConfiguration<Dim> m_inputDomainConfig; + unsigned int m_fftFlags; + bool m_includeOutputPeriodicBds, m_mirrorOutputPeriodicBds; + }; + + + + /* Implementation */ + + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>::FftPoissonSolver(const std::array<T,Dim> &p_domainSize, const domain::DomainConfiguration<Dim> &p_inputDomainConfig, + unsigned int p_fftFlags, + bool p_includeOutputPeriodicBds, bool p_mirrorOutputPeriodicBds): + m_domainSize(p_domainSize), m_inputDomainConfig(p_inputDomainConfig), + m_fftFlags(p_fftFlags), + m_includeOutputPeriodicBds(p_includeOutputPeriodicBds), m_mirrorOutputPeriodicBds(p_mirrorOutputPeriodicBds) { + } + + + /* Mutators */ + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::setDomainSize(const std::array<T,Dim>& p_domainSize) { + m_domainSize = p_domainSize; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::setInputDomainConfiguration(const domain::DomainConfiguration<Dim>& p_inputDomainConfig) { + m_inputDomainConfig = p_inputDomainConfig; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::includeOutputPeriodicBoundaries(bool p_val) { + m_includeOutputPeriodicBds = p_val; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::excludeOutputPeriodicBoundaries() { + return this->includeOutputPeriodicBoundaries(false); + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::setFftFlags(unsigned int p_flags) { + m_fftFlags = p_flags; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::appendFftFlags(unsigned int p_flags) { + m_fftFlags |= p_flags; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::enableOutputPeriodicBoundariesMirroring(bool p_val) { + m_mirrorOutputPeriodicBds = p_val; + return *this; + } + template <typename T, std::size_t Dim> + FftPoissonSolver<T,Dim>& FftPoissonSolver<T,Dim>::disableOutputPeriodicBoundariesMirroring() { + return this->enableOutputPeriodicBoundariesMirroring(false); + } + + /* Accessors */ + template <typename T, std::size_t Dim> + std::array<T,Dim> FftPoissonSolver<T,Dim>::domainSize() const { + return m_domainSize; + } + template <typename T, std::size_t Dim> + domain::DomainConfiguration<Dim> FftPoissonSolver<T,Dim>::inputDomainConfig() const { + return m_inputDomainConfig; + } + template <typename T, std::size_t Dim> + unsigned int FftPoissonSolver<T,Dim>::fftFlags() const { + return m_fftFlags; + } + template <typename T, std::size_t Dim> + bool FftPoissonSolver<T,Dim>::includeOutputPeriodicBds() const { + return m_includeOutputPeriodicBds; + } + template <typename T, std::size_t Dim> + bool FftPoissonSolver<T,Dim>::mirrorOutputPeriodicBds() const { + return m_mirrorOutputPeriodicBds; + } + + /* Apply operator */ + template <typename T, std::size_t Dim> + void FftPoissonSolver<T,Dim>::apply( + hysop::const_multi_array_ref<T,Dim> p_input, + hysop::multi_array_ref<T,Dim> p_output) const { + + + fft::Planner<T,Dim> planner; + planner.setFixedAxeWaveNumberPows(std::complex<T>(0,0)); + { + std::array<int, Dim> order; + for (std::size_t d=0; d<Dim; d++) + order[d] = (m_inputDomainConfig[d].first == domain::Boundary::NONE ? 0 : 2); + bool success = planner.plan(p_input, p_output, m_inputDomainConfig, order, m_domainSize, m_fftFlags, + m_includeOutputPeriodicBds, m_mirrorOutputPeriodicBds); + if(!success) + throw std::runtime_error("Failed to plan transforms !"); + } + const T normalisationFactor = planner.normalisationFactor(); + //std::cout << planner << std::endl; + + planner.executeForwardTransform(); + { + AccumulatorIndex<std::complex<T>,Dim> idx; + idx.setAccumulatorSource(planner.waveNumbersPows()).setAccumulatorFunction(std::plus<std::complex<T>>()); + + if(planner.transformType() == fft::FftTransformType::FFT_R2R) { + multi_array_view<T,Dim> view = planner.transformedRealData(); + //view.print("PRE-RDATA"); + + idx.reset(view.shape()); + while(!idx.atMaxId()) { + T filter = idx.accumulatedVal().real(); + filter = (std::fpclassify(filter)==FP_ZERO ? T(0) : (T(1)/filter)*(T(1)/normalisationFactor)); + view(idx.ids()) *= filter; + ++idx; + } + + //view.print("POST-RDATA"); + } + else if(planner.transformType() == fft::FftTransformType::FFT_R2C) { + multi_array_ref<std::complex<T>,Dim> ref = planner.transformedComplexData(); + std::complex<T> *data = ref.data(); + //ref.print("PRE-CDATA"); + + idx.reset(ref.shape()); + while(!idx.atMaxId()) { + std::complex<T> filter = idx.accumulatedVal(); + filter = ((std::fpclassify(filter.real())==FP_ZERO) && (std::fpclassify(filter.imag())==FP_ZERO) ? + std::complex<T>(0,0) : (T(1)/filter)*(T(1)/normalisationFactor)); + data[idx()] *= filter; + ++idx; + } + + //ref.print("POST-CDATA"); + } + } + planner.executeBackwardTransform(); + } + + } /* end of namespace solver */ +} /* end of namespace hysop */ + + + +#endif /* end of include guard: HYSOP_FFTPOISSONSOLVER_H */ diff --git a/HySoP/src/hysop++/src/solver/poissonSolver.h b/HySoP/src/hysop++/src/solver/poissonSolver.h new file mode 100644 index 000000000..61b16440a --- /dev/null +++ b/HySoP/src/hysop++/src/solver/poissonSolver.h @@ -0,0 +1,27 @@ + +#ifndef HYSOP_POISSONSOLVER_H +#define HYSOP_POISSONSOLVER_H + +#include "data/multi_array/multi_array.h" + +namespace hysop { + namespace solver { + template <typename T, std::size_t Dim> + class PoissonSolver { + + public: + virtual void apply(hysop::const_multi_array_ref<T,Dim> input, + hysop::multi_array_ref<T,Dim> output) const = 0; + + void operator()(hysop::const_multi_array_ref<T,Dim> input, + hysop::multi_array_ref<T,Dim> output) { + this->apply(input,output); + } + + }; + + } /* end of namespace solver */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_POISSONSOLVER_H */ + diff --git a/HySoP/src/hysop++/src/utils/constants.h b/HySoP/src/hysop++/src/utils/constants.h new file mode 100644 index 000000000..e2bad74ba --- /dev/null +++ b/HySoP/src/hysop++/src/utils/constants.h @@ -0,0 +1,16 @@ + +#ifndef CONSTANTS_H +#define CONSTANTS_H + +#include <cmath> +#include "utils/types.h" + +namespace hysop { + namespace constants { + static constexpr hysop::types::complex I(0.0L,1.0L); + static constexpr hysop::types::complex Z(0.0L,0.0L); + static constexpr long double pi = acosl(-1.0L); + } +} + +#endif /* end of include guard: CONSTANTS_H */ diff --git a/HySoP/src/hysop++/src/utils/default.h b/HySoP/src/hysop++/src/utils/default.h new file mode 100644 index 000000000..efeacfb8d --- /dev/null +++ b/HySoP/src/hysop++/src/utils/default.h @@ -0,0 +1,20 @@ + +#ifndef HYSOP_DEFAULT_H +#define HYSOP_DEFAULT_H + +#include "data/memory/minimalAllocator.h" +#include "data/memory/fftwAllocator.h" + +namespace hysop { + namespace _default { + + template <typename T> + using allocator = hysop::data::memory::MinimalAllocator<T>; + + template <typename T> + using fft_allocator = hysop::data::memory::FftwAllocator<T>; + + } /* end of namespace _default */ +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_DEFAULT_H */ diff --git a/HySoP/src/hysop++/src/utils/defines.h b/HySoP/src/hysop++/src/utils/defines.h new file mode 100644 index 000000000..b456d75fa --- /dev/null +++ b/HySoP/src/hysop++/src/utils/defines.h @@ -0,0 +1,46 @@ + +#ifndef HYSOP_DEFINES_H +#define HYSOP_DEFINES_H + +#include <stdexcept> +#include <string> + +#if __cplusplus >= 201103L +#define HAS_CXX11 +#endif + +#if __cplusplus >= 201402L +#define HAS_CXX14 +#endif + +#define CAT(X,Y) CAT2(X,Y) +#define CAT2(X,Y) X##Y +#define CAT_2 CAT +#define CAT_3(X,Y,Z) CAT(X,CAT(Y,Z)) +#define CAT_4(A,X,Y,Z) CAT(A,CAT_3(X,Y,Z)) +#define CAT_5(A,B,X,Y,Z) CAT_3(A,B,CAT_3(X,Y,Z)) + +#define NOT_IMPLEMENTED_YET { \ + throw std::runtime_error("Function not implemented yet in " + std::string(__FILE__) + ":" + std::to_string(__LINE__) + "."); \ +} + +//Linux console colors +#define RESET "\033[0m" +#define BLACK "\033[30m" /* Black */ +#define RED "\033[31m" /* Red */ +#define GREEN "\033[32m" /* Green */ +#define YELLOW "\033[33m" /* Yellow */ +#define BLUE "\033[34m" /* Blue */ +#define MAGENTA "\033[35m" /* Magenta */ +#define CYAN "\033[36m" /* Cyan */ +#define WHITE "\033[37m" /* White */ +#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */ +#define BOLDRED "\033[1m\033[31m" /* Bold Red */ +#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */ +#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */ +#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */ +#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */ +#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */ +#define BOLDWHITE "\033[1m\033[37m" /* Bold White */ + +#endif /* end of include guard: HYSOP_DEFINES_H */ diff --git a/HySoP/src/hysop++/src/utils/types.h b/HySoP/src/hysop++/src/utils/types.h new file mode 100644 index 000000000..58dcdbc28 --- /dev/null +++ b/HySoP/src/hysop++/src/utils/types.h @@ -0,0 +1,59 @@ + +#ifndef HYSOP_TYPES_H +#define HYSOP_TYPES_H + +#include <complex> +#include <array> + +#include "utils/utils.h" +#include "utils/default.h" + +namespace hysop { + + /* forward declare external types */ + namespace data { + template <typename T, std::size_t Dim, typename Allocator> + class multi_array; + template <typename T, std::size_t Dim> + class multi_array_ref; + template <typename T, std::size_t Dim> + class multi_array_view; + template <typename T, std::size_t Dim> + class const_multi_array_view; + template <typename T, std::size_t Dim> + class const_multi_array_ref; + } /* end of namespace data */ + + + namespace types { + typedef double real; + typedef std::complex<real> complex; + } /* end of namespace types */ + + + /* expose those types to namespace hysop */ + template <std::size_t Dim> + using Shape = std::array<std::size_t, Dim>; + + template <std::size_t Dim> + using Offset = std::array<std::ptrdiff_t, Dim>; + + + template <typename T, std::size_t Dim, typename Allocator = hysop::_default::allocator<T>> + using multi_array = hysop::data::multi_array<T,Dim,Allocator>; + + template <typename T, std::size_t Dim> + using multi_array_view = hysop::data::multi_array_view<T,Dim>; + + template <typename T, std::size_t Dim> + using const_multi_array_view = hysop::data::const_multi_array_view<T,Dim>; + + template <typename T, std::size_t Dim> + using multi_array_ref = hysop::data::multi_array_ref<T,Dim>; + + template <typename T, std::size_t Dim> + using const_multi_array_ref = hysop::data::const_multi_array_ref<T,Dim>; + +} /* end of namespace hysop */ + +#endif /* end of include guard: HYSOP_TYPES_H */ diff --git a/HySoP/src/hysop++/src/utils/utils.cpp b/HySoP/src/hysop++/src/utils/utils.cpp new file mode 100644 index 000000000..501243272 --- /dev/null +++ b/HySoP/src/hysop++/src/utils/utils.cpp @@ -0,0 +1,17 @@ + +#include "utils/utils.h" + +namespace std { + std::ostream& operator<<(std::ostream& os, const fftw_iodim& iodim) { + os << "[n=" << iodim.n << ", is=" << iodim.is << ", os=" << iodim.os << "]"; + return os; + } +} + +namespace hysop { + namespace utils { + + } +} + + diff --git a/HySoP/src/hysop++/src/utils/utils.h b/HySoP/src/hysop++/src/utils/utils.h new file mode 100644 index 000000000..066516c28 --- /dev/null +++ b/HySoP/src/hysop++/src/utils/utils.h @@ -0,0 +1,196 @@ + +#ifndef HYSOP_UTILS_H +#define HYSOP_UTILS_H + +#include <fftw3.h> +#include <array> +#include <vector> +#include <tuple> +#include <iostream> +#include <limits> + +#include "detail/index_seq.h" +#include "boost/multi_array.hpp" + +namespace hysop { + namespace utils { + + template <typename... T> + void printTuple(std::ostream& os, const std::tuple<T...>& tuple); + + template <typename T> + bool areEqual(const T &lhs, const T &rhs); + template <typename T> + bool areNotEqual(const T &lhs, const T &rhs); + + + /* boost related utilities */ + template <std::size_t NumDims> + boost::detail::multi_array::index_gen<NumDims, NumDims> buildView(); + + template <std::size_t NumDims> + boost::detail::multi_array::index_gen<NumDims, NumDims> buildIndices( + const std::array<boost::multi_array_types::index_range, NumDims> &p_ranges); + + template <std::size_t NumRanges> + boost::detail::multi_array::extent_gen<NumRanges> buildExtents( + const std::array<std::size_t, NumRanges> &p_shape); + + + /* Implementation */ + + template <typename Tuple, int... I> + void printTupleImpl(std::ostream& os, const Tuple& tuple, hysop::detail::index_seq<I...>) { + const int dummy[sizeof...(I)] = { (os << std::get<I>(tuple) << ",", 0)... }; + os << std::get<sizeof...(I)>(tuple); + } + template <typename... T> + void printTuple(std::ostream& os, const std::tuple<T...>& tuple) { + os << "("; + printTupleImpl(os,tuple, hysop::detail::index_seq_gen<sizeof...(T)-1>()); + os << ")"; + } + + template <typename T, typename> + bool areEqualImpl(const T& lhs, const T& rhs) { + return lhs == rhs; + } + template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type* = nullptr> + bool areEqualImpl(const T& lhs, const T& rhs) { + return (std::abs(rhs - lhs) <= std::numeric_limits<T>::epsilon() * std::max(std::abs(lhs), std::abs(rhs))); + } + + template <typename T> + bool areEqual(const T &lhs, const T &rhs) { + return areEqualImpl<T>(lhs,rhs); + } + template <typename T> + bool areNotEqual(const T &lhs, const T &rhs) { + return !areEqualImpl<T>(lhs,rhs); + } + + + /* boost related utilities */ + template <std::size_t NumDims> + struct BuildViewImpl { + static_assert(NumDims >= 1, "NumDims cannot be < 1"); + boost::detail::multi_array::index_gen<NumDims, NumDims> build() const { + return BuildViewImpl<NumDims-1>().build()[boost::multi_array_types::index_range()]; + } + }; + template <> + struct BuildViewImpl<1> { + boost::detail::multi_array::index_gen<1,1> build() const { + return boost::multi_array_types::index_gen()[boost::multi_array_types::index_range()]; + } + }; + template <std::size_t NumDims> + boost::detail::multi_array::index_gen<NumDims, NumDims> buildView() { + return BuildViewImpl<NumDims>().build(); + } + + template <std::size_t NumDims, std::size_t K=NumDims> + struct BuildIndicesImpl { + static_assert(NumDims >= 1, "NumDims cannot be < 1"); + const std::array<boost::multi_array_types::index_range, NumDims> &m_ranges; + BuildIndicesImpl(const std::array<boost::multi_array_types::index_range, NumDims> &p_ranges): m_ranges(p_ranges) {} + boost::detail::multi_array::index_gen<K,K> build() const { + return BuildIndicesImpl<NumDims,K-1>(m_ranges).build()[m_ranges[K-1]]; + } + }; + template <std::size_t NumDims> + struct BuildIndicesImpl<NumDims,1> { + const std::array<boost::multi_array_types::index_range, NumDims> &m_ranges; + BuildIndicesImpl(const std::array<boost::multi_array_types::index_range, NumDims> &p_ranges): m_ranges(p_ranges) {} + boost::detail::multi_array::index_gen<1,1> build() const { + return boost::multi_array_types::index_gen()[m_ranges[0]]; + } + }; + template <std::size_t NumDims> + boost::detail::multi_array::index_gen<NumDims, NumDims> buildIndices( + const std::array<boost::multi_array_types::index_range, NumDims> &p_ranges) { + return BuildIndicesImpl<NumDims>(p_ranges).build(); + } + + template <std::size_t NumRanges, std::size_t K=NumRanges> + struct BuildExtentImpl { + static_assert(NumRanges >= 1, "NumDims cannot be < 1"); + const std::array<std::size_t,NumRanges>& m_shape; + BuildExtentImpl(const std::array<std::size_t, NumRanges>& p_shape): m_shape(p_shape) {} + boost::detail::multi_array::extent_gen<K> build() const { + return BuildExtentImpl<NumRanges,K-1>(m_shape).build()[m_shape[K-1]]; + } + }; + template <std::size_t NumRanges> + struct BuildExtentImpl<NumRanges,1> { + const std::array<std::size_t,NumRanges>& m_shape; + BuildExtentImpl(const std::array<std::size_t, NumRanges>& p_shape): m_shape(p_shape) {} + boost::detail::multi_array::extent_gen<1> build() const { + return boost::multi_array_types::extent_gen()[m_shape[0]]; + } + }; + template <std::size_t NumRanges> + boost::detail::multi_array::extent_gen<NumRanges> buildExtents( + const std::array<std::size_t, NumRanges> &p_shape) { + return BuildExtentImpl<NumRanges>(p_shape).build(); + } + } +} + + +/* quick and dirty fix to allow non namespace dependant operators << for std containers */ +namespace std { + + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const std::array<T,Dim>& array); + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const boost::array<T,Dim>& array); + + template <typename T> + std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector); + + template <typename... T> + std::ostream& operator<<(std::ostream& os, const std::tuple<T...>& tuple); + + std::ostream& operator<<(std::ostream& os, const fftw_iodim& iodim); + + + + /* Implementation */ + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const std::array<T,Dim>& array) { + os << "["; + for (std::size_t i = 0; i < Dim-1; i++) + os << array[i] << ","; + os << array[Dim-1]; + os << "]"; + return os; + } + template <typename T, std::size_t Dim> + std::ostream& operator<<(std::ostream& os, const boost::array<T,Dim>& array) { + os << "["; + for (std::size_t i = 0; i < Dim-1; i++) + os << array[i] << ","; + os << array[Dim-1]; + os << "]"; + return os; + } + template <typename T> + std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector) { + os << "["; + if(!vector.empty()) { + for (std::size_t i = 0; i < vector.size()-1; i++) + os << vector[i] << ","; + os << vector[vector.size()-1]; + } + os << "]"; + return os; + } + template <typename... T> + std::ostream& operator<<(std::ostream& os, const std::tuple<T...>& tuple) { + hysop::utils::printTuple(os,tuple); + return os; + } +} + +#endif /* end of include guard: HYSOP_UTILS_H */ diff --git a/HySoP/src/hysop++/tests/CMakeLists.txt b/HySoP/src/hysop++/tests/CMakeLists.txt new file mode 100644 index 000000000..e16d2d63c --- /dev/null +++ b/HySoP/src/hysop++/tests/CMakeLists.txt @@ -0,0 +1,7 @@ + +include("${CMAKE_SOURCE_DIR}/CMake/GoogleTestHelper.cmake") + +add_subdirectory("testPolynoms") +add_subdirectory("testPlanner") +add_subdirectory("testDiffSolver") +add_subdirectory("testPoissonSolver") diff --git a/HySoP/src/hysop++/tests/testDiffSolver/CMakeLists.txt b/HySoP/src/hysop++/tests/testDiffSolver/CMakeLists.txt new file mode 100644 index 000000000..eda82de15 --- /dev/null +++ b/HySoP/src/hysop++/tests/testDiffSolver/CMakeLists.txt @@ -0,0 +1,12 @@ + +file(GLOB CPP_SRCS *.cpp) +set(SRCS ${CPP_SRCS}) + +get_filename_component(test_name ${CMAKE_CURRENT_SOURCE_DIR} NAME) +add_executable(${test_name} ${SRCS}) + +target_link_libraries(${test_name} ${HYSOP_CXX_LIBRARY}) +target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${CXX_EXT_LIBS}) + +add_test("${test_name}" "${test_name}") + diff --git a/HySoP/src/hysop++/tests/testDiffSolver/main.cpp b/HySoP/src/hysop++/tests/testDiffSolver/main.cpp new file mode 100644 index 000000000..b6141085e --- /dev/null +++ b/HySoP/src/hysop++/tests/testDiffSolver/main.cpp @@ -0,0 +1,10 @@ + +#include "gtest/gtest.h" + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + int ret = RUN_ALL_TESTS(); + + return ret; +} diff --git a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp new file mode 100644 index 000000000..7e66b594c --- /dev/null +++ b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp @@ -0,0 +1,221 @@ + +#include "testDiffSolver.h" + +#include <cstdlib> + +#include "domain/domain.h" +#include "solver/fftDiffSolver.h" +#include "data/multi_array/multi_array.h" +#include "utils/constants.h" +#include "fft/extension.h" + +using namespace hysop; +using namespace hysop::domain; + +static constexpr std::size_t nExtensions = 4 ; +static constexpr std::size_t nExtensionsPair = 7 ; +static constexpr fft::Extension ext[nExtensions] = +{ fft::Extension::NONE, fft::Extension::ODD, fft::Extension::EVEN, fft::Extension::PERIODIC }; +static constexpr std::pair<fft::Extension,fft::Extension> pext[nExtensionsPair] { + std::make_pair(ext[3],ext[3]), //periodic-periodic + std::make_pair(ext[3],ext[3]), //periodic-periodic + std::make_pair(ext[2],ext[1]), //even-odd + std::make_pair(ext[1],ext[2]), //odd-even + std::make_pair(ext[2],ext[2]), //even-even + std::make_pair(ext[1],ext[1]), //odd-odd + std::make_pair(ext[0],ext[0]), //none-none +}; + +static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L }; + +template <typename T> +std::function<T(T)> func(std::size_t k) { + switch(k) { + case 0: return [=](T x) {return cos(T(freqs[0])*x);}; + case 1: return [=](T x) {return sin(T(freqs[1])*x);}; + case 2: return [=](T x) {return cos(T(freqs[2])*x);}; + case 3: return [=](T x) {return sin(T(freqs[3])*x);}; + case 4: return [=](T x) {return cos(T(freqs[4])*x);}; + case 5: return [=](T x) {return sin(T(freqs[5])*x);}; + default: return[=](T x) { return T(1); }; + } +} + +template <typename T> +std::function<T(T)> derivative(std::size_t k, int order) { + bool even = (k%2==0); + std::size_t p, offset; + T sign, coeff; + if(k>5) { + if(order != 0) + throw std::runtime_error("Non zero order !"); + return func<T>(k); + } + else if(even) { /* cos func */ + p = (order%2==0 ? k : k+1); + sign = std::pow(T(-1),(order+1)/2); + coeff = std::pow(freqs[k], order); + } + else { /* sin func */ + p = (order%2==0 ? k : k-1); + sign = std::pow(T(-1),order/2); + coeff = std::pow(freqs[k], order); + } + return [=](T x) { return sign*coeff*(func<T>(p)(x)); }; +} + +template <typename T, std::size_t Dim, bool verbose=false> +void test(std::size_t p_maxOrder, bool includePeriodicBds=false) { + Shape<Dim> shape; + typename Domain<T,Dim>::DomainSize domainSize; + Domain<T,Dim> ref, inBuffer, outBuffer; + + Domain<T,Dim>& in = inBuffer; + Domain<T,Dim>& out = outBuffer; + + std::array<int,Dim> order; + + shape.fill(8); + domainSize.fill(2*hysop::constants::pi); + + T eps = std::numeric_limits<T>::epsilon(); + const std::size_t N = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<std::size_t>()); + + ref.resize(domainSize).reshape(shape); + in = ref; + out = ref; + + Shape<Dim> maxOrder, testCases; + maxOrder.fill(p_maxOrder+1); + testCases.fill(nExtensionsPair); + Index<Dim> orderId(maxOrder); + Index<Dim> testCaseId; + std::size_t testCaseCount; + while(!(++orderId).atMaxId()) { + std::cout << " ::Order::" << orderId.ids() << (verbose ? "\n" : ""); + + std::array<T,3> meanDists; + meanDists.fill(0); + testCaseId.reset(testCases); + testCaseCount = testCaseId.maxId(); + while(!testCaseId.atMaxId()) { + std::copy(orderId.ids().begin(),orderId.ids().end(), order.begin()); + + /* generate transform configuration */ + std::array<std::pair<fft::Extension,fft::Extension>, Dim> extConfig; + for (std::size_t k=0; k<Dim; k++) { + std::size_t id = testCaseId[k]; + extConfig[k] = pext[id]; + if(pext[id].first==fft::Extension::NONE) + order[k] = 0; + } + fft::FftDomainConfiguration<Dim> domainConfig(extConfig, includePeriodicBds); + + const std::size_t orderSum = std::accumulate(order.begin(), order.end(), 0); + if(orderSum == 0) { + testCaseCount--; + ++testCaseId; + continue; + } + T orderPow = std::pow(T(10),T(orderSum)); + if(std::is_same<T,long double>::value) /* just in case long doubles are not hardware supported... */ + orderPow *= 1e3; + const auto criteria = std::make_tuple(orderPow*eps*N,orderPow*eps*sqrt(N),2*orderPow*eps); + + const auto f = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = func<T>(testCaseId[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= func<T>(testCaseId[d])(x[d]); + return val; + }; + const auto d = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = derivative<T>(testCaseId[0],order[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= derivative<T>(testCaseId[d],order[d])(x[d]); + return val; + }; + { + ref.resetDomainConfiguration(domainConfig.boundariesConfiguration()); + in = ref; + out = ref; + + in.apply(f); + ref.apply(d); + out.data().apply([](T& v){ v=T(0);}); + } + + solver::FftDiffSolver<T,Dim> solver(domainSize, domainConfig, FFTW_MEASURE, includePeriodicBds, includePeriodicBds); + solver.apply(in.data(), out.data(), order); + + std::stringstream ss; + ss << "["; + for (std::size_t k=0; k<Dim-1; k++) + ss << extConfig[k].first << "/" << extConfig[k].second << ","; + ss << extConfig[Dim-1].first << "/" << extConfig[Dim-1].second; + ss << "]"; + + const auto dist = out.distance(ref); + const bool pass = (std::get<0>(dist) < std::get<0>(criteria)) + && (std::get<1>(dist) < std::get<1>(criteria)) + && (std::get<2>(dist) < std::get<2>(criteria)); + + if((pass && verbose) || !pass) { + std::cout << (pass ? GREEN : RED); + std::cout << "\t" << std::setw(Dim*15) << ss.str() << " => " << (pass ? "OK" : "KO") + << " " << RESET << std::scientific << std::setprecision(2) << dist << std::endl; + } + if(!pass) { + //in.print("IN"); + //ref.print("REF"); + //out.print("OUT"); + std::cout << "Test failed => Criteria was: " << criteria << std::endl; + } + + meanDists[0] += std::get<0>(dist); + meanDists[1] += std::get<1>(dist); + meanDists[2] += std::get<2>(dist); + + EXPECT_TRUE(pass); + + ++testCaseId; + } + for (std::size_t k = 0; k < 3; k++) + meanDists[k] /= T(testCaseCount); + std::cout << "=> mean distances over " << std::scientific << std::setprecision(1) << std::setw(4) + << testCaseCount << " testcases: " << meanDists; + for (std::size_t k = 0; k < 3; k++) + meanDists[k] = std::round(meanDists[k]/eps); + std::cout << " ~= " << std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; + } +} + +TEST_F(DiffSolverTest, FloatDerivatives) { + std::cout << std::endl; + std::cout << "== TEST 1D - float ==" << std::endl; + test<float,1,false>(5); + std::cout << "== TEST 2D - float ==" << std::endl; + test<float,2,false>(3); + std::cout << "== TEST 3D - float ==" << std::endl; + test<float,3,false>(1); +} + +TEST_F(DiffSolverTest, DoubleDerivatives) { + std::cout << std::endl; + std::cout << "== TEST 1D - double ==" << std::endl; + test<double,1,false>(5); + std::cout << "== TEST 2D - double ==" << std::endl; + test<double,2,false>(3); + std::cout << "== TEST 3D - double ==" << std::endl; + test<double,3,false>(1); +} + +TEST_F(DiffSolverTest, LongDoubleDerivatives) { + std::cout << std::endl; + std::cout << "== TEST 1D - long double ==" << std::endl; + test<long double,1,false>(5); + std::cout << "== TEST 2D - long double ==" << std::endl; + test<long double,2,false>(3); + std::cout << "== TEST 3D - long double ==" << std::endl; + test<long double,3,false>(1); +} + diff --git a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.h b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.h new file mode 100644 index 000000000..b51090915 --- /dev/null +++ b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.h @@ -0,0 +1,14 @@ + +#include "gtest/gtest.h" +#include "solver/fftDiffSolver.h" +#include "fft/extension.h" + +using namespace hysop; + +class DiffSolverTest : public ::testing::Test { + protected: + DiffSolverTest() {} + void SetUp() {} + void TearDown() {} + virtual ~DiffSolverTest() {} +}; diff --git a/HySoP/src/hysop++/tests/testPlanner/CMakeLists.txt b/HySoP/src/hysop++/tests/testPlanner/CMakeLists.txt new file mode 100644 index 000000000..eda82de15 --- /dev/null +++ b/HySoP/src/hysop++/tests/testPlanner/CMakeLists.txt @@ -0,0 +1,12 @@ + +file(GLOB CPP_SRCS *.cpp) +set(SRCS ${CPP_SRCS}) + +get_filename_component(test_name ${CMAKE_CURRENT_SOURCE_DIR} NAME) +add_executable(${test_name} ${SRCS}) + +target_link_libraries(${test_name} ${HYSOP_CXX_LIBRARY}) +target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${CXX_EXT_LIBS}) + +add_test("${test_name}" "${test_name}") + diff --git a/HySoP/src/hysop++/tests/testPlanner/main.cpp b/HySoP/src/hysop++/tests/testPlanner/main.cpp new file mode 100644 index 000000000..b6141085e --- /dev/null +++ b/HySoP/src/hysop++/tests/testPlanner/main.cpp @@ -0,0 +1,10 @@ + +#include "gtest/gtest.h" + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + int ret = RUN_ALL_TESTS(); + + return ret; +} diff --git a/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp b/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp new file mode 100644 index 000000000..539d03426 --- /dev/null +++ b/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp @@ -0,0 +1,180 @@ + +#include "testPlanner.h" + +#include "data/multi_array/multi_array.h" +#include "domain/domain.h" +#include "utils/constants.h" +#include "fft/planner.h" +#include "fft/extension.h" + + +using namespace hysop; +using namespace hysop::domain; + +static constexpr std::size_t nExtensions = 4 ; +static constexpr std::size_t nExtensionsPair = 6 ; +static constexpr fft::Extension ext[nExtensions] = +{ fft::Extension::NONE, fft::Extension::ODD, fft::Extension::EVEN, fft::Extension::PERIODIC }; +static constexpr std::pair<fft::Extension,fft::Extension> pext[nExtensionsPair] { + std::make_pair(ext[0],ext[0]), //none-none + std::make_pair(ext[1],ext[1]), //odd-odd + std::make_pair(ext[1],ext[2]), //odd-even + std::make_pair(ext[1],ext[2]), //even-odd + std::make_pair(ext[2],ext[2]), //even-even + std::make_pair(ext[3],ext[3]), //periodic-periodic +}; + +template <typename T, std::size_t Dim, bool verbose=false> +void test(bool inplace, bool includePeriodicBds); + +TEST_F(PlannerTest, InplaceFloatTransforms) { + std::cout << std::endl; + std::cout << "== TEST 1D - float ==\t"; + test<float,1>(false,true); + std::cout << "== TEST 2D - float ==\t"; + test<float,2>(false,true); + std::cout << "== TEST 3D - float ==\t"; + test<float,3>(false,true); +} + +TEST_F(PlannerTest, InplaceDoubleTransforms) { + std::cout << std::endl; + std::cout << "== TEST 1D - double ==\t"; + test<double,1>(true,true); + std::cout << "== TEST 2D - double ==\t"; + test<double,2>(true,true); + std::cout << "== TEST 3D - double ==\t"; + test<double,3>(true,true); +} + +TEST_F(PlannerTest, InplaceLongDoubleTransforms) { + std::cout << std::endl; + std::cout << "== TEST 1D - long double ==\t"; + test<long double,1>(false,false); + std::cout << "== TEST 2D - long double ==\t"; + test<long double,2>(false,false); + std::cout << "== TEST 3D - long double ==\t"; + test<long double,3>(false,false); +} + +template <typename T, std::size_t Dim, bool verbose> +void test(bool inplace, bool includePeriodicBds) { + Shape<Dim> shape; + typename Domain<T,Dim>::DomainSize domainSize; + Domain<T,Dim> ref, inBuffer, outBuffer; + + Domain<T,Dim>& in = inBuffer; + Domain<T,Dim>& out = (inplace ? inBuffer : outBuffer); + + fft::Planner<T,Dim> planner; + std::array<int,Dim> order; + + const std::size_t nPoints = 16; + shape.fill(nPoints); + domainSize.fill(1.0); + order.fill(2); + + const T eps = std::numeric_limits<T>::epsilon(); + const std::size_t N = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<std::size_t>()); + const auto criteria = std::make_tuple(50*eps*N,sqrt(50)*eps*N,700*eps); + + ref.resize(domainSize).reshape(shape); + in = ref; + out = ref; + + Shape<Dim> testCases; + testCases.fill(nExtensionsPair); + Index<Dim> testCaseId(testCases); + std::array<T,3> meanDists{0}; + while(!testCaseId.atMaxId()) { + /* generate transform configuration */ + std::array<std::pair<fft::Extension,fft::Extension>, Dim> extConfig; + for (std::size_t k=0; k<Dim; k++) { + std::size_t id = testCaseId[k]; + extConfig[k] = pext[id]; + } + fft::FftDomainConfiguration<Dim> domainConfig(extConfig, includePeriodicBds); + + const auto f = [&](T &val, const hysop::Index<Dim>& idx) { + val = static_cast<T>(rand())/static_cast<T>(RAND_MAX); + for (std::size_t d=0; d<Dim; d++) { + if(idx[d]==0) { + if(extConfig[d].first == fft::Extension::ODD) { + val=T(0); + return; + } + else if(extConfig[d].first == fft::Extension::PERIODIC) + val=T(0.42); + } + else if(std::size_t(idx[d]) == idx.dim()[d]-1) { + if(extConfig[d].second == fft::Extension::ODD) { + val=T(0); + return; + } + else if(extConfig[d].second == fft::Extension::PERIODIC && includePeriodicBds) + val=T(0.42); + } + } + }; + + if(includePeriodicBds) + ref.resetDomainConfiguration(domainConfig.boundariesConfiguration()); + + /* fill reference and copy into input buffer */ + ref.data().apply(f); + in = ref; + + /* plan transforms and check if planning succeeded */ + assert(planner.plan(in.data(), out.data(), domainConfig, order, domainSize, FFTW_MEASURE, + includePeriodicBds, includePeriodicBds) || testCaseId()==0); + + /* execute forward and backward inplace transforms */ + planner.executeForwardTransform(); + { + if(planner.transformType() == fft::FftTransformType::FFT_R2C) + planner.transformedComplexData().apply([&](std::complex<T>& val) { val /= planner.normalisationFactor(); }); + else if(planner.transformType() == fft::FftTransformType::FFT_R2R) + planner.transformedRealData().apply([&](T& val) { val /= planner.normalisationFactor(); }); + } + planner.executeBackwardTransform(); + + std::stringstream ss; + ss << "["; + for (std::size_t k=0; k<Dim-1; k++) + ss << extConfig[k].first << "/" << extConfig[k].second << ","; + ss << extConfig[Dim-1].first << "/" << extConfig[Dim-1].second; + ss << "]"; + + const auto dist = out.distance(ref); + const bool pass = (std::get<0>(dist) < std::get<0>(criteria)) + && (std::get<1>(dist) < std::get<1>(criteria)) + && (std::get<2>(dist) < std::get<2>(criteria)); + + if((pass && verbose) || !pass) { + std::cout << (pass ? GREEN : RED); + std::cout << "\t" << std::setw(Dim*15) << ss.str() << " => " << (pass ? "OK" : "KO") + << " " << RESET << std::scientific << std::setprecision(2) << dist << std::endl; + } + if(!pass) { + if(!inplace) + in.print("IN"); + ref.print("REF"); + out.print("OUT"); + std::cout << planner << std::endl; + exit(EXIT_FAILURE); + } + + meanDists[0] += std::get<0>(dist); + meanDists[1] += std::get<1>(dist); + meanDists[2] += std::get<2>(dist); + + ++testCaseId; + } + for (std::size_t k = 0; k < 3; k++) + meanDists[k] /= T(testCaseId.maxId()); + std::cout << "Mean distances over " << std::scientific << std::setprecision(1) << std::setw(4) + << testCaseId.maxId() << " testcases: " << meanDists; + for (std::size_t k = 0; k < 3; k++) + meanDists[k] = std::round(meanDists[k]/eps); + std::cout << " ~= " << std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; +} diff --git a/HySoP/src/hysop++/tests/testPlanner/testPlanner.h b/HySoP/src/hysop++/tests/testPlanner/testPlanner.h new file mode 100644 index 000000000..5be17aa72 --- /dev/null +++ b/HySoP/src/hysop++/tests/testPlanner/testPlanner.h @@ -0,0 +1,20 @@ + +#include "gtest/gtest.h" +#include "fft/planner.h" +#include "fft/extension.h" + +using T = double; +constexpr std::size_t Dim = 1; + +using namespace hysop; + +class PlannerTest : public ::testing::Test { + protected: + PlannerTest() {} + void SetUp() {} + void TearDown() {} + virtual ~PlannerTest() {} + + public: + fft::Planner<T,Dim> planner; +}; diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/CMakeLists.txt b/HySoP/src/hysop++/tests/testPoissonSolver/CMakeLists.txt new file mode 100644 index 000000000..eda82de15 --- /dev/null +++ b/HySoP/src/hysop++/tests/testPoissonSolver/CMakeLists.txt @@ -0,0 +1,12 @@ + +file(GLOB CPP_SRCS *.cpp) +set(SRCS ${CPP_SRCS}) + +get_filename_component(test_name ${CMAKE_CURRENT_SOURCE_DIR} NAME) +add_executable(${test_name} ${SRCS}) + +target_link_libraries(${test_name} ${HYSOP_CXX_LIBRARY}) +target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${CXX_EXT_LIBS}) + +add_test("${test_name}" "${test_name}") + diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/main.cpp b/HySoP/src/hysop++/tests/testPoissonSolver/main.cpp new file mode 100644 index 000000000..b6141085e --- /dev/null +++ b/HySoP/src/hysop++/tests/testPoissonSolver/main.cpp @@ -0,0 +1,10 @@ + +#include "gtest/gtest.h" + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + int ret = RUN_ALL_TESTS(); + + return ret; +} diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp new file mode 100644 index 000000000..e144e167f --- /dev/null +++ b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp @@ -0,0 +1,186 @@ + +#include "testPoissonSolver.h" + +#include <cstdlib> + +#include "domain/domain.h" +#include "solver/fftPoissonSolver.h" +#include "data/multi_array/multi_array.h" +#include "utils/constants.h" +#include "domain/boundary.h" + +using namespace hysop; +using namespace hysop::domain; + +static constexpr std::size_t nBoundaries = 4; +static constexpr std::size_t nBoundaryPairs = 7; +static constexpr domain::Boundary bds[nBoundaries] = +{ domain::Boundary::NONE, domain::Boundary::HOMOGENEOUS_NEUMANN, domain::Boundary::HOMOGENEOUS_DIRICHLET, domain::Boundary::PERIODIC }; +static constexpr std::pair<domain::Boundary,domain::Boundary> pbds[nBoundaryPairs] { + std::make_pair(bds[3],bds[3]), //periodic-periodic + std::make_pair(bds[3],bds[3]), //periodic-periodic + std::make_pair(bds[2],bds[1]), //even-odd + std::make_pair(bds[1],bds[2]), //odd-even + std::make_pair(bds[2],bds[2]), //even-even + std::make_pair(bds[1],bds[1]), //odd-odd + std::make_pair(bds[0],bds[0]), //none-none +}; + +static constexpr long double freqs[6] = { 1.0L, 1.0L, 0.75L, 0.75L, 0.50L, 0.50L }; + +template <typename T> +std::function<T(T)> func(std::size_t k) { + switch(k) { + case 0: return [=](T x) {return cos(T(freqs[0])*x);}; + case 1: return [=](T x) {return sin(T(freqs[1])*x);}; + case 2: return [=](T x) {return cos(T(freqs[2])*x);}; + case 3: return [=](T x) {return sin(T(freqs[3])*x);}; + case 4: return [=](T x) {return cos(T(freqs[4])*x);}; + case 5: return [=](T x) {return sin(T(freqs[5])*x);}; + default: return[=](T x) { return T(1); }; + } +} + +template <typename T, std::size_t Dim, bool verbose=false> +void test(bool includePeriodicBds=false) { + Shape<Dim> shape; + typename Domain<T,Dim>::DomainSize domainSize; + Domain<T,Dim> ref, inBuffer, outBuffer; + + Domain<T,Dim>& in = inBuffer; + Domain<T,Dim>& out = outBuffer; + + shape.fill(16); + domainSize.fill(2*hysop::constants::pi); + + T eps = std::numeric_limits<T>::epsilon(); + const std::size_t N = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<std::size_t>()); + + ref.resize(domainSize).reshape(shape); + in = ref; + out = ref; + + Shape<Dim> testCases; + testCases.fill(nBoundaryPairs); + Index<Dim> testCaseId(testCases); + std::array<T,3> meanDists{0}; + std::size_t testCaseCount = testCaseId.maxId()-1; + + if(verbose) + std::cout << std::endl; + + while(testCaseId() != testCaseId.maxId()-1) { + + /* generate transform configuration */ + std::size_t orderSum = 0; + std::array<std::pair<domain::Boundary,domain::Boundary>, Dim> bdsConfig; + T W2sum = T(0); + for (std::size_t k=0; k<Dim; k++) { + std::size_t id = testCaseId[k]; + bdsConfig[k] = pbds[id]; + if(bdsConfig[k].first != domain::Boundary::NONE) { + W2sum += freqs[id]*freqs[id]; + orderSum+=2; + } + } + domain::DomainConfiguration<Dim> domainConfig(bdsConfig, includePeriodicBds); + + T orderPow = std::pow(T(10),T(orderSum)); + if(std::is_same<T,long double>::value) /* just in case long doubles are not hardware supported... */ + orderPow *= 1e3; + const auto criteria = std::make_tuple(orderPow*eps*N,orderPow*eps*sqrt(N),2*orderPow*eps); + + const auto phi = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + T val = func<T>(testCaseId[0])(x[0]); + for (std::size_t d=1; d < Dim; d++) + val *= func<T>(testCaseId[d])(x[d]); + return val; + }; + const auto f = [&](const typename Domain<T,Dim>::SpaceVariable &x) { + return -W2sum*phi(x); + }; + + { + ref.resetDomainConfiguration(domainConfig); + in = ref; + out = ref; + + in.apply(f); + ref.apply(phi); + out.data().apply([](T& v){ v=T(0);}); + } + + solver::FftPoissonSolver<T,Dim> solver(domainSize, domainConfig, FFTW_MEASURE, includePeriodicBds, includePeriodicBds); + solver.apply(in.data(), out.data()); + + std::stringstream ss; + ss << "["; + for (std::size_t k=0; k<Dim-1; k++) + ss << bdsConfig[k].first << "/" << bdsConfig[k].second << ","; + ss << bdsConfig[Dim-1].first << "/" << bdsConfig[Dim-1].second; + ss << "]"; + + const auto dist = out.distance(ref); + const bool pass = (std::get<0>(dist) < std::get<0>(criteria)) + && (std::get<1>(dist) < std::get<1>(criteria)) + && (std::get<2>(dist) < std::get<2>(criteria)); + + if((pass && verbose) || !pass) { + std::cout << (pass ? GREEN : RED); + std::cout << "\t" << std::setw(Dim*15) << ss.str() << " => " << (pass ? "OK" : "KO") + << " " << RESET << std::scientific << std::setprecision(2) << dist << std::endl; + } + if(!pass) { + //in.print("IN"); + //ref.print("REF"); + //out.print("OUT"); + std::cout << "\t\tTest Failed... Criteria was " << criteria << "." << std::endl; + } + EXPECT_TRUE(pass); + + meanDists[0] += std::get<0>(dist); + meanDists[1] += std::get<1>(dist); + meanDists[2] += std::get<2>(dist); + + ++testCaseId; + } + + for (std::size_t k = 0; k < 3; k++) + meanDists[k] /= T(testCaseCount); + std::cout << "\t=> mean distances over " << std::scientific << std::setprecision(1) << std::setw(4) + << testCaseCount << " testcases: " << meanDists; + for (std::size_t k = 0; k < 3; k++) + meanDists[k] = std::round(meanDists[k]/eps); + std::cout << " ~= " << std::fixed << std::setprecision(0) << meanDists << " eps" << std::endl; +} + + +TEST_F(PoissonSolverTest, FloatPoissonSolver) { + std::cout << std::endl; + std::cout << "== TEST 1D - float =="; + test<float,1,false>(true); + std::cout << "== TEST 2D - float =="; + test<float,2,false>(); + std::cout << "== TEST 3D - float =="; + test<float,3,false>(); +} + +TEST_F(PoissonSolverTest, DoublePoissonSolver) { + std::cout << std::endl; + std::cout << "== TEST 1D - double =="; + test<double,1,false>(); + std::cout << "== TEST 2D - double =="; + test<double,2,false>(); + std::cout << "== TEST 3D - double =="; + test<double,3,false>(); +} + +TEST_F(PoissonSolverTest, LongDoublePoissonSolver) { + std::cout << std::endl; + std::cout << "== TEST 1D - long double =="; + test<long double,1,false>(); + std::cout << "== TEST 2D - long double =="; + test<long double,2,false>(); + std::cout << "== TEST 3D - long double =="; + test<long double,3,false>(); +} diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.h b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.h new file mode 100644 index 000000000..a5115abef --- /dev/null +++ b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.h @@ -0,0 +1,14 @@ + +#include "gtest/gtest.h" +#include "solver/fftPoissonSolver.h" +#include "fft/extension.h" + +using namespace hysop; + +class PoissonSolverTest : public ::testing::Test { + protected: + PoissonSolverTest() {} + void SetUp() {} + void TearDown() {} + virtual ~PoissonSolverTest() {} +}; diff --git a/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt b/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt new file mode 100644 index 000000000..802cfd2e2 --- /dev/null +++ b/HySoP/src/hysop++/tests/testPolynoms/CMakeLists.txt @@ -0,0 +1,12 @@ + +file(GLOB CPP_SRCS *.cpp) +set(SRCS ${CPP_SRCS}) + +get_filename_component(test_name ${CMAKE_CURRENT_SOURCE_DIR} NAME) +add_executable(${test_name} ${SRCS}) + +target_link_libraries(${test_name} ${HYSOP_LIBRARY}) +target_link_libraries(${test_name} ${GTEST_LIBRARIES} ${EXT_LIBRARIES}) + +add_test("${test_name}" "${test_name}") + diff --git a/HySoP/src/hysop++/tests/testPolynoms/main.cpp b/HySoP/src/hysop++/tests/testPolynoms/main.cpp new file mode 100644 index 000000000..b6141085e --- /dev/null +++ b/HySoP/src/hysop++/tests/testPolynoms/main.cpp @@ -0,0 +1,10 @@ + +#include "gtest/gtest.h" + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + int ret = RUN_ALL_TESTS(); + + return ret; +} diff --git a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp new file mode 100644 index 000000000..4dd35271a --- /dev/null +++ b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp @@ -0,0 +1,103 @@ + +#include "testPolynoms.h" + +using namespace hysop::maths; + +template <typename T, std::size_t Dim> +void evalTest(std::size_t p_size, std::size_t p_samples); + +TEST_F(PolynomialTest, EvalTest1D) { + const std::size_t order=10; + const std::size_t samples=1024; + + evalTest<bool,1>(order,samples); + + evalTest<char,1>(order,samples); + evalTest<short,1>(order,samples); + evalTest<int,1>(order,samples); + evalTest<long int,1>(order,samples); + evalTest<long long int,1>(order,samples); + + evalTest<unsigned char,1>(order,samples); + evalTest<unsigned short,1>(order,samples); + evalTest<unsigned int,1>(order,samples); + evalTest<unsigned long int,1>(order,samples); + evalTest<unsigned long long int,1>(order,samples); + + evalTest<float,1>(order,samples); + evalTest<double,1>(order,samples); + evalTest<long double,1>(order,samples); + + evalTest<std::size_t,1>(order,samples); + evalTest<std::ptrdiff_t,1>(order,samples); +} + +TEST_F(PolynomialTest, EvalTest2D) { + const std::size_t order=10; + const std::size_t samples=128; + + evalTest<float,2>(order,samples); + evalTest<double,2>(order,samples); + evalTest<long double,2>(order,samples); +} + +TEST_F(PolynomialTest, EvalTest3D) { + const std::size_t order=10; + const std::size_t samples=16; + + evalTest<float,3>(order,samples); + evalTest<double,3>(order,samples); + evalTest<long double,3>(order,samples); +} + +TEST_F(PolynomialTest, EvalTest4D) { + const std::size_t order=10; + const std::size_t samples=4; + + evalTest<float,4>(order,samples); + evalTest<double,4>(order,samples); + evalTest<long double,4>(order,samples); +} + +template <typename T, std::size_t Dim> +void evalTest(const std::size_t p_size, const std::size_t p_samples) { + Polynomial<T,Dim> P; + Index<Dim> polyIdx; + { + Shape<Dim> polyShape; + polyShape.fill(p_size); + P.reshape(polyShape).applyToCoefficients([](T& ak, const Index<Dim>& idx){ + ak = T(idx())/T(idx.maxId()); + }); + } + + std::array<T,Dim> X; + T dX; + { + const T a = T(0); + const T b = T(1); + dX = (b-a)/(p_samples-1); + } + + Shape<Dim> sampleShape; + sampleShape.fill(p_samples); + Index<Dim> sampleIdx(sampleShape); + while(!sampleIdx.atMaxId()) { + for (std::size_t d=0; d < Dim; d++) + X[d] = sampleIdx[d]*dX; + T lhs, rhs; + lhs = P(X); + rhs = T(0); + polyIdx.reset(P.shape()); + while(!polyIdx.atMaxId()) { + T val = T(1); + for (std::size_t d=0; d<Dim; d++) + val *= std::pow(X[d],polyIdx[d]); + rhs += T(polyIdx())/T(polyIdx.maxId())*val; + ++polyIdx; + } + ASSERT_LE(std::abs(rhs-lhs),std::pow(10,Dim)*std::numeric_limits<T>::epsilon()); + ++sampleIdx; + } +} + diff --git a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.h b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.h new file mode 100644 index 000000000..d7dd0bd6a --- /dev/null +++ b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.h @@ -0,0 +1,13 @@ + +#include "gtest/gtest.h" +#include "maths/polynomial.h" + +using namespace hysop; + +class PolynomialTest : public ::testing::Test { + protected: + PolynomialTest() {} + void SetUp() {} + void TearDown() {} + virtual ~PolynomialTest() {} +}; -- GitLab