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/simd/valueproxy.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 Proxy, class ValueType>
      struct ProxyOperations
      {
        using value_type = ValueType;
 
        // 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>                                    \
        friend auto operator OP(Proxy&& l, const Vc::Vector<T, Abi> &r) \
          -> decltype(std::declval<value_type>() OP r)                  \
        {                                                               \
          return value_type(l) 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>      \
        friend auto                                                     \
        operator OP(Proxy&& l, const Vc::SimdArray<T, n, Vec, m> &r)    \
          -> decltype(std::declval<value_type>() OP r)                  \
        {                                                               \
          return value_type(l) 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(static_cast<const Proxy &>(*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(static_cast<const Proxy &>(*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
      };
 
 
      template<class V>
      using Proxy = Impl::ValueProxy<V, ProxyOperations>;
 
    } // 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> {
  };
 
} // 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