simd/vc.hh

Go to the documentation of this file.

// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_COMMON_SIMD_VC_HH
#define DUNE_COMMON_SIMD_VC_HH
 
#include <Vc/Vc>
 
#include <cstddef>
#include <type_traits>
#include <utility>
 
#include <dune/common/indices.hh>
#include <dune/common/simd/base.hh>
#include <dune/common/simd/defaults.hh> // for anyFalse()
#include <dune/common/simd/loop.hh>
#include <dune/common/typetraits.hh>
 
namespace Dune {
  namespace Simd {
 
    namespace VcImpl {
 
      template<class V, class SFINAE = void>
      struct IsMask : std::false_type {};
 
      template<typename T, typename A>
      struct IsMask<Vc::Mask<T, A> > : std::true_type {};
 
      template<typename T, std::size_t n, typename V, std::size_t m>
      struct IsMask<Vc::SimdMaskArray<T, n, V, m> > : std::true_type {};
 
      template<class V, class SFINAE = void>
      struct IsVector : IsMask<V> {};
 
      template<typename T, typename A>
      struct IsVector<Vc::Vector<T, A> > : std::true_type {};
 
      template<typename T, std::size_t n, typename V, std::size_t m>
      struct IsVector<Vc::SimdArray<T, n, V, m> > : std::true_type {};
 
      template<typename T> struct IsVectorizable : std::false_type {};
      template<> struct IsVectorizable<double>        : std::true_type {};
      template<> struct IsVectorizable<float>         : std::true_type {};
      template<> struct IsVectorizable<std::int32_t>  : std::true_type {};
      template<> struct IsVectorizable<std::uint32_t> : std::true_type {};
      template<> struct IsVectorizable<std::int16_t>  : std::true_type {};
      template<> struct IsVectorizable<std::uint16_t> : std::true_type {};
 
 
      template<class V>
      class Proxy
      {
        static_assert(std::is_same<V, std::decay_t<V> >::value, "Class Proxy "
                      "may only be instantiated with unqualified types");
      public:
        using value_type = typename V::value_type;
 
      private:
        static_assert(std::is_arithmetic<value_type>::value,
                      "Only arithmetic types are supported");
        V &vec_;
        std::size_t idx_;
 
      public:
        Proxy(std::size_t idx, V &vec)
          : vec_(vec), idx_(idx)
        { }
 
        Proxy(const Proxy&) = delete;
        // allow move construction so we can return proxies from functions
        Proxy(Proxy&&) = default;
 
        operator value_type() const { return vec_[idx_]; }
 
        // assignment operators
#define DUNE_SIMD_VC_ASSIGNMENT(OP)                              \
        template<class T,                                        \
                 class = decltype(std::declval<value_type&>() OP \
                                  autoCopy(std::declval<T>()) )> \
        Proxy operator OP(T &&o) &&                              \
        {                                                        \
          vec_[idx_] OP autoCopy(std::forward<T>(o));            \
          return { idx_, vec_ };                                 \
        }
        DUNE_SIMD_VC_ASSIGNMENT(=);
        DUNE_SIMD_VC_ASSIGNMENT(*=);
        DUNE_SIMD_VC_ASSIGNMENT(/=);
        DUNE_SIMD_VC_ASSIGNMENT(%=);
        DUNE_SIMD_VC_ASSIGNMENT(+=);
        DUNE_SIMD_VC_ASSIGNMENT(-=);
        DUNE_SIMD_VC_ASSIGNMENT(<<=);
        DUNE_SIMD_VC_ASSIGNMENT(>>=);
        DUNE_SIMD_VC_ASSIGNMENT(&=);
        DUNE_SIMD_VC_ASSIGNMENT(^=);
        DUNE_SIMD_VC_ASSIGNMENT(|=);
#undef DUNE_SIMD_VC_ASSIGNMENT
 
        // unary (prefix) operators
        template<class T = value_type,
                 class = std::enable_if_t<!std::is_same<T, bool>::value> >
        Proxy operator++() { ++(vec_[idx_]); return *this; }
        template<class T = value_type,
                 class = std::enable_if_t<!std::is_same<T, bool>::value> >
        Proxy operator--() { --(vec_[idx_]); return *this; }
 
        // postfix operators
        template<class T = value_type,
                 class = std::enable_if_t<!std::is_same<T, bool>::value> >
        value_type operator++(int) { return vec_[idx_]++; }
        template<class T = value_type,
                 class = std::enable_if_t<!std::is_same<T, bool>::value> >
        value_type operator--(int) { return vec_[idx_]--; }
 
 
        // swap on proxies swaps the proxied vector entries.  As such, it
        // applies to rvalues of proxies too, not just lvalues
        friend void swap(const Proxy &a, const Proxy &b) {
          // don't use swap() ourselves -- not supported by Vc 1.3.0 (but is
          // supported by Vc 1.3.2)
          value_type tmp = std::move(a.vec_[a.idx_]);
          a.vec_[a.idx_] = std::move(b.vec_[b.idx_]);
          b.vec_[b.idx_] = std::move(tmp);
        }
        friend void swap(value_type &a, const Proxy &b) {
          // don't use swap() ourselves -- not supported by Vc 1.3.0 (but is
          // supported by Vc 1.3.2)
          value_type tmp = std::move(a);
          a = std::move(b.vec_[b.idx_]);
          b.vec_[b.idx_] = std::move(tmp);
        }
        friend void swap(const Proxy &a, value_type &b) {
          // don't use swap() ourselves -- not supported by Vc 1.3.0 (but is
          // supported by Vc 1.3.2)
          value_type tmp = std::move(a.vec_[a.idx_]);
          a.vec_[a.idx_] = std::move(b);
          b = std::move(tmp);
        }
 
        // binary operators
        //
        // Normally, these are provided by the conversion operator in
        // combination with C++'s builtin binary operators.  Other classes
        // that need to provide the binary operators themselves should either
        // 1. deduce the "foreign" operand type independently, i.e. use
        //      template<class... Args, class Foreign>
        //      auto operator@(MyClass<Args...>, Foreign);
        //    or
        // 2. not deduce anything from the foreign argument, i.e.
        //      template<class... Args>
        //      auto operator@(MyClass<Args...>,
        //                     typename MyClass<Args...>::value_type);
        //    or
        //      template<class T, class... Args>
        //      struct MyClass {
        //        auto operator@(T);
        //      }
        //    or
        //      template<class T, class... Args>
        //      struct MyClass {
        //        friend auto operator@(MyClass, T);
        //      }
        //
        // This allows either for an exact match (in the case of option 1.) or
        // for conversions to be applied to the foreign argument (options 2.).
        // In contrast, allowing some of the template parameters being deduced
        // from the self argument also being deduced from the foreign argument
        // will likely lead to ambiguous deduction when the foreign argument is
        // a proxy:
        //   template<class T, class... Args>
        //   auto operator@(MyClass<T, Args...>, T);
        // One class that suffers from this problem is std::complex.
        //
        // Note that option 1. is a bit dangerous, as the foreign argument is
        // catch-all.  This seems tempting in the case of a proxy class, as
        // the operator could just be forwarded to the proxied object with the
        // foreign argument unchanged, immediately creating interoperability
        // with arbitrary foreign classes.  However, if the foreign class also
        // choses option 1., this will result in ambiguous overloads, and there
        // is no clear guide to decide which class should provide the overload
        // and which should not.
        //
        // Fortunately, deferring to the conversion and the built-in operators
        // mostly works in the case of this proxy class, because only built-in
        // types can be proxied anyway.  Unfortunately, the Vc vectors and
        // arrays suffer from a slightly different problem.  They chose option
        // 1., but they can't just accept the argument type they are given,
        // since they need to somehow implement the operation in terms of
        // intrinsics.  So they check the argument whether it is one of the
        // expected types, and remove the operator from the overload set if it
        // isn't via SFINAE.  Of course, this proxy class is not one of the
        // expected types, even though it would convert to them...
        //
        // So what we have to do here, unfortunately, is to provide operators
        // for the Vc types explicitly, and hope that there won't be some Vc
        // version that gets the operators right, thus creating ambiguous
        // overloads.  Well, if guess it will be #ifdef time if it comes to
        // that.
#define DUNE_SIMD_VC_BINARY(OP)                                         \
        template<class T, class Abi>                                    \
        friend auto operator OP(const Vc::Vector<T, Abi> &l, Proxy&& r) \
          -> decltype(l OP std::declval<value_type>())                  \
        {                                                               \
          return l OP value_type(r);                                    \
        }                                                               \
        template<class T, class Abi>                                    \
        auto operator OP(const Vc::Vector<T, Abi> &r) &&                \
          -> decltype(std::declval<value_type>() OP r)                  \
        {                                                               \
          return value_type(*this) OP r;                                \
        }                                                               \
        template<class T, std::size_t n, class Vec, std::size_t m>      \
        friend auto                                                     \
        operator OP(const Vc::SimdArray<T, n, Vec, m> &l, Proxy&& r)    \
          -> decltype(l OP std::declval<value_type>())                  \
        {                                                               \
          return l OP value_type(r);                                    \
        }                                                               \
        template<class T, std::size_t n, class Vec, std::size_t m>      \
        auto operator OP(const Vc::SimdArray<T, n, Vec, m> &r) &&       \
          -> decltype(std::declval<value_type>() OP r)                  \
        {                                                               \
          return value_type(*this) OP r;                                \
        }
 
        DUNE_SIMD_VC_BINARY(*);
        DUNE_SIMD_VC_BINARY(/);
        DUNE_SIMD_VC_BINARY(%);
        DUNE_SIMD_VC_BINARY(+);
        DUNE_SIMD_VC_BINARY(-);
        DUNE_SIMD_VC_BINARY(<<);
        DUNE_SIMD_VC_BINARY(>>);
        DUNE_SIMD_VC_BINARY(&);
        DUNE_SIMD_VC_BINARY(^);
        DUNE_SIMD_VC_BINARY(|);
        DUNE_SIMD_VC_BINARY(<);
        DUNE_SIMD_VC_BINARY(>);
        DUNE_SIMD_VC_BINARY(<=);
        DUNE_SIMD_VC_BINARY(>=);
        DUNE_SIMD_VC_BINARY(==);
        DUNE_SIMD_VC_BINARY(!=);
#undef DUNE_SIMD_VC_BINARY
 
        // this is needed to implement broadcast construction from proxy as
        // the unadorned assignment operator cannot be a non-member
        template<class T, class Abi,
                 class = std::enable_if_t<std::is_convertible<value_type,
                                                              T>::value> >
        operator Vc::Vector<T, Abi>() &&
        {
          return value_type(*this);
        }
        template<class T, std::size_t n, class Vec, std::size_t m,
                 class = std::enable_if_t<std::is_convertible<value_type,
                                                              T>::value> >
        operator Vc::SimdArray<T, n, Vec, m>() &&
        {
          return value_type(*this);
        }
 
#define DUNE_SIMD_VC_ASSIGN(OP)                                         \
        template<class T, class Abi>                                    \
        friend auto operator OP(Vc::Vector<T, Abi> &l, Proxy&& r)       \
          -> decltype(l OP std::declval<value_type>())                  \
        {                                                               \
          return l OP value_type(r);                                    \
        }
 
        DUNE_SIMD_VC_ASSIGN(*=);
        DUNE_SIMD_VC_ASSIGN(/=);
        DUNE_SIMD_VC_ASSIGN(%=);
        DUNE_SIMD_VC_ASSIGN(+=);
        DUNE_SIMD_VC_ASSIGN(-=);
        DUNE_SIMD_VC_ASSIGN(&=);
        DUNE_SIMD_VC_ASSIGN(^=);
        DUNE_SIMD_VC_ASSIGN(|=);
        // The shift assignment would not be needed for Vc::Vector since it
        // has overloads for `int` rhs and the proxy can convert to that --
        // except that there is also overloads for Vector, and because of the
        // conversion operator needed to support unadorned assignments, the
        // proxy can convert to that, too.
        DUNE_SIMD_VC_ASSIGN(<<=);
        DUNE_SIMD_VC_ASSIGN(>>=);
#undef DUNE_SIMD_VC_ASSIGN
      };
 
    } // namespace VcImpl
 
    namespace Overloads {
 
 
      template<class V>
      struct ScalarType<V, std::enable_if_t<VcImpl::IsVector<V>::value> >
      {
        using type = typename V::value_type;
      };
 
 
      template<class V>
      struct RebindType<Simd::Scalar<V>, V,
                        std::enable_if_t<VcImpl::IsVector<V>::value> >
      {
        using type = V;
      };
 
 
      template<class V>
      struct RebindType<bool, V, std::enable_if_t<VcImpl::IsVector<V>::value &&
                                                  !VcImpl::IsMask<V>::value>>
      {
        using type = typename V::mask_type;
      };
 
 
      template<class M>
      struct RebindType<Scalar<typename M::Vector>, M,
                        std::enable_if_t<VcImpl::IsMask<M>::value>>
      {
        using type = typename M::Vector;
      };
 
 
      template<class S, class M>
      struct RebindType<S, M,
                        std::enable_if_t<
                          VcImpl::IsMask<M>::value &&
                          VcImpl::IsVectorizable<S>::value &&
                          !std::is_same<S, Scalar<typename M::Vector> >::value
                          > >
      {
        using type = Vc::fixed_size_simd<S, int(Simd::lanes<M>())>;
      };
 
 
      template<class S, class V>
      struct RebindType<S, V,
                        std::enable_if_t<VcImpl::IsVector<V>::value &&
                                         !VcImpl::IsMask<V>::value &&
                                         VcImpl::IsVectorizable<S>::value &&
                                         !std::is_same<S, Scalar<V> >::value> >
      {
        using type = Vc::fixed_size_simd<S, int(Simd::lanes<V>())>;
      };
 
 
      template<class S, class V>
      struct RebindType<S, V,
                        std::enable_if_t<VcImpl::IsVector<V>::value &&
                                         !VcImpl::IsVectorizable<S>::value &&
                                         !std::is_same<S, bool>::value &&
                                         !std::is_same<S, Scalar<V> >::value> >
      {
        using type = LoopSIMD<S, Simd::lanes<V>()>;
      };
 
 
      template<class V>
      struct LaneCount<V, std::enable_if_t<VcImpl::IsVector<V>::value> >
        : public index_constant<V::size()>
      { };
 
      template<class V>
      VcImpl::Proxy<V> lane(ADLTag<5, VcImpl::IsVector<V>::value>,
                            std::size_t l, V &v)
      {
        return { l, v };
      }
 
      template<class V>
      Scalar<V> lane(ADLTag<5, VcImpl::IsVector<V>::value>,
                     std::size_t l, const V &v)
      {
        return v[l];
      }
 
      /*
       * The hack with the SFINAE is necessary, because if I use just
       * Scalar<V> as the return type, the compiler still errors out if V is
       * an lvalue-reference T&.  You'd think he'd notice that he can't
       * instantiate this declaration for this template parameter, and would
       * simply remove it from the overload set, but no...
       */
      template<class V,
               class = std::enable_if_t<!std::is_reference<V>::value> >
      Scalar<V> lane(ADLTag<5, VcImpl::IsVector<V>::value>,
                     std::size_t l, V &&v)
      {
        return std::forward<V>(v)[l];
      }
 
      template<class V>
      V cond(ADLTag<5, VcImpl::IsVector<V>::value &&
                       !VcImpl::IsMask<V>::value>,
             const Mask<V> &mask, const V &ifTrue, const V &ifFalse)
      {
        return Vc::iif(mask, ifTrue, ifFalse);
      }
 
      /*
       * Kludge because iif seems to be unimplemented for masks
       */
      template<class V>
      V cond(ADLTag<5, VcImpl::IsMask<V>::value>,
             const V &mask, const V &ifTrue, const V &ifFalse)
      {
        return (mask && ifTrue) || (!mask && ifFalse);
      }
 
      template<class V>
      auto max(ADLTag<5, VcImpl::IsVector<V>::value &&
                         !VcImpl::IsMask<V>::value>,
               const V &v1, const V &v2)
      {
        return Simd::cond(v1 < v2, v2, v1);
      }
 
      template<class M>
      M max(ADLTag<5, VcImpl::IsMask<M>::value>,
            const M &m1, const M &m2)
      {
        // for SimdMaskArray, Vc returns the derived fixed_size_simd_mask
        return m1 || m2;
      }
 
      template<class V>
      auto min(ADLTag<5, VcImpl::IsVector<V>::value &&
                         !VcImpl::IsMask<V>::value>,
               const V &v1, const V &v2)
      {
        return Simd::cond(v1 < v2, v1, v2);
      }
 
      template<class M>
      M min(ADLTag<5, VcImpl::IsMask<M>::value>,
            const M &m1, const M &m2)
      {
        // for SimdMaskArray, Vc returns the derived fixed_size_simd_mask
        return m1 && m2;
      }
 
      template<class M>
      bool anyTrue (ADLTag<5, VcImpl::IsMask<M>::value>, const M &mask)
      {
        return Vc::any_of(mask);
      }
 
      template<class M>
      bool allTrue (ADLTag<5, VcImpl::IsMask<M>::value>, const M &mask)
      {
        return Vc::all_of(mask);
      }
 
      // nothing like anyFalse() in Vc, so let defaults.hh handle it
 
      template<class M>
      bool allFalse(ADLTag<5, VcImpl::IsMask<M>::value>, const M &mask)
      {
        return Vc::none_of(mask);
      }
 
      template<class V>
      auto max(ADLTag<5, VcImpl::IsVector<V>::value &&
                         !VcImpl::IsMask<V>::value>,
               const V &v)
      {
        return v.max();
      }
 
      template<class M>
      bool max(ADLTag<5, VcImpl::IsMask<M>::value>, const M &mask)
      {
        return Vc::any_of(mask);
      }
 
      template<class V>
      auto min(ADLTag<5, VcImpl::IsVector<V>::value &&
                         !VcImpl::IsMask<V>::value>,
               const V &v)
      {
        return v.min();
      }
 
      template<class M>
      bool min(ADLTag<5, VcImpl::IsMask<M>::value>, const M &mask)
      {
        return !Vc::any_of(!mask);
      }
 
      template<class S1, class V2>
      auto maskAnd(ADLTag<5, std::is_same<Mask<S1>, bool>::value &&
                             VcImpl::IsVector<V2>::value>,
                   const S1 &s1, const V2 &v2)
      {
        return Simd::Mask<V2>(Simd::mask(s1)) && Simd::mask(v2);
      }
 
      template<class V1, class S2>
      auto maskAnd(ADLTag<5, VcImpl::IsVector<V1>::value &&
                             std::is_same<Mask<S2>, bool>::value>,
                   const V1 &v1, const S2 &s2)
      {
        return Simd::mask(v1) && Simd::Mask<V1>(Simd::mask(s2));
      }
 
      template<class S1, class V2>
      auto maskOr(ADLTag<5, std::is_same<Mask<S1>, bool>::value &&
                            VcImpl::IsVector<V2>::value>,
                   const S1 &s1, const V2 &v2)
      {
        return Simd::Mask<V2>(Simd::mask(s1)) || Simd::mask(v2);
      }
 
      template<class V1, class S2>
      auto maskOr(ADLTag<5, VcImpl::IsVector<V1>::value &&
                            std::is_same<Mask<S2>, bool>::value>,
                   const V1 &v1, const S2 &s2)
      {
        return Simd::mask(v1) || Simd::Mask<V1>(Simd::mask(s2));
      }
 
 
    } // namespace Overloads
 
  } // namespace Simd
 
  /*
   * Specialize IsNumber for Vc::SimdArray and Vc::Vector to be able to use
   * it as a scalar in DenseMatrix etc.
   */
  template <typename T, std::size_t N, class V, size_t Wt>
  struct IsNumber<Vc::SimdArray<T, N, V, Wt>>
    : public std::integral_constant<bool, IsNumber<T>::value> {
  };
 
  template <typename T, typename Abi>
  struct IsNumber<Vc::Vector<T, Abi>>
    : public std::integral_constant<bool, IsNumber<T>::value> {
  };
 
  template<class V>
  struct AutonomousValueType<Simd::VcImpl::Proxy<V> > :
    AutonomousValueType<typename Simd::VcImpl::Proxy<V>::value_type> {};
 
} // namespace Dune
 
namespace Vc_VERSIONED_NAMESPACE {
 
 
  template<class T, int N>
  fixed_size_simd<T, N> max(const fixed_size_simd<T, N> &a,
                            const fixed_size_simd<T, N> &b)
  {
    using A = SimdArray<T, std::size_t(N)>;
    return max(static_cast<const A &>(a), static_cast<const A &>(b));
  }
 
 
  template<class T, int N>
  fixed_size_simd<T, N> min(const fixed_size_simd<T, N> &a,
                            const fixed_size_simd<T, N> &b)
  {
    using A = SimdArray<T, std::size_t(N)>;
    return min(static_cast<const A &>(a), static_cast<const A &>(b));
  }
 
} // namespace Vc_VERSIONED_NAMESPACE
 
#endif // DUNE_COMMON_SIMD_VC_HH