espressomd · fweik · Apr 11, 2019 · Apr 8, 2019 · Apr 8, 2019 · Apr 9, 2019
diff --git a/src/core/electrostatics_magnetostatics/p3m_gpu_common.hpp b/src/core/electrostatics_magnetostatics/p3m_gpu_common.hpp
diff --git a/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu b/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu
@@ -48,7 +48,12 @@
 #include "EspressoSystemInterface.hpp"
 #include "global.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "p3m_gpu_common.hpp"
+
+#include <utils/math/int_pow.hpp>
+using Utils::int_pow;
+#include <utils/math/sinc.hpp>
+#include <utils/math/sqr.hpp>
+using Utils::sqr;
 
 #if defined(OMPI_MPI_H) || defined(_MPI_H)
 #error CU-file includes mpi.h! This should not happen!
@@ -258,13 +263,13 @@ __device__ void static Aliasing_sums_ik(const P3MGpuData p, int NX, int NY,
 
   for (MX = -P3M_BRILLOUIN; MX <= P3M_BRILLOUIN; MX++) {
     NMX = ((NX > p.mesh[0] / 2) ? NX - p.mesh[0] : NX) + p.mesh[0] * MX;
-    S1 = int_pow<2 * cao>(csinc(Meshi[0] * NMX));
+    S1 = int_pow<2 * cao>(Utils::sinc(Meshi[0] * NMX));
     for (MY = -P3M_BRILLOUIN; MY <= P3M_BRILLOUIN; MY++) {
       NMY = ((NY > p.mesh[1] / 2) ? NY - p.mesh[1] : NY) + p.mesh[1] * MY;
-      S2 = S1 * int_pow<2 * cao>(csinc(Meshi[1] * NMY));
+      S2 = S1 * int_pow<2 * cao>(Utils::sinc(Meshi[1] * NMY));
       for (MZ = -P3M_BRILLOUIN; MZ <= P3M_BRILLOUIN; MZ++) {
         NMZ = ((NZ > p.mesh[2] / 2) ? NZ - p.mesh[2] : NZ) + p.mesh[2] * MZ;
-        S3 = S2 * int_pow<2 * cao>(csinc(Meshi[2] * NMZ));
+        S3 = S2 * int_pow<2 * cao>(Utils::sinc(Meshi[2] * NMZ));
 
         NM2 = sqr(NMX * Leni[0]) + sqr(NMY * Leni[1]) + sqr(NMZ * Leni[2]);
         *Nenner += S3;

diff --git a/src/core/electrostatics_magnetostatics/p3m_gpu_error_cuda.cu b/src/core/electrostatics_magnetostatics/p3m_gpu_error_cuda.cu
@@ -31,7 +31,11 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #include "cuda_utils.hpp"
 #include "p3m_gpu_error.hpp"
 
-#include "p3m_gpu_common.hpp"
+#include "utils/math/int_pow.hpp"
+using Utils::int_pow;
+#include "utils/math/sinc.hpp"
+#include "utils/math/sqr.hpp"
+using Utils::sqr;
 
 #if defined(OMPI_MPI_H) || defined(_MPI_H)
 #error CU-file includes mpi.h! This should not happen!
@@ -95,8 +99,9 @@ __global__ void p3m_k_space_error_gpu_kernel_ik(int3 mesh, double3 meshi,
                       p3m_analytic_cotangent_sum<cao>(ny, meshi.y);
     const double ex = exp(-(PI * alpha_L_i) * (PI * alpha_L_i) * n2);
     const double ex2 = sqr(ex);
-    const double U2 = int_pow<2 * cao>(
-        csinc(meshi.x * nx) * csinc(meshi.y * ny) * csinc(meshi.z * nz));
+    const double U2 =
+        int_pow<2 * cao>(Utils::sinc(meshi.x * nx) * Utils::sinc(meshi.y * ny) *
+                         Utils::sinc(meshi.z * nz));
     auto const alias1 = ex2 / n2;
     auto const d = alias1 - sqr(U2 * ex / cs) / n2;
 
@@ -142,8 +147,8 @@ __global__ void p3m_k_space_error_gpu_kernel_ad(const int3 mesh,
 
           ex2 = sqr(ex);
 
-          U2 = pow((double)csinc(meshi.x * nmx) * csinc(meshi.y * nmy) *
-                       csinc(meshi.z * nmz),
+          U2 = pow((double)Utils::sinc(meshi.x * nmx) *
+                       Utils::sinc(meshi.y * nmy) * Utils::sinc(meshi.z * nmz),
                    2.0 * cao);
 
           alias1 += ex2 / n2;
@@ -201,8 +206,8 @@ __global__ void p3m_k_space_error_gpu_kernel_ik_i(const int3 mesh,
 
           ex2 = sqr(ex);
 
-          U2 = pow((double)csinc(meshi.x * nmx) * csinc(meshi.y * nmy) *
-                       csinc(meshi.z * nmz),
+          U2 = pow((double)Utils::sinc(meshi.x * nmx) *
+                       Utils::sinc(meshi.y * nmy) * Utils::sinc(meshi.z * nmz),
                    2.0 * cao);
 
           alias1 += ex2 / n2;
@@ -264,8 +269,8 @@ __global__ void p3m_k_space_error_gpu_kernel_ad_i(const int3 mesh,
 
           ex2 = sqr(ex);
 
-          U2 = pow((double)csinc(meshi.x * nmx) * csinc(meshi.y * nmy) *
-                       csinc(meshi.z * nmz),
+          U2 = pow((double)Utils::sinc(meshi.x * nmx) *
+                       Utils::sinc(meshi.y * nmy) * Utils::sinc(meshi.z * nmz),
                    2.0 * cao);
 
           alias1 += ex2 / n2;

diff --git a/src/core/utils/Span.hpp b/src/core/utils/Span.hpp
@@ -26,7 +26,6 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #include <type_traits>
 
 #include "device_qualifier.hpp"
-#include "type_traits.hpp"
 
 namespace Utils {
 namespace detail {
@@ -72,11 +71,8 @@ template <class T> class Span {
       typename std::enable_if<detail::has_data<T, U>::value, U>::type;
 
 public:
-  DEVICE_QUALIFIER
   Span() = default;
-  DEVICE_QUALIFIER
   Span(const Span &) = default;
-  DEVICE_QUALIFIER
   Span &operator=(const Span &) = default;
 
   DEVICE_QUALIFIER

diff --git a/src/core/utils/interpolation/detail/ll_and_dist.hpp b/src/core/utils/interpolation/detail/ll_and_dist.hpp
@@ -23,14 +23,12 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <array>
 #include <cmath>
+#include <type_traits>
 #include <utility>
 
-#include "utils/type_traits.hpp"
-
 namespace Utils {
 namespace Interpolation {
 namespace detail {
-using std::enable_if_t;
 
 struct Block {
   /* Index of the lower left corner of the assignment cube */
@@ -42,7 +40,7 @@ struct Block {
 template <size_t order, typename = void> struct ll_and_dist_;
 
 template <size_t order>
-struct ll_and_dist_<order, enable_if_t<(order % 2) == 1>> {
+struct ll_and_dist_<order, std::enable_if_t<(order % 2) == 1>> {
   Block operator()(const Vector3d &pos, const Vector3d &grid_spacing,
                    const Vector3d &offset) const {
     Vector3d dist;
@@ -61,7 +59,7 @@ struct ll_and_dist_<order, enable_if_t<(order % 2) == 1>> {
 };
 
 template <size_t order>
-struct ll_and_dist_<order, enable_if_t<(order % 2) == 0>> {
+struct ll_and_dist_<order, std::enable_if_t<(order % 2) == 0>> {
   Block operator()(const Vector3d &pos, const Vector3d &grid_spacing,
                    const Vector3d &offset) const {
     Vector3d dist;

diff --git a/src/core/utils/math/int_pow.hpp b/src/core/utils/math/int_pow.hpp
@@ -23,28 +23,32 @@
 #ifndef UTILS_MATH_INT_POW_HPP
 #define UTILS_MATH_INT_POW_HPP
 
-#include "utils/type_traits.hpp"
+#include "utils/device_qualifier.hpp"
+
+#include <type_traits>
 
 namespace Utils {
 namespace detail {
 template <class T, unsigned n, class = void> struct int_pow_impl {
-  constexpr T operator()(T x) const {
+  DEVICE_QUALIFIER constexpr T operator()(T x) const {
     return x * int_pow_impl<T, (n - 1) / 2>{}(x * x);
   }
 };
 
 /* Specializaton for n even */
 template <class T, unsigned n>
 struct int_pow_impl<T, n, std::enable_if_t<n % 2 == 0>> {
-  constexpr T operator()(T x) const { return int_pow_impl<T, n / 2>{}(x * x); }
+  DEVICE_QUALIFIER constexpr T operator()(T x) const {
+    return int_pow_impl<T, n / 2>{}(x * x);
+  }
 };
 
 template <class T> struct int_pow_impl<T, 1> {
-  constexpr T operator()(T x) const { return x; }
+  DEVICE_QUALIFIER constexpr T operator()(T x) const { return x; }
 };
 
 template <class T> struct int_pow_impl<T, 0> {
-  constexpr T operator()(T x) const { return T{1}; }
+  DEVICE_QUALIFIER constexpr T operator()(T x) const { return T{1}; }
 };
 } // namespace detail
 
@@ -55,7 +59,7 @@ template <class T> struct int_pow_impl<T, 0> {
  * at compile time. It uses exponentiation by
  * squaring to construct a efficient function.
  */
-template <unsigned n, typename T> constexpr T int_pow(T x) {
+template <unsigned n, typename T> DEVICE_QUALIFIER constexpr T int_pow(T x) {
   return detail::int_pow_impl<T, n>{}(x);
 }
 } // namespace Utils

diff --git a/src/core/utils/math/sinc.hpp b/src/core/utils/math/sinc.hpp
@@ -23,6 +23,9 @@
 #ifndef UTILS_MATH_SINC_HPP
 #define UTILS_MATH_SINC_HPP
 
+#include "utils/constants.hpp"
+#include "utils/device_qualifier.hpp"
+
 #include <boost/math/constants/constants.hpp>
 
 #include <cmath>
@@ -41,10 +44,10 @@ namespace Utils {
  * also save time, since it reduces the number of function calls to
  * sin().
  */
-template <typename T> T sinc(T d) {
+template <typename T> DEVICE_QUALIFIER T sinc(T d) {
   const constexpr T epsi = 0.1;
 
-  const auto PId = boost::math::constants::pi<T>() * d;
+  const auto PId = PI * d;
 
   if (std::abs(d) > epsi)
     return std::sin(PId) / PId;

diff --git a/src/core/utils/math/sqr.hpp b/src/core/utils/math/sqr.hpp
@@ -19,10 +19,11 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #ifndef UTILS_MATH_SQR_HPP
 #define UTILS_MATH_SQR_HPP
 
+#include "utils/device_qualifier.hpp"
+
 namespace Utils {
 /** Calculates the SQuaRe of x */
-template <typename T> constexpr T sqr(T x) { return x * x; }
-
+template <typename T> DEVICE_QUALIFIER constexpr T sqr(T x) { return x * x; }
 } // namespace Utils
 
 #endif
diff --git a/src/core/utils/type_traits.hpp b/src/core/utils/type_traits.hpp
@@ -38,18 +38,6 @@ struct function_remove_const<R(Args...) const> {
   using type = R(Args...);
 };
 
-template <class T>
-using function_remove_const_t = typename function_remove_const<T>::type;
-
-/**
- * @brief True iff T is an instantiation of Template.
- */
-template <typename T, template <typename...> class Template>
-struct is_instance_of : public std::false_type {};
-
-template <typename... T, template <typename...> class Template>
-struct is_instance_of<Template<T...>, Template> : public std::true_type {};
-
 template <class...> struct conjunction : std::true_type {};
 template <class B1> struct conjunction<B1> : B1 {};
 template <class B1, class... Bn>