Reduction

cmake/public/gridtools_setup_targets.cmake

@@ -238,6 +238,9 @@ macro(_gt_setup_targets _config_mode clang_cuda_mode)
     _gt_add_library(${_config_mode} stencil_naive)
     target_link_libraries(${_gt_namespace}stencil_naive INTERFACE ${_gt_namespace}gridtools)
 
+    _gt_add_library(${_config_mode} reduction_naive)
+    target_link_libraries(${_gt_namespace}reduction_naive INTERFACE ${_gt_namespace}gridtools)
+
     set(_required_nlohmann_json_version "3.7.3")
     include(get_nlohmann_json)
     get_nlohmann_json(${_required_nlohmann_json_version})
@@ -248,6 +251,7 @@ macro(_gt_setup_targets _config_mode clang_cuda_mode)
     endif()
 
     set(GT_STENCILS naive)
+    set(GT_REDUCTIONS naive)
     set(GT_STORAGES cpu_kfirst cpu_ifirst)
     set(GT_GCL_ARCHS)
 
@@ -260,6 +264,10 @@ macro(_gt_setup_targets _config_mode clang_cuda_mode)
         target_link_libraries(${_gt_namespace}stencil_gpu_horizontal INTERFACE ${_gt_namespace}gridtools _gridtools_cuda)
         list(APPEND GT_STENCILS gpu_horizontal)
 
+        _gt_add_library(${_config_mode} reduction_gpu)
+        target_link_libraries(${_gt_namespace}reduction_gpu INTERFACE ${_gt_namespace}gridtools _gridtools_cuda)
+        list(APPEND GT_REDUCTIONS gpu)
+
         if(MPI_CXX_FOUND)
             option(GT_GCL_GPU "Disable if your MPI implementation is not CUDA-aware" ON)
             if(GT_GCL_GPU)

include/gridtools/common/ct_dispatch.hpp

@@ -0,0 +1,200 @@
+/*
+ * GridTools
+ *
+ * Copyright (c) 2014-2019, ETH Zurich
+ * All rights reserved.
+ *
+ * Please, refer to the LICENSE file in the root directory.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+
+/***
+ *  `ct_dispatch` performs compile time dispatch on runtime value
+ *
+ *  Usage:
+ *
+ *  Say you have a function that accepts an integer in compile time:
+ *
+ *  template <size_t N> int foo() { ... }
+ *
+ *  And you need something like:
+ *
+ *  auto bar(int n) {
+ *      switch(n) {
+ *      case 0:
+ *          return foo<0>();
+ *      case 1:
+ *          return foo<1>();
+ *      case 2:
+ *          return foo<2>();
+ *      case 3:
+ *          return foo<3>();
+ *      }
+ *  }
+ *
+ *  You can use `ct_dispatch` here to reduce the boilerplate:
+ *
+ *  auto bar(int n) {
+ *      return ct_dispatch<4>([](auto n) {
+ *          return foo<decltype(n)::value>();)
+ *      }, n);
+ *  }
+ *
+ */
+
+#include <cassert>
+#include <cstdlib>
+#include <type_traits>
+#include <utility>
+
+namespace gridtools {
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 1, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n == 0);
+        return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 2, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 2);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 3, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 3);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 4, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 4);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        case 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 3>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 5, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 5);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        case 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        case 3:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 3>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 4>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 6, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 6);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        case 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        case 3:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 3>());
+        case 4:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 4>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 5>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 7, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 7);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        case 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        case 3:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 3>());
+        case 4:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 4>());
+        case 5:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 5>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 6>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<Lim == 8, int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < 8);
+        switch (n) {
+        case 0:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 0>());
+        case 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 1>());
+        case 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 2>());
+        case 3:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 3>());
+        case 4:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 4>());
+        case 5:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 5>());
+        case 6:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 6>());
+        default:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, 7>());
+        }
+    }
+
+    template <size_t Lim, class Sink, std::enable_if_t<(Lim > 8), int> = 0>
+    auto ct_dispatch(Sink &&sink, size_t n) {
+        assert(n < Lim);
+        switch (n) {
+        case Lim - 8:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 8>());
+        case Lim - 7:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 7>());
+        case Lim - 6:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 6>());
+        case Lim - 5:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 5>());
+        case Lim - 4:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 4>());
+        case Lim - 3:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 3>());
+        case Lim - 2:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 2>());
+        case Lim - 1:
+            return std::forward<Sink>(sink)(std::integral_constant<size_t, Lim - 1>());
+        default:
+            return ct_dispatch<Lim - 8>(std::forward<Sink>(sink), n);
+        }
+    }
+} // namespace gridtools

include/gridtools/common/hugepage_alloc.hpp

@@ -177,7 +177,7 @@ namespace gridtools {
                 return hugepage_mode::disabled;
             if (std::strcmp(env_value, "explicit") == 0)
                 return hugepage_mode::explicit_allocation;
-            std::fprintf(stderr, "warning: env variable GT_HUGEPAGE_MODE set to invalid value '%s'\n", env_value);
+            fprintf(stderr, "warning: env variable GT_HUGEPAGE_MODE set to invalid value '%s'\n", env_value);
             return hugepage_mode::transparent;
         }
 

include/gridtools/common/numeric.hpp

@@ -0,0 +1,54 @@
+/*
+ * GridTools
+ *
+ * Copyright (c) 2014-2019, ETH Zurich
+ * All rights reserved.
+ *
+ * Please, refer to the LICENSE file in the root directory.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+
+#pragma once
+
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <numeric>
+#include <type_traits>
+
+namespace gridtools {
+#if __cplusplus < 201703
+    template <class T>
+    constexpr T gcd(T m, T n) {
+        static_assert(!std::is_signed<T>::value, "");
+        return n == 0 ? m : gcd<T>(n, m % n);
+    }
+
+    template <class T, class U>
+    constexpr std::common_type_t<T, U> gcd(T m, U n) {
+        static_assert(std::is_integral<T>() && std::is_integral<U>(), "Arguments to gcd must be integer types");
+        static_assert(!std::is_same<std::remove_cv_t<T>, bool>(), "First argument to gcd cannot be bool");
+        static_assert(!std::is_same<std::remove_cv_t<U>, bool>(), "Second argument to gcd cannot be bool");
+        using res_t = std::common_type_t<T, U>;
+        using ures_t = std::make_unsigned_t<res_t>;
+        return static_cast<res_t>(gcd(static_cast<ures_t>(std::abs(m)), static_cast<ures_t>(std::abs(n))));
+    }
+
+    template <class T, class U>
+    constexpr std::common_type_t<T, U> lcm(T m, U n) {
+        static_assert(std::is_integral<T>() && std::is_integral<U>(), "Arguments to lcm must be integer types");
+        static_assert(!std::is_same<std::remove_cv_t<T>, bool>(), "First argument to gcd cannot be bool");
+        static_assert(!std::is_same<std::remove_cv_t<U>, bool>(), "Second argument to gcd cannot be bool");
+        if (m == 0 || n == 0)
+            return 0;
+        using res_t = std::common_type_t<T, U>;
+        res_t a = std::abs(m) / gcd(m, n);
+        res_t b = std::abs(n);
+        assert(std::numeric_limits<res_t>::max() / a > b);
+        return a * b;
+    }
+#else
+    using std::gcd;
+    using std::lcm;
+#endif
+} // namespace gridtools

include/gridtools/reduction.hpp

@@ -0,0 +1,13 @@
+/*
+ * GridTools
+ *
+ * Copyright (c) 2014-2019, ETH Zurich
+ * All rights reserved.
+ *
+ * Please, refer to the LICENSE file in the root directory.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+#pragma once
+
+#include "reduction/frontend.hpp"
+#include "reduction/functions.hpp"

include/gridtools/reduction/frontend.hpp

@@ -0,0 +1,89 @@
+/*
+ * GridTools
+ *
+ * Copyright (c) 2014-2019, ETH Zurich
+ * All rights reserved.
+ *
+ * Please, refer to the LICENSE file in the root directory.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+#pragma once
+
+#include <cassert>
+#include <cstdlib>
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+#include "../common/tuple.hpp"
+#include "../common/tuple_util.hpp"
+#include "../meta.hpp"
+#include "../sid/allocator.hpp"
+#include "../storage/traits.hpp"
+
+namespace gridtools {
+    namespace reduction {
+        namespace frontend_impl_ {
+            template <class Sizes>
+            inline auto zeros(Sizes const &sizes) {
+                return tuple_util::transform([](auto &&) { return integral_constant<int_t, 0>(); }, sizes);
+            }
+
+            template <class Sizes>
+            using zeros_type = decltype(zeros(std::declval<Sizes const &>()));
+
+            template <class Backend, class Origin, class Strides, class StridesKind, class Sizes>
+            struct reducible {
+                std::shared_ptr<void> m_alloc;
+                Origin m_origin;
+                size_t m_size;
+                Strides m_strides;
+                Sizes m_sizes;
+
+                template <class F>
+                auto reduce(F f) const {
+                    assert(m_size);
+                    return reduction_reduce(Backend(), f, m_origin(), m_size);
+                }
+
+                friend Strides sid_get_strides(reducible const &obj) { return obj.m_strides; }
+                friend Origin sid_get_origin(reducible const &obj) { return {obj.m_origin}; }
+                friend zeros_type<Sizes> sid_get_lower_bounds(reducible const &obj) { return zeros(obj.m_sizes); }
+                friend Sizes sid_get_upper_bounds(reducible const &obj) { return obj.m_sizes; }
+            };
+
+            template <class Backend, class Origin, class Strides, class StridesKind, class Sizes>
+            StridesKind sid_get_strides_kind(reducible<Backend, Origin, Strides, StridesKind, Sizes> const &);
+
+            template <class Backend, class StorageTraits, class Id = void, class T, class... Dims>
+            auto make_reducible(T const &initial_value, Dims... dims) {
+                auto alloc = sid::host_device::make_cached_allocator(
+                    [](size_t size) { return storage::traits::allocate<StorageTraits, char>(size); });
+                auto lengths = tuple_util::make<tuple>(dims...);
+                auto info = storage::traits::make_info<StorageTraits, T>(lengths);
+                auto strides = info.native_strides();
+                size_t data_size = info.length();
+                size_t rounded_size = reduction_round_size(Backend(), data_size);
+                size_t allocation_size = reduction_allocation_size(Backend(), rounded_size);
+                auto origin = allocate(alloc, meta::lazy::id<T>(), allocation_size);
+                reduction_fill(Backend(),
+                    initial_value,
+                    origin(),
+                    data_size,
+                    rounded_size,
+                    allocation_size,
+                    storage::traits::has_holes<StorageTraits, T>(lengths));
+                return reducible<Backend,
+                    decltype(origin),
+                    decltype(strides),
+                    storage::traits::strides_kind<StorageTraits, T, decltype(lengths), Id>,
+                    decltype(lengths)>{std::make_shared<decltype(alloc)>(std::move(alloc)),
+                    std::move(origin),
+                    rounded_size,
+                    std::move(strides),
+                    std::move(lengths)};
+            }
+        } // namespace frontend_impl_
+        using frontend_impl_::make_reducible;
+    } // namespace reduction
+} // namespace gridtools