Add pure xtensor fft implementation

spectre-ns · Apr 24, 2024 · 1d78377 · 1d78377
1 parent a17f3de
commit 1d78377
Show file tree

Hide file tree

Showing 10 changed files with 268 additions and 0 deletions.
diff --git a/.gitignore b/.gitignore
@@ -62,3 +62,4 @@ __pycache__
 
 # Generated files
 *.pc
+.vscode/settings.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -140,6 +140,7 @@ set(XTENSOR_HEADERS
     ${XTENSOR_INCLUDE_DIR}/xtensor/xfixed.hpp
     ${XTENSOR_INCLUDE_DIR}/xtensor/xfunction.hpp
     ${XTENSOR_INCLUDE_DIR}/xtensor/xfunctor_view.hpp
+    ${XTENSOR_INCLUDE_DIR}/xtensor/xfft.hpp
     ${XTENSOR_INCLUDE_DIR}/xtensor/xgenerator.hpp
     ${XTENSOR_INCLUDE_DIR}/xtensor/xhistogram.hpp
     ${XTENSOR_INCLUDE_DIR}/xtensor/xindex_view.hpp
@@ -199,6 +200,7 @@ target_link_libraries(xtensor INTERFACE xtl)
 
 OPTION(XTENSOR_ENABLE_ASSERT "xtensor bound check" OFF)
 OPTION(XTENSOR_CHECK_DIMENSION "xtensor dimension check" OFF)
+OPTION(XTENSOR_FORCE_TEMPORARY_MEMORY_IN_ASSIGNMENTS "xtensor force the use of temporary memory when assigning instead of an automatic overlap check" ON)
 OPTION(BUILD_TESTS "xtensor test suite" OFF)
 OPTION(BUILD_BENCHMARK "xtensor benchmark" OFF)
 OPTION(DOWNLOAD_GTEST "build gtest from downloaded sources" OFF)
@@ -219,6 +221,10 @@ if(XTENSOR_CHECK_DIMENSION)
     add_definitions(-DXTENSOR_ENABLE_CHECK_DIMENSION)
 endif()
 
+if(XTENSOR_FORCE_TEMPORARY_MEMORY_IN_ASSIGNMENTS)
+    add_definitions(-DXTENSOR_FORCE_TEMPORARY_MEMORY_IN_ASSIGNMENTS)
+endif()
+
 if(DEFAULT_COLUMN_MAJOR)
     add_definitions(-DXTENSOR_DEFAULT_LAYOUT=layout_type::column_major)
 endif()

diff --git a/docs/source/api/container_index.rst b/docs/source/api/container_index.rst
@@ -33,3 +33,4 @@ xexpression API is actually implemented in ``xstrided_container`` and ``xcontain
    xindex_view
    xfunctor_view
    xrepeat
+   xfft
diff --git a/include/xtensor/xbroadcast.hpp b/include/xtensor/xbroadcast.hpp
@@ -118,6 +118,29 @@ namespace xt
         return linear_end(c.expression());
     }
 
+    /*************************************
+     * overlapping_memory_checker_traits *
+     *************************************/
+
+    template <class E>
+    struct overlapping_memory_checker_traits<
+        E,
+        std::enable_if_t<!has_memory_address<E>::value && is_specialization_of<xbroadcast, E>::value>>
+    {
+        static bool check_overlap(const E& expr, const memory_range& dst_range)
+        {
+            if (expr.size() == 0)
+            {
+                return false;
+            }
+            else
+            {
+                using ChildE = std::decay_t<decltype(expr.expression())>;
+                return overlapping_memory_checker_traits<ChildE>::check_overlap(expr.expression(), dst_range);
+            }
+        }
+    };
+
     /**
      * @class xbroadcast
      * @brief Broadcasted xexpression to a specified shape.

diff --git a/include/xtensor/xfunction.hpp b/include/xtensor/xfunction.hpp
@@ -162,6 +162,42 @@ namespace xt
     {
     };
 
+    /*************************************
+     * overlapping_memory_checker_traits *
+     *************************************/
+
+    template <class E>
+    struct overlapping_memory_checker_traits<
+        E,
+        std::enable_if_t<!has_memory_address<E>::value && is_specialization_of<xfunction, E>::value>>
+    {
+        template <std::size_t I = 0, class... T, std::enable_if_t<(I == sizeof...(T)), int> = 0>
+        static bool check_tuple(const std::tuple<T...>&, const memory_range&)
+        {
+            return false;
+        }
+
+        template <std::size_t I = 0, class... T, std::enable_if_t<(I < sizeof...(T)), int> = 0>
+        static bool check_tuple(const std::tuple<T...>& t, const memory_range& dst_range)
+        {
+            using ChildE = std::decay_t<decltype(std::get<I>(t))>;
+            return overlapping_memory_checker_traits<ChildE>::check_overlap(std::get<I>(t), dst_range)
+                   || check_tuple<I + 1>(t, dst_range);
+        }
+
+        static bool check_overlap(const E& expr, const memory_range& dst_range)
+        {
+            if (expr.size() == 0)
+            {
+                return false;
+            }
+            else
+            {
+                return check_tuple(expr.arguments(), dst_range);
+            }
+        }
+    };
+
     /*************
      * xfunction *
      *************/

diff --git a/include/xtensor/xgenerator.hpp b/include/xtensor/xgenerator.hpp
@@ -76,6 +76,21 @@ namespace xt
         using size_type = std::size_t;
     };
 
+    /*************************************
+     * overlapping_memory_checker_traits *
+     *************************************/
+
+    template <class E>
+    struct overlapping_memory_checker_traits<
+        E,
+        std::enable_if_t<!has_memory_address<E>::value && is_specialization_of<xgenerator, E>::value>>
+    {
+        static bool check_overlap(const E&, const memory_range&)
+        {
+            return false;
+        }
+    };
+
     /**
      * @class xgenerator
      * @brief Multidimensional function operating on indices.

diff --git a/include/xtensor/xmath.hpp b/include/xtensor/xmath.hpp
@@ -338,6 +338,7 @@ namespace xt
         XTENSOR_UNARY_MATH_FUNCTOR(isfinite);
         XTENSOR_UNARY_MATH_FUNCTOR(isinf);
         XTENSOR_UNARY_MATH_FUNCTOR(isnan);
+        XTENSOR_UNARY_MATH_FUNCTOR(conj);
     }
 
 #undef XTENSOR_UNARY_MATH_FUNCTOR

diff --git a/include/xtensor/xsemantic.hpp b/include/xtensor/xsemantic.hpp
@@ -217,6 +217,29 @@ namespace xt
     template <class E, class R = void>
     using disable_xcontainer_semantics = typename std::enable_if<!has_container_semantics<E>::value, R>::type;
 
+
+    template <class D>
+    class xview_semantic;
+
+    template <class E>
+    struct overlapping_memory_checker_traits<
+        E,
+        std::enable_if_t<!has_memory_address<E>::value && is_crtp_base_of<xview_semantic, E>::value>>
+    {
+        static bool check_overlap(const E& expr, const memory_range& dst_range)
+        {
+            if (expr.size() == 0)
+            {
+                return false;
+            }
+            else
+            {
+                using ChildE = std::decay_t<decltype(expr.expression())>;
+                return overlapping_memory_checker_traits<ChildE>::check_overlap(expr.expression(), dst_range);
+            }
+        }
+    };
+
     /**
      * @class xview_semantic
      * @brief Implementation of the xsemantic_base interface for
@@ -598,8 +621,22 @@ namespace xt
     template <class E>
     inline auto xsemantic_base<D>::operator=(const xexpression<E>& e) -> derived_type&
     {
+#ifdef XTENSOR_FORCE_TEMPORARY_MEMORY_IN_ASSIGNMENTS
         temporary_type tmp(e);
         return this->derived_cast().assign_temporary(std::move(tmp));
+#else
+        auto&& this_derived = this->derived_cast();
+        auto memory_checker = make_overlapping_memory_checker(this_derived);
+        if (memory_checker.check_overlap(e.derived_cast()))
+        {
+            temporary_type tmp(e);
+            return this_derived.assign_temporary(std::move(tmp));
+        }
+        else
+        {
+            return this->assign(e);
+        }
+#endif
     }
 
     /**************************************

diff --git a/include/xtensor/xutils.hpp b/include/xtensor/xutils.hpp
@@ -119,6 +119,20 @@ namespace xt
         using type = T;
     };
 
+    /***************************************
+     * is_specialization_of implementation *
+     ***************************************/
+
+    template <template <class...> class TT, class T>
+    struct is_specialization_of : std::false_type
+    {
+    };
+
+    template <template <class...> class TT, class... Ts>
+    struct is_specialization_of<TT, TT<Ts...>> : std::true_type
+    {
+    };
+
     /*******************************
      * remove_class implementation *
      *******************************/
@@ -860,6 +874,139 @@ namespace xt
     {
     };
 
+    /*************************************
+     * overlapping_memory_checker_traits *
+     *************************************/
+
+    template <class T, class Enable = void>
+    struct has_memory_address : std::false_type
+    {
+    };
+
+    template <class T>
+    struct has_memory_address<T, void_t<decltype(std::addressof(*std::declval<T>().begin()))>> : std::true_type
+    {
+    };
+
+    struct memory_range
+    {
+        // Checking pointer overlap is more correct in integer values,
+        // for more explanation check https://devblogs.microsoft.com/oldnewthing/20170927-00/?p=97095
+        const uintptr_t m_first = 0;
+        const uintptr_t m_last = 0;
+
+        explicit memory_range() = default;
+
+        template <class T>
+        explicit memory_range(T* first, T* last)
+            : m_first(reinterpret_cast<uintptr_t>(last < first ? last : first))
+            , m_last(reinterpret_cast<uintptr_t>(last < first ? first : last))
+        {
+        }
+
+        template <class T>
+        bool overlaps(T* first, T* last) const
+        {
+            if (first <= last)
+            {
+                return reinterpret_cast<uintptr_t>(first) <= m_last
+                       && reinterpret_cast<uintptr_t>(last) >= m_first;
+            }
+            else
+            {
+                return reinterpret_cast<uintptr_t>(last) <= m_last
+                       && reinterpret_cast<uintptr_t>(first) >= m_first;
+            }
+        }
+    };
+
+    template <class E, class Enable = void>
+    struct overlapping_memory_checker_traits
+    {
+        static bool check_overlap(const E&, const memory_range&)
+        {
+            return true;
+        }
+    };
+
+    template <class E>
+    struct overlapping_memory_checker_traits<E, std::enable_if_t<has_memory_address<E>::value>>
+    {
+        static bool check_overlap(const E& expr, const memory_range& dst_range)
+        {
+            if (expr.size() == 0)
+            {
+                return false;
+            }
+            else
+            {
+                return dst_range.overlaps(std::addressof(*expr.begin()), std::addressof(*expr.rbegin()));
+            }
+        }
+    };
+
+    struct overlapping_memory_checker_base
+    {
+        memory_range m_dst_range;
+
+        explicit overlapping_memory_checker_base() = default;
+
+        explicit overlapping_memory_checker_base(memory_range dst_memory_range)
+            : m_dst_range(std::move(dst_memory_range))
+        {
+        }
+
+        template <class E>
+        bool check_overlap(const E& expr) const
+        {
+            if (!m_dst_range.m_first || !m_dst_range.m_last)
+            {
+                return false;
+            }
+            else
+            {
+                return overlapping_memory_checker_traits<E>::check_overlap(expr, m_dst_range);
+            }
+        }
+    };
+
+    template <class Dst, class Enable = void>
+    struct overlapping_memory_checker : overlapping_memory_checker_base
+    {
+        explicit overlapping_memory_checker(const Dst&)
+            : overlapping_memory_checker_base()
+        {
+        }
+    };
+
+    template <class Dst>
+    struct overlapping_memory_checker<Dst, std::enable_if_t<has_memory_address<Dst>::value>>
+        : overlapping_memory_checker_base
+    {
+        explicit overlapping_memory_checker(const Dst& aDst)
+            : overlapping_memory_checker_base(
+                [&]()
+                {
+                    if (aDst.size() == 0)
+                    {
+                        return memory_range();
+                    }
+                    else
+                    {
+                        return memory_range(std::addressof(*aDst.begin()), std::addressof(*aDst.rbegin()));
+                    }
+                }()
+            )
+        {
+        }
+    };
+
+    template <class Dst>
+    auto make_overlapping_memory_checker(const Dst& a_dst)
+    {
+        return overlapping_memory_checker<Dst>(a_dst);
+    }
+
     /********************
      * rebind_container *
      ********************/

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -186,6 +186,7 @@ set(XTENSOR_TESTS
     test_xdynamic_view.cpp
     test_xfunctor_adaptor.cpp
     test_xfixed.cpp
+    test_xfft.cpp
     test_xhistogram.cpp
     test_xpad.cpp
     test_xindex_view.cpp