From b7679ed67fb8593ebd0479acd96382f8cfbd46e2 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Fri, 26 Feb 2016 20:11:18 +0100
Subject: [PATCH] Added missing main source (planner.cpp), added rpath and
macro define parsing in setup.py.in, extended swig coverage
---
.gitignore | 2 +-
HySoP/.gitignore | 5 +
HySoP/CMakeLists.txt | 4 +-
HySoP/setup.py.in | 7 +-
HySoP/src/hysop++/main/diffSolver.cpp | 4 +-
HySoP/src/hysop++/main/planner.cpp | 197 ++++
HySoP/src/hysop++/main/poissonSolver.cpp | 16 +-
.../data/multi_array/multi_array_defines.h | 6 +-
.../src/data/multi_array/multi_array_impl.h | 12 +-
HySoP/src/hysop++/src/domain/boundary.cpp | 2 +-
HySoP/src/hysop++/src/domain/boundary.h | 8 +-
HySoP/src/hysop++/src/domain/domain.h | 14 +-
.../hysop++/src/domain/domainConfiguration.h | 28 +-
HySoP/src/hysop++/src/fft/extension.cpp | 2 +-
HySoP/src/hysop++/src/fft/extension.h | 2 +-
.../hysop++/src/fft/fftDomainConfiguration.h | 8 +-
HySoP/src/hysop++/src/fft/fftw3.cpp | 40 +
HySoP/src/hysop++/src/fft/fftw3.h | 839 ++++++++++--------
HySoP/src/hysop++/src/fft/planner.h | 22 +-
HySoP/src/hysop++/src/maths/quad_maths.h | 8 +-
HySoP/src/hysop++/src/utils/types.h | 18 +-
.../tests/testDiffSolver/testDiffSolver.cpp | 4 +-
.../hysop++/tests/testPlanner/testPlanner.cpp | 4 +-
.../testPoissonSolver/testPoissonSolver.cpp | 4 +-
.../tests/testPolynoms/testPolynoms.cpp | 4 +-
HySoP/swig/constants.i | 7 -
HySoP/swig/cpp2hysop.i | 5 +-
HySoP/swig/{ => fftw}/fftw.i | 3 +-
HySoP/swig/hysop++/domain.i | 38 +
HySoP/swig/hysop++/fft.i | 91 ++
HySoP/swig/hysop++/hysop++.i | 5 +
HySoP/swig/hysop++/utils.i | 29 +
HySoP/swig/start.i | 4 +
HySoP/swig/types.i | 7 -
34 files changed, 970 insertions(+), 479 deletions(-)
create mode 100644 HySoP/.gitignore
create mode 100644 HySoP/src/hysop++/main/planner.cpp
create mode 100644 HySoP/src/hysop++/src/fft/fftw3.cpp
delete mode 100644 HySoP/swig/constants.i
rename HySoP/swig/{ => fftw}/fftw.i (62%)
create mode 100644 HySoP/swig/hysop++/domain.i
create mode 100644 HySoP/swig/hysop++/fft.i
create mode 100644 HySoP/swig/hysop++/hysop++.i
create mode 100644 HySoP/swig/hysop++/utils.i
delete mode 100644 HySoP/swig/types.i
diff --git a/.gitignore b/.gitignore
index aa4dfb4ff..9a8bdc053 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,4 @@
HySoP/hysop/__init__.py
.#*
HySoP/hysop/.pyflymakercc
-.DS_Store
\ No newline at end of file
+.DS_Store
diff --git a/HySoP/.gitignore b/HySoP/.gitignore
new file mode 100644
index 000000000..c4bde9e9b
--- /dev/null
+++ b/HySoP/.gitignore
@@ -0,0 +1,5 @@
+build/
+debug/
+release/
+.cache/
+hysop/f2hysop.pyf
diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index d5c833b5a..c79770845 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -252,10 +252,10 @@ set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}")
set(CXX_FLAGS "${CMAKE_CXX_FLAGS}")
set(CXX_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}")
-set(CXX_EXT_INCLUDES ${PYTHON_INCLUDE_DIR} ${FFTW_INCLUDES})
+set(CXX_EXT_INCLUDES ${PYTHON_INCLUDE_DIR} ${FFTW_INCLUDE_DIRS})
set(CXX_EXT_LIBS ${PYTHON_LIBRARIES} ${FFTW_LIBRARIES})
set(CXX_EXT_LIB_DIRS ${FFTW_LIBRARY_DIRS})
-set(CXX_EXTRA_DEFINES ${FFTW_DEFINES})
+set(CXX_EXTRA_DEFINES ${FFTW_DEFINES} -DHAS_EXTERN_TEMPLATES)
#swig package name (lib name generated by swig)
set(CPP_2_HYSOP "cpp2hysop")
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 3f14efafa..b68a69ce7 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -28,7 +28,7 @@ def parseCMakeDefines(var):
if defines == None:
return None
- # regex to match defines like -DMACRO_NAME = MACRO_VALUE
+ # regex to match compiler defines like -DMACRO_NAME or -DMACRO_NAME = MACRO_VALUE
p = re.compile('\s*(?:-D)?\s*(\w+)(?:\s*=\s*(\w+))?\s*')
res = list()
@@ -116,8 +116,8 @@ def create_swig_extension(name, inc_dirs, src_dirs=None, sources=None):
name = 'hysop._' + name
swig_opts = ['-I' + swig_dir,
- '-c++', '-modern']
- #-outdir', '/Users/Franck/toto',
+ '-O', '-Wextra', '-Werror',
+ '-c++', '-extranative', '-safecstrings']
extern_includes = parseCMakeVar("@CXX_EXT_INCLUDES@")
if(extern_includes != None):
@@ -147,6 +147,7 @@ def create_swig_extension(name, inc_dirs, src_dirs=None, sources=None):
library_dirs=library_dirs,
libraries=libraries,
define_macros=define_macros,
+ runtime_library_dirs=library_dirs,
extra_compile_args=extra_compile_args,
extra_link_args=extra_link_args)
diff --git a/HySoP/src/hysop++/main/diffSolver.cpp b/HySoP/src/hysop++/main/diffSolver.cpp
index 29548ea56..3d0316e3f 100644
--- a/HySoP/src/hysop++/main/diffSolver.cpp
+++ b/HySoP/src/hysop++/main/diffSolver.cpp
@@ -69,7 +69,7 @@ std::function<T(T)> derivative(std::size_t k, int order) {
template <typename T, std::size_t Dim, bool verbose=false>
void test(std::size_t p_maxOrder, bool includePeriodicBds=false) {
- Shape<Dim> shape;
+ typename Shape<Dim>::type shape;
typename Domain<T,Dim>::DomainSize domainSize;
Domain<T,Dim> ref, inBuffer, outBuffer;
@@ -88,7 +88,7 @@ void test(std::size_t p_maxOrder, bool includePeriodicBds=false) {
in = ref;
out = ref;
- Shape<Dim> maxOrder, testCases;
+ typename Shape<Dim>::type maxOrder, testCases;
maxOrder.fill(p_maxOrder+1);
testCases.fill(nExtensionsPair);
Index<Dim> orderId(maxOrder);
diff --git a/HySoP/src/hysop++/main/planner.cpp b/HySoP/src/hysop++/main/planner.cpp
new file mode 100644
index 000000000..287c6713c
--- /dev/null
+++ b/HySoP/src/hysop++/main/planner.cpp
@@ -0,0 +1,197 @@
+
+#include "maths/quad_maths.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);
+
+
+int main(int argc, const char *argv[]) {
+
+#ifdef FFTW_HAS_FFTW3F
+ 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);
+#endif
+
+#ifdef FFTW_HAS_FFTW3D
+ 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);
+#endif
+
+#ifdef FFTW_HAS_FFTW3L
+ 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);
+#endif
+
+#ifdef FFTW_HAS_FFTW3Q
+ std::cout << std::endl;
+ std::cout << "== TEST 1D - __float128 ==\t";
+ test<__float128,1>(false,false);
+ std::cout << "== TEST 2D - __float128 ==\t";
+ test<__float128,2>(false,false);
+ std::cout << "== TEST 3D - __float128 ==\t";
+ test<__float128,3>(false,false);
+#endif
+
+ return 0;
+}
+
+template <typename T, std::size_t Dim, bool verbose>
+void test(bool inplace, bool includePeriodicBds) {
+ typename Shape<Dim>::type 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;
+
+ typename Shape<Dim>::type 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 */
+ bool status = planner.plan(in.data(), out.data(), domainConfig, order, domainSize, FFTW_MEASURE,
+ includePeriodicBds, includePeriodicBds);
+ assert(status || 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++/main/poissonSolver.cpp b/HySoP/src/hysop++/main/poissonSolver.cpp
index dfb829a49..e95a8cd50 100644
--- a/HySoP/src/hysop++/main/poissonSolver.cpp
+++ b/HySoP/src/hysop++/main/poissonSolver.cpp
@@ -45,19 +45,19 @@ std::function<T(T)> func(std::size_t k) {
std::string bdsToStr(domain::Boundary bd) {
switch(bd) {
- case(NONE): return "None ";
- case(PERIODIC) : return "Periodic ";
- case(HOMOGENEOUS_NEUMANN): return "Hom_Neum. ";
- case(HOMOGENEOUS_DIRICHLET): return "Hom_Diric.";
- case(NEUMANN): return "Neumann ";
- case(DIRICHLET): return "Dirichlet ";
+ case(Boundary::NONE): return "None ";
+ case(Boundary::PERIODIC) : return "Periodic ";
+ case(Boundary::HOMOGENEOUS_NEUMANN): return "Hom_Neum. ";
+ case(Boundary::HOMOGENEOUS_DIRICHLET): return "Hom_Diric.";
+ case(Boundary::NEUMANN): return "Neumann ";
+ case(Boundary::DIRICHLET): return "Dirichlet ";
}
return "";
}
template <typename T, std::size_t Dim, bool verbose=false>
void test(bool includePeriodicBds=false) {
- Shape<Dim> shape;
+ typename Shape<Dim>::type shape;
typename Domain<T,Dim>::DomainSize domainSize;
Domain<T,Dim> ref, inBuffer, outBuffer;
@@ -74,7 +74,7 @@ void test(bool includePeriodicBds=false) {
in = ref;
out = ref;
- Shape<Dim> testCases;
+ typename Shape<Dim>::type testCases;
testCases.fill(nBoundaryPairs);
Index<Dim> testCaseId(testCases);
std::array<T,3> meanDists{0};
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
index 75e51a227..cc1b03018 100644
--- a/HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h
+++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_defines.h
@@ -60,7 +60,7 @@
/* public const common interfaces for views and referencess */
#define PUBLIC_COMMON_CONST_INTERFACE(TYPENAME) \
/* shape related */ \
- Shape<Dim> shape() const; \
+ typename Shape<Dim>::type shape() const; \
\
/* print data */ \
const TYPENAME& print(const std::string& name, std::ostream& os = std::cout, \
@@ -146,8 +146,8 @@
#define COMMON_CONST_IMPLEMENTATION(TYPENAME,TEMPLATES) \
/* shape related */ \
TEMPLATES \
- Shape<Dim> TYPENAME::shape() const { \
- Shape<Dim> shape; \
+ typename Shape<Dim>::type TYPENAME::shape() const { \
+ typename Shape<Dim>::type 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]); \
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
index 870595368..efb1b1c3b 100644
--- a/HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h
+++ b/HySoP/src/hysop++/src/data/multi_array/multi_array_impl.h
@@ -21,7 +21,7 @@ namespace hysop {
public:
multi_array(const extents_gen<Dim>& extents = extents_gen<Dim>());
- multi_array(const Shape<Dim>& shape);
+ multi_array(const typename Shape<Dim>::type& shape);
multi_array(const multi_array& other);
multi_array(multi_array&& other);
@@ -52,10 +52,10 @@ namespace hysop {
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);
+ multi_array& reshape(const typename Shape<Dim>::type& shape);
protected:
- static extents_gen<Dim> shapeToExtents(const Shape<Dim> &shape);
+ static extents_gen<Dim> shapeToExtents(const typename Shape<Dim>::type &shape);
};
@@ -67,7 +67,7 @@ namespace hysop {
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):
+ multi_array<T,Dim,Allocator>::multi_array(const typename Shape<Dim>::type& shape):
boost_multi_array<T,Dim,Allocator>(shapeToExtents(shape)) {}
template <typename T, std::size_t Dim, typename Allocator>
@@ -165,7 +165,7 @@ namespace hysop {
/* 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) {
+ typename multi_array<T,Dim,Allocator>::template extents_gen<Dim> multi_array<T,Dim,Allocator>::shapeToExtents(const typename Shape<Dim>::type &shape) {
return utils::buildExtents(shape);
}
@@ -175,7 +175,7 @@ namespace hysop {
/* 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) {
+ multi_array<T,Dim,Allocator>& multi_array<T,Dim,Allocator>::reshape(const typename Shape<Dim>::type& shape) {
boost::array<int,Dim> extents;
for (std::size_t d = 0; d < Dim; d++)
extents[d] = static_cast<int>(shape[d]);
diff --git a/HySoP/src/hysop++/src/domain/boundary.cpp b/HySoP/src/hysop++/src/domain/boundary.cpp
index 4cb3f491b..4a3f862af 100644
--- a/HySoP/src/hysop++/src/domain/boundary.cpp
+++ b/HySoP/src/hysop++/src/domain/boundary.cpp
@@ -14,7 +14,7 @@ namespace hysop {
};
const char* toStringBoundary(Boundary bd) {
- return boundary_strings[bd+1];
+ return boundary_strings[static_cast<int>(bd)+1];
}
std::ostream& operator<<(std::ostream& os, const Boundary& bd) {
diff --git a/HySoP/src/hysop++/src/domain/boundary.h b/HySoP/src/hysop++/src/domain/boundary.h
index d0aa55438..9320085e0 100644
--- a/HySoP/src/hysop++/src/domain/boundary.h
+++ b/HySoP/src/hysop++/src/domain/boundary.h
@@ -1,12 +1,12 @@
-#ifndef BOUNDARY_H
-#define BOUNDARY_H
+#ifndef HYSOP_BOUNDARY_H
+#define HYSOP_BOUNDARY_H
#include <iostream>
namespace hysop {
namespace domain {
- enum Boundary : int {
+ enum class Boundary : int {
NONE = -1,
HOMOGENEOUS_DIRICHLET = 0,
HOMOGENEOUS_NEUMANN = 1,
@@ -20,4 +20,4 @@ namespace hysop {
}
}
-#endif /* end of include guard: BOUNDARY_H */
+#endif /* end of include guard: HYSOP_BOUNDARY_H */
diff --git a/HySoP/src/hysop++/src/domain/domain.h b/HySoP/src/hysop++/src/domain/domain.h
index 9a9b0d27c..70d860ece 100644
--- a/HySoP/src/hysop++/src/domain/domain.h
+++ b/HySoP/src/hysop++/src/domain/domain.h
@@ -28,7 +28,7 @@ namespace hysop {
m_shape{0}, m_dataShape{0}, m_leftDataOffset{0},
m_domainConfig(), m_domainSize{0}, m_spaceStep{0}, m_data() {}
- Domain(const Shape<Dim>& p_shape, const DomainConfiguration<Dim>& p_domainConfig, const DomainSize& p_domainSize) :
+ Domain(const typename Shape<Dim>::type& 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}, m_data() {
@@ -39,7 +39,7 @@ namespace hysop {
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) {
+ Domain<T,Dim>& reshape(const typename Shape<Dim>::type& p_domainShape) {
m_shape = p_domainShape;
m_domainConfig.getShapeConfiguration(m_shape, m_dataShape, m_leftDataOffset);
m_data.reshape(m_dataShape);
@@ -68,11 +68,11 @@ namespace hysop {
return fft::FftDomainConfiguration<Dim>(m_domainConfig);
}
- const Shape<Dim>& shape() const { return m_shape; }
- const Shape<Dim>& dataShape() const { return m_dataShape; }
+ const typename Shape<Dim>::type& shape() const { return m_shape; }
+ const typename Shape<Dim>::type& 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 typename Offset<Dim>::type& leftDataOffset() const { return m_leftDataOffset; }
const DomainConfiguration<Dim> boundaryConfig() const { return m_domainConfig; }
const hysop::multi_array<T,Dim>& data() const { return m_data; }
@@ -109,8 +109,8 @@ namespace hysop {
}
protected:
- Shape<Dim> m_shape, m_dataShape;
- Offset<Dim> m_leftDataOffset;
+ typename Shape<Dim>::type m_shape, m_dataShape;
+ typename Offset<Dim>::type m_leftDataOffset;
DomainConfiguration<Dim> m_domainConfig;
DomainSize m_domainSize;
diff --git a/HySoP/src/hysop++/src/domain/domainConfiguration.h b/HySoP/src/hysop++/src/domain/domainConfiguration.h
index 42af207fa..3f1c7ded7 100644
--- a/HySoP/src/hysop++/src/domain/domainConfiguration.h
+++ b/HySoP/src/hysop++/src/domain/domainConfiguration.h
@@ -22,8 +22,8 @@ namespace hysop {
class DomainConfiguration {
public:
- using BoundaryPair = std::pair<Boundary,Boundary>;
- using BoundaryArray = std::array<BoundaryPair,Dim>;
+ typedef std::pair<Boundary,Boundary> BoundaryPair;
+ typedef std::array<BoundaryPair,Dim> BoundaryArray;
public:
DomainConfiguration(const BoundaryArray& p_boundaries = defaultDomainBoundaries(),
//bool p_includeDirichletBoundaries = true,
@@ -40,14 +40,14 @@ namespace hysop {
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 {
+ void getShapeConfiguration(const typename Shape<Dim>::type &p_fullShape, typename Shape<Dim>::type &p_realShape, typename Offset<Dim>::type &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);
+ bool hasPeriodicRightOffset = (pair.second == Boundary::PERIODIC);
std::size_t periodicRightOffset = hasPeriodicRightOffset && !this->includePeriodicBoundaries();
std::size_t leftOffset = 0;
std::size_t rightOffset = periodicRightOffset;
@@ -62,19 +62,21 @@ namespace hysop {
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))
+ if((pair.first == Boundary::PERIODIC) ^ (pair.second == Boundary::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>());
- }
+ #ifndef SWIG
+ 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>());
+ }
+ #endif
protected:
BoundaryArray m_boundaries;
diff --git a/HySoP/src/hysop++/src/fft/extension.cpp b/HySoP/src/hysop++/src/fft/extension.cpp
index eeacb339b..6f2a67e5e 100644
--- a/HySoP/src/hysop++/src/fft/extension.cpp
+++ b/HySoP/src/hysop++/src/fft/extension.cpp
@@ -12,7 +12,7 @@ namespace hysop {
};
const char* toStringExtension(Extension ext) {
- return extension_strings[ext+1];
+ return extension_strings[static_cast<int>(ext)+1];
}
std::ostream& operator<<(std::ostream& os, const Extension& ext) {
diff --git a/HySoP/src/hysop++/src/fft/extension.h b/HySoP/src/hysop++/src/fft/extension.h
index 815f25a65..5aaba36f2 100644
--- a/HySoP/src/hysop++/src/fft/extension.h
+++ b/HySoP/src/hysop++/src/fft/extension.h
@@ -7,7 +7,7 @@
namespace hysop {
namespace fft {
- enum Extension : int {
+ enum class Extension : int {
NONE=-1,
EVEN=0,
ODD=1,
diff --git a/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h b/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h
index f3ccea896..57f014d63 100644
--- a/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h
+++ b/HySoP/src/hysop++/src/fft/fftDomainConfiguration.h
@@ -19,10 +19,10 @@ namespace hysop {
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>;
+ typedef typename domain::DomainConfiguration<Dim>::BoundaryPair BoundaryPair;
+ typedef typename domain::DomainConfiguration<Dim>::BoundaryArray BoundaryArray;
+ typedef std::pair<fft::Extension,fft::Extension> ExtensionPair;
+ typedef std::array<ExtensionPair,Dim> ExtensionArray;
public:
FftDomainConfiguration(const ExtensionArray& p_extensions, bool p_includePeriodicBoundaries) :
m_extensions(p_extensions), m_includePeriodicBoundaries(p_includePeriodicBoundaries) {
diff --git a/HySoP/src/hysop++/src/fft/fftw3.cpp b/HySoP/src/hysop++/src/fft/fftw3.cpp
new file mode 100644
index 000000000..25a877c87
--- /dev/null
+++ b/HySoP/src/hysop++/src/fft/fftw3.cpp
@@ -0,0 +1,40 @@
+
+#include "fft/fftw3.h"
+
+namespace hysop {
+ namespace fft {
+ template struct Fftw3<float>;
+ template struct Fftw3<double>;
+ template struct Fftw3<long double>;
+ #ifdef HAS_QUADMATHS
+ template struct Fftw3<__float128>;
+ #endif
+ }
+}
+
+/* needed for swig to avoid undefined reference link error */
+#if !defined(FFTW_HAS_FFTW3F_THREADS) || !defined(FFTW_HAS_FFTW3F_OMP)
+ int fftwf_init_threads() { return 0; }
+ void fftwf_plan_with_nthreads(int) {}
+ void fftwf_cleanup_threads() {}
+#endif
+
+#if !defined(FFTW_HAS_FFTW3D_THREADS) || !defined(FFTW_HAS_FFTW3D_OMP)
+ int fftw_init_threads() { return 0; }
+ void fftw_plan_with_nthreads(int) {}
+ void fftw_cleanup_threads() {}
+#endif
+
+#if !defined(FFTW_HAS_FFTW3L_THREADS) || !defined(FFTW_HAS_FFTW3L_OMP)
+ int fftwl_init_threads() { return 0; }
+ void fftwl_plan_with_nthreads(int) {}
+ void fftwl_cleanup_threads() {}
+#endif
+
+#ifdef HAS_QUADMATHS
+ #if !defined(FFTW_HAS_FFTW3Q_THREADS) || !defined(FFTW_HAS_FFTW3Q_OMP)
+ int fftwq_init_threads() { return 0; }
+ void fftwq_plan_with_nthreads(int) {}
+ void fftwq_cleanup_threads() {}
+ #endif
+#endif
diff --git a/HySoP/src/hysop++/src/fft/fftw3.h b/HySoP/src/hysop++/src/fft/fftw3.h
index f08bca38b..56c5c1eac 100644
--- a/HySoP/src/hysop++/src/fft/fftw3.h
+++ b/HySoP/src/hysop++/src/fft/fftw3.h
@@ -13,7 +13,6 @@
/* 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
@@ -36,384 +35,459 @@
#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); \
- } \
+#define FFTW_DEFINE_CXX_API(X, Y, REAL, has_thread_support) \
+template <> \
+ struct Fftw3<REAL> { \
+ typedef REAL R; \
+ typedef Y(complex) C; \
+ typedef Y(plan) plan; \
+ typedef Y(iodim) iodim; \
+ typedef Y(iodim64) iodim64; \
+ typedef Y(r2r_kind) r2r_kind; \
+ typedef Y(read_char_func) read_char_func; \
+ typedef Y(write_char_func) write_char_func; \
+ typedef std::complex<REAL> stdC; \
+ \
+ 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); \
+ } \
+ \
+ template <typename T = void> \
+ typename std::enable_if<has_thread_support, T>::type \
+ X(plan_with_nthreads)(int nthreads) const { \
+ Y(plan_with_nthreads)(nthreads); \
+ } \
+ template <typename T = int> \
+ typename std::enable_if<has_thread_support, T>::type \
+ X(init_threads)() const { \
+ return Y(init_threads)(); \
+ } \
+ template <typename T=void> \
+ typename std::enable_if<has_thread_support, T>::type \
+ 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); \
+ } \
};
+/* Constants */
namespace hysop {
namespace fft {
+ /* float */
+ #ifdef FFTW_HAS_FFTW3F
+ static constexpr bool fftw_has_float_support = true;
+ #ifdef FFTW_HAS_FFTW3F_MPI
+ static constexpr bool fftw_has_float_mpi_support = true;
+ #else
+ static constexpr bool fftw_has_float_mpi_support = false;
+ #endif
+ #if defined(FFTW_HAS_FFTW3F_THREADS) || defined(FFTW_HAS_FFTW3F_OMP)
+ static constexpr bool fftw_has_float_thread_support = true;
+ #else
+ static constexpr bool fftw_has_float_thread_support = false;
+ #endif
+ #else
+ static constexpr bool fftw_has_float_support = false;
+ static constexpr bool fftw_has_float_thread_support = false;
+ static constexpr bool fftw_has_float_mpi_support = false;
+ #endif
- 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; };
+ /* double */
+ #ifdef FFTW_HAS_FFTW3D
+ static constexpr bool fftw_has_double_support = true;
+ #ifdef FFTW_HAS_FFTW3D_MPI
+ static constexpr bool fftw_has_double_mpi_support = true;
+ #else
+ static constexpr bool fftw_has_double_mpi_support = false;
+ #endif
+ #if defined(FFTW_HAS_FFTW3D_THREADS) || defined(FFTW_HAS_FFTW3D_OMP)
+ static constexpr bool fftw_has_double_thread_support = true;
+ #else
+ static constexpr bool fftw_has_double_thread_support = false;
+ #endif
+ #else
+ static constexpr bool fftw_has_double_support = false;
+ static constexpr bool fftw_has_double_thread_support = false;
+ static constexpr bool fftw_has_double_mpi_support = false;
+ #endif
- template <typename T> struct is_fftw_supported_complex_type { static constexpr bool value = false; };
- template <> struct is_fftw_supported_complex_type<std::complex<float>> { static constexpr bool value = true; };
- template <> struct is_fftw_supported_complex_type<std::complex<double>> { static constexpr bool value = true; };
- template <> struct is_fftw_supported_complex_type<std::complex<long double>> { static constexpr bool value = true; };
- template <> struct is_fftw_supported_complex_type<fftwf_complex> { static constexpr bool value = true; };
- template <> struct is_fftw_supported_complex_type<fftw_complex> { static constexpr bool value = true; };
- template <> struct is_fftw_supported_complex_type<fftwl_complex> { static constexpr bool value = true; };
+ /* long double */
+ #ifdef FFTW_HAS_FFTW3L
+ static constexpr bool fftw_has_long_double_support = true;
+ #ifdef FFTW_HAS_FFTW3L_MPI
+ static constexpr bool fftw_has_long_double_mpi_support = true;
+ #else
+ static constexpr bool fftw_has_long_double_mpi_support = false;
+ #endif
+ #if defined(FFTW_HAS_FFTW3L_THREADS) || defined(FFTW_HAS_FFTW3L_OMP)
+ static constexpr bool fftw_has_long_double_thread_support = true;
+ #else
+ static constexpr bool fftw_has_long_double_thread_support = false;
+ #endif
+ #else
+ static constexpr bool fftw_has_long_double_support = false;
+ static constexpr bool fftw_has_long_double_thread_support = false;
+ static constexpr bool fftw_has_long_double_mpi_support = false;
+ #endif
+
+ /* __float128 */
+ #ifdef FFTW_HAS_FFTW3Q
+ static constexpr bool fftw_has_quad_float_support = true;
+ #ifdef FFTW_HAS_FFTW3Q_MPI
+ static constexpr bool fftw_has_quad_float_mpi_support = true;
+ #else
+ static constexpr bool fftw_has_quad_float_mpi_support = false;
+ #endif
+ #if defined(FFTW_HAS_FFTW3Q_THREADS) || defined(FFTW_HAS_FFTW3Q_OMP)
+ static constexpr bool fftw_has_quad_float_thread_support = true;
+ #else
+ static constexpr bool fftw_has_quad_float_thread_support = false;
+ #endif
+ #else
+ static constexpr bool fftw_has_quad_float_support = false;
+ static constexpr bool fftw_has_quad_float_thread_support = false;
+ static constexpr bool fftw_has_quad_float_mpi_support = false;
+ #endif
+ }
+}
+
+/* Wrappers */
+namespace hysop {
+ namespace fft {
template <typename T>
- struct Fftw3 {
+ struct Fftw3 {
Fftw3() {
throw std::runtime_error(
"Can only use Fftw3 wrapper with types {float, double, long double, __float128} ! "
@@ -421,18 +495,33 @@ namespace hysop {
);
}
};
-
+ template <typename T> struct is_fftw_supported_type { static constexpr bool value = false; };
+ template <typename T> struct is_fftw_supported_complex_type { static constexpr bool value = false; };
+
+
/* Generate Ffftw<> template specialisations */
- 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)
+ FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_FLOAT, float, hysop::fft::fftw_has_float_thread_support)
+ FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_DOUBLE, double, hysop::fft::fftw_has_double_thread_support)
+ FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_LONG_DOUBLE, long double, hysop::fft::fftw_has_long_double_thread_support)
+
+ 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; };
+
+ template <> struct is_fftw_supported_complex_type<std::complex<float>> { static constexpr bool value = true; };
+ template <> struct is_fftw_supported_complex_type<std::complex<double>> { static constexpr bool value = true; };
+ template <> struct is_fftw_supported_complex_type<std::complex<long double>> { static constexpr bool value = true; };
+ template <> struct is_fftw_supported_complex_type<fftwf_complex> { static constexpr bool value = true; };
+ template <> struct is_fftw_supported_complex_type<fftw_complex> { static constexpr bool value = true; };
+ template <> struct is_fftw_supported_complex_type<fftwl_complex> { static constexpr bool value = true; };
#ifdef HAS_QUADMATHS
+ FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_QUAD, __float128, hysop::fft::fftw_has_quad_float_thread_support)
template <> struct is_fftw_supported_type<__float128> { static constexpr bool value = true; };
template <> struct is_fftw_supported_complex_type<std::complex<__float128>> { static constexpr bool value = true; };
template <> struct is_fftw_supported_complex_type<fftwq_complex> { static constexpr bool value = true; };
- FFTW_DEFINE_CXX_API(FFTW_MANGLE_CLASS, FFTW_MANGLE_QUAD, __float128)
#endif
+
}
}
diff --git a/HySoP/src/hysop++/src/fft/planner.h b/HySoP/src/hysop++/src/fft/planner.h
index f9b5f9e0e..75b8ca852 100644
--- a/HySoP/src/hysop++/src/fft/planner.h
+++ b/HySoP/src/hysop++/src/fft/planner.h
@@ -148,7 +148,7 @@ namespace hysop {
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});
+ m_complexBuffer.reshape(typename Shape<Dim>::type{0});
for(fftw_plan plan : m_forward_R2R_plans)
this->fftw_destroy_plan(plan);
for(fftw_plan plan : m_backward_R2R_plans)
@@ -208,7 +208,7 @@ namespace hysop {
m_mirrorOutputPeriodicBoundaries = includeOutputPeriodicBoundaries && mirrorOutputPeriodicBoundaries;
TransformArray forwardTransforms, backwardTransforms;
- Shape<Dim> realBufferShape, complexBufferShape;
+ typename Shape<Dim>::type 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;
@@ -268,10 +268,10 @@ namespace hysop {
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);
+ const std::size_t firstExtOdd = (inputExtensions.first == Extension::ODD);
+ const std::size_t secondExtOdd = (inputExtensions.second == Extension::ODD);
+ const std::size_t firstExtEven = (inputExtensions.first == Extension::EVEN);
+ const std::size_t secondExtEven = (inputExtensions.second == Extension::EVEN);
forward_transform_size[d] -= (firstExtOdd + secondExtOdd);
forward_input_offset[d] = std::ptrdiff_t(firstExtOdd);
@@ -602,7 +602,7 @@ namespace hysop {
this->fftw_execute(*plan);
if(m_mirrorOutputPeriodicBoundaries) {
- const Shape<Dim> rshape = m_realBuffer.shape();
+ const typename Shape<Dim>::type rshape = m_realBuffer.shape();
T *data = m_realBuffer.rdata();
for(const int axe : m_periodicAxes) {
const int N = rshape[axe];
@@ -678,14 +678,14 @@ namespace hysop {
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)
+ if(ed.first == Extension::EVEN) {
+ if(ed.second == Extension::EVEN)
return fft::Transform(FFTW_REDFT00);
else
return fft::Transform(FFTW_REDFT01);
}
- else if(ed.first == ODD) {
- if(ed.second == EVEN)
+ else if(ed.first == Extension::ODD) {
+ if(ed.second == Extension::EVEN)
return fft::Transform(FFTW_RODFT01);
else
return fft::Transform(FFTW_RODFT00);
diff --git a/HySoP/src/hysop++/src/maths/quad_maths.h b/HySoP/src/hysop++/src/maths/quad_maths.h
index a411d8041..97ae1fe0b 100644
--- a/HySoP/src/hysop++/src/maths/quad_maths.h
+++ b/HySoP/src/hysop++/src/maths/quad_maths.h
@@ -179,15 +179,15 @@ namespace std {
__complex128 X{x.real(),x.imag()};
return cabsq(X);
}
- inline __float128 carg(std::complex<__float128> x) {
+ inline __float128 arg(std::complex<__float128> x) {
__complex128 X{x.real(),x.imag()};
return cargq(X);
}
- inline __float128 cimag(std::complex<__float128> x) {
+ inline __float128 imag(std::complex<__float128> x) {
__complex128 X{x.real(),x.imag()};
return cimagq(X);
}
- inline __float128 creal(std::complex<__float128> x) {
+ inline __float128 real(std::complex<__float128> x) {
__complex128 X{x.real(),x.imag()};
return crealq(X);
}
@@ -288,7 +288,7 @@ namespace std {
}
inline std::complex< __float128 > pow(std::complex< __float128> x , __float128 y) {
- __float128 R = powq(abs(x), y);
+ __float128 R = powq(std::abs(x), y);
__float128 phi = atanq(x.imag()/x.real());
return std::complex<__float128 >(R*cosq(y*phi), R*sinq(y*phi));
}
diff --git a/HySoP/src/hysop++/src/utils/types.h b/HySoP/src/hysop++/src/utils/types.h
index 58dcdbc28..ccbd3c345 100644
--- a/HySoP/src/hysop++/src/utils/types.h
+++ b/HySoP/src/hysop++/src/utils/types.h
@@ -31,14 +31,18 @@ namespace hysop {
} /* end of namespace types */
- /* expose those types to namespace hysop */
- template <std::size_t Dim>
- using Shape = std::array<std::size_t, Dim>;
+/* expose the folowwing types to namespace hysop */
+
+/* swig does not support alias templates... */
+ template <std::size_t Dim>
+ struct Shape {
+ typedef std::array<std::size_t, Dim> type;
+ };
+ template <std::size_t Dim>
+ struct Offset {
+ typedef std::array<std::ptrdiff_t, Dim> type;
+ };
- 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>;
diff --git a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
index 9144b1fdc..55fec62f5 100644
--- a/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
+++ b/HySoP/src/hysop++/tests/testDiffSolver/testDiffSolver.cpp
@@ -70,7 +70,7 @@ std::function<T(T)> derivative(std::size_t k, int order) {
template <typename T, std::size_t Dim, bool verbose=false>
void test(std::size_t p_maxOrder, bool includePeriodicBds=false) {
- Shape<Dim> shape;
+ typename Shape<Dim>::type shape;
typename Domain<T,Dim>::DomainSize domainSize;
Domain<T,Dim> ref, inBuffer, outBuffer;
@@ -89,7 +89,7 @@ void test(std::size_t p_maxOrder, bool includePeriodicBds=false) {
in = ref;
out = ref;
- Shape<Dim> maxOrder, testCases;
+ typename Shape<Dim>::type maxOrder, testCases;
maxOrder.fill(p_maxOrder+1);
testCases.fill(nExtensionsPair);
Index<Dim> orderId(maxOrder);
diff --git a/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp b/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp
index 09311df6a..aab308e98 100644
--- a/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp
+++ b/HySoP/src/hysop++/tests/testPlanner/testPlanner.cpp
@@ -78,7 +78,7 @@ TEST_F(PlannerTest, InplaceQuadDoubleTransforms) {
template <typename T, std::size_t Dim, bool verbose>
void test(bool inplace, bool includePeriodicBds) {
- Shape<Dim> shape;
+ typename Shape<Dim>::type shape;
typename Domain<T,Dim>::DomainSize domainSize;
Domain<T,Dim> ref, inBuffer, outBuffer;
@@ -101,7 +101,7 @@ void test(bool inplace, bool includePeriodicBds) {
in = ref;
out = ref;
- Shape<Dim> testCases;
+ typename Shape<Dim>::type testCases;
testCases.fill(nExtensionsPair);
Index<Dim> testCaseId(testCases);
std::array<T,3> meanDists{0};
diff --git a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
index 9711a4cba..58d87d544 100644
--- a/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
+++ b/HySoP/src/hysop++/tests/testPoissonSolver/testPoissonSolver.cpp
@@ -47,7 +47,7 @@ std::function<T(T)> func(std::size_t k) {
template <typename T, std::size_t Dim, bool verbose=false>
void test(bool includePeriodicBds=false) {
- Shape<Dim> shape;
+ typename Shape<Dim>::type shape;
typename Domain<T,Dim>::DomainSize domainSize;
Domain<T,Dim> ref, inBuffer, outBuffer;
@@ -64,7 +64,7 @@ void test(bool includePeriodicBds=false) {
in = ref;
out = ref;
- Shape<Dim> testCases;
+ typename Shape<Dim>::type testCases;
testCases.fill(nBoundaryPairs);
Index<Dim> testCaseId(testCases);
std::array<T,3> meanDists{0};
diff --git a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
index e4001fa85..0849af6a4 100644
--- a/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
+++ b/HySoP/src/hysop++/tests/testPolynoms/testPolynoms.cpp
@@ -76,7 +76,7 @@ void evalTest(const std::size_t p_size, const std::size_t p_samples) {
Polynomial<T,Dim> P;
Index<Dim> polyIdx;
{
- Shape<Dim> polyShape;
+ typename Shape<Dim>::type polyShape;
polyShape.fill(p_size);
P.reshape(polyShape).applyToCoefficients([](T& ak, const Index<Dim>& idx){
ak = T(idx())/T(idx.maxId());
@@ -91,7 +91,7 @@ void evalTest(const std::size_t p_size, const std::size_t p_samples) {
dX = (b-a)/(p_samples-1);
}
- Shape<Dim> sampleShape;
+ typename Shape<Dim>::type sampleShape;
sampleShape.fill(p_samples);
Index<Dim> sampleIdx(sampleShape);
while(!sampleIdx.atMaxId()) {
diff --git a/HySoP/swig/constants.i b/HySoP/swig/constants.i
deleted file mode 100644
index 929176357..000000000
--- a/HySoP/swig/constants.i
+++ /dev/null
@@ -1,7 +0,0 @@
-
-// -*- C++ -*-
-%{
- #include "utils/constants.h"
-%}
-
-%include "utils/constants.h"
diff --git a/HySoP/swig/cpp2hysop.i b/HySoP/swig/cpp2hysop.i
index abc7b82a9..4a2ead8af 100644
--- a/HySoP/swig/cpp2hysop.i
+++ b/HySoP/swig/cpp2hysop.i
@@ -3,7 +3,6 @@
%include start.i
-%include constants.i
-
-%include fftw.i
+%include fftw/fftw.i
+%include hysop++/hysop++.i
diff --git a/HySoP/swig/fftw.i b/HySoP/swig/fftw/fftw.i
similarity index 62%
rename from HySoP/swig/fftw.i
rename to HySoP/swig/fftw/fftw.i
index 912c3b467..28ab6fe64 100644
--- a/HySoP/swig/fftw.i
+++ b/HySoP/swig/fftw/fftw.i
@@ -1,6 +1,7 @@
+
// -*- C++ -*-
%{
- #include "fftw.hpp"
+ #include "fftw.hpp"
%}
%include "fftw.hpp"
diff --git a/HySoP/swig/hysop++/domain.i b/HySoP/swig/hysop++/domain.i
new file mode 100644
index 000000000..8e996563f
--- /dev/null
+++ b/HySoP/swig/hysop++/domain.i
@@ -0,0 +1,38 @@
+%module domain
+
+// -*- C++ -*-
+%{
+ #include "domain/boundary.h"
+ #include "domain/domainConfiguration.h"
+%}
+
+/* domain boundary */
+%ignore hysop::domain::operator<<;
+%include "domain/boundary.h"
+
+
+/* domain configuration */
+%ignore hysop::domain::operator<<;
+%ignore hysop::domain::DomainConfiguration::operator[];
+%include "domain/domainConfiguration.h"
+
+%define INSTANTIATE_DOMAIN_CONFIG_IMPL(CLASS_NAME,CLASS_TEMPLATE)
+ %template(CLASS_NAME) CLASS_TEMPLATE;
+ %extend CLASS_TEMPLATE {
+ CLASS_TEMPLATE::BoundaryPair __get_item__(unsigned int k) const {
+ return $self->operator[](k);
+ }
+ const char* __str__() const {
+ std::stringstream ss;
+ ss << *($self);
+ return ss.str().c_str();
+ }
+ }
+%enddef
+
+%define INSTANTIATE_DOMAIN_CONFIG(INT)
+ INSTANTIATE_DOMAIN_CONFIG_IMPL(DomainConfiguration ## INT ## D, hysop::domain::DomainConfiguration<INT>)
+%enddef
+
+%formacro(INSTANTIATE_DOMAIN_CONFIG, INSTANTIATED_DIMENSIONS)
+
diff --git a/HySoP/swig/hysop++/fft.i b/HySoP/swig/hysop++/fft.i
new file mode 100644
index 000000000..accd97cd9
--- /dev/null
+++ b/HySoP/swig/hysop++/fft.i
@@ -0,0 +1,91 @@
+
+%module fft
+
+// -*- C++ -*-
+%{
+ #include "fftw3.h"
+ #include "fft/fftw3.h"
+ #include "fft/extension.h"
+ #include "fft/transform.h"
+ #include "fft/fftDomainConfiguration.h"
+%}
+
+/* fftw3.h */
+#if !defined(FFTW_HAS_FFTW3F_THREADS) || !defined(FFTW_HAS_FFTW3F_OMP)
+ %ignore fftwf_init_threads();
+ %ignore fftwf_plan_with_nthreads(int);
+ %ignore fftwf_cleanup_threads();
+#endif
+
+#if !defined(FFTW_HAS_FFTW3D_THREADS) || !defined(FFTW_HAS_FFTW3D_OMP)
+ %ignore fftw_init_threads();
+ %ignore fftw_plan_with_nthreads(int);
+ %ignore fftw_cleanup_threads();
+#endif
+
+#if !defined(FFTW_HAS_FFTW3L_THREADS) || !defined(FFTW_HAS_FFTW3L_OMP)
+ %ignore fftwl_init_threads();
+ %ignore fftwl_plan_with_nthreads(int);
+ %ignore fftwl_cleanup_threads();
+#endif
+
+#ifdef HAS_QUADMATHS
+ #if !defined(FFTW_HAS_FFTW3Q_THREADS) || !defined(FFTW_HAS_FFTW3Q_OMP)
+ %ignore fftwq_init_threads();
+ %ignore fftwq_plan_with_nthreads(int);
+ %ignore fftwq_cleanup_threads();
+ #endif
+#endif
+%rename(is_) fftw_iodim_do_not_use_me::is;
+%rename(is_) fftw_iodim64_do_not_use_me::is;
+%include "fftw3.h"
+
+
+/* fftw3 c++ wrappers */
+%include "fft/fftw3.h"
+%template(Fftw3f) hysop::fft::Fftw3<float>;
+%template(Fftw3d) hysop::fft::Fftw3<double>;
+%template(Fftw3l) hysop::fft::Fftw3<long double>;
+%template(Fftw3q) hysop::fft::Fftw3<__float128>;
+
+
+/* fft transforms */
+%ignore hysop::fft::operator<<;
+%include "fft/transform.h"
+%extend hysop::fft::Transform {
+ const char* __str__() const {
+ return $self->toString().c_str();
+ }
+}
+
+
+/* fft extensions */
+%ignore hysop::fft::operator<<;
+%include "fft/extension.h"
+
+
+/* fft domain configuration */
+%ignore hysop::fft::operator<<;
+%ignore hysop::fft::FftDomainConfiguration::operator[];
+%include "fft/fftDomainConfiguration.h"
+
+%define INSTANTIATE_FFT_DOMAIN_CONFIG_IMPL(CLASS_NAME,CLASS_TEMPLATE...)
+ %template(CLASS_NAME) CLASS_TEMPLATE;
+ %extend CLASS_TEMPLATE {
+ CLASS_TEMPLATE::ExtensionPair __get_item__(unsigned int k) const {
+ return $self->operator[](k);
+ }
+ const char* __str__() const {
+ std::stringstream ss;
+ ss << *($self);
+ return ss.str().c_str();
+ }
+ }
+%enddef
+
+%define INSTANTIATE_FFT_DOMAIN_CONFIG(INT)
+ INSTANTIATE_FFT_DOMAIN_CONFIG_IMPL(FftDomainConfiguration ## INT ## D, hysop::fft::FftDomainConfiguration<INT>)
+%enddef
+
+%formacro(INSTANTIATE_FFT_DOMAIN_CONFIG, INSTANTIATED_DIMENSIONS)
+
diff --git a/HySoP/swig/hysop++/hysop++.i b/HySoP/swig/hysop++/hysop++.i
new file mode 100644
index 000000000..a681efa69
--- /dev/null
+++ b/HySoP/swig/hysop++/hysop++.i
@@ -0,0 +1,5 @@
+
+/* hysop++ */
+%include "utils.i"
+%include "domain.i"
+%include "fft.i"
diff --git a/HySoP/swig/hysop++/utils.i b/HySoP/swig/hysop++/utils.i
new file mode 100644
index 000000000..541179ea9
--- /dev/null
+++ b/HySoP/swig/hysop++/utils.i
@@ -0,0 +1,29 @@
+
+%module utils
+
+// -*- C++ -*-
+%{
+ #include "utils/types.h"
+ #include "utils/constants.h"
+%}
+
+/* types */
+%warnfilter(342) hysop::multi_array;
+%warnfilter(342) hysop::multi_array_view;
+%warnfilter(342) hysop::const_multi_array_view;
+%warnfilter(342) hysop::multi_array_ref;
+%warnfilter(342) hysop::const_multi_array_ref;
+%ignore hysop::multi_array;
+%ignore hysop::multi_array_view;
+%ignore hysop::const_multi_array_view;
+%ignore hysop::multi_array_ref;
+%ignore hysop::const_multi_array_ref;
+%include "utils/types.h"
+
+%define INSTANTIATE_SHAPE(INT) %template(Shape ## INT ## D) hysop::Shape<INT>; %enddef
+%define INSTANTIATE_OFFSET(INT) %template(Offset ## INT ## D) hysop::Offset<INT>; %enddef
+%formacro(INSTANTIATE_SHAPE, INSTANTIATED_DIMENSIONS)
+%formacro(INSTANTIATE_OFFSET, INSTANTIATED_DIMENSIONS)
+
+/* constants */
+%include "utils/constants.h"
diff --git a/HySoP/swig/start.i b/HySoP/swig/start.i
index ee33d7837..24aa1a26a 100644
--- a/HySoP/swig/start.i
+++ b/HySoP/swig/start.i
@@ -8,3 +8,7 @@
%init %{
import_array();
%}
+
+
+//usefull macros used to instanciate templates
+%define INSTANTIATED_DIMENSIONS 1,2,3,4,5,6,7 %enddef
diff --git a/HySoP/swig/types.i b/HySoP/swig/types.i
deleted file mode 100644
index 6a2ef1049..000000000
--- a/HySoP/swig/types.i
+++ /dev/null
@@ -1,7 +0,0 @@
-
-// -*- C++ -*-
-%{
- #include "utils/types.h"
-%}
-
-%include "utils/types.h"
--
GitLab