From 6107467227abfe235001d870929463d312da51df Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 18 May 2018 16:27:58 +0200
Subject: [PATCH] Implement VectorizedArray for Altivec/VSX (IBM Power) in
 double precision.

Undefine vector, pixel, bool from altivec.h
---
 include/deal.II/base/vectorization.h | 209 +++++++++++++++++++++++++++
 1 file changed, 209 insertions(+)

diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h
index 084980e077..d6ea23541b 100644
--- a/include/deal.II/base/vectorization.h
+++ b/include/deal.II/base/vectorization.h
@@ -56,6 +56,16 @@
 #endif
 
 
+#ifdef __ALTIVEC__
+#include <altivec.h>
+
+// altivec.h defines vector, pixel, bool, but we do not use them, so undefine
+// them before they make trouble
+#undef vector
+#undef pixel
+#undef bool
+#endif
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -2882,6 +2892,205 @@ vectorized_transpose_and_store(const bool                    add_into,
 #endif // if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
 
 
+#if defined(__ALTIVEC__) && defined(__VSX__)
+
+template <>
+class VectorizedArray<double>
+{
+public:
+  /**
+   * This gives the number of vectors collected in this class.
+   */
+  static const unsigned int n_array_elements = 2;
+
+  /**
+   * This function assigns a scalar to this class.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray &
+  operator = (const double x)
+  {
+    data = vec_splats(x);
+    return *this;
+  }
+
+  /**
+   * Access operator. The component must be either 0 or 1.
+   */
+  DEAL_II_ALWAYS_INLINE double &
+  operator [] (const unsigned int comp)
+  {
+    AssertIndexRange (comp, 2);
+    return *(reinterpret_cast<double *>(&data)+comp);
+  }
+
+  /**
+   * Constant access operator.
+   */
+  DEAL_II_ALWAYS_INLINE const double &
+  operator [] (const unsigned int comp) const
+  {
+    AssertIndexRange (comp, 2);
+    return *(reinterpret_cast<const double *>(&data)+comp);
+  }
+
+  /**
+   * Addition.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray &
+  operator += (const VectorizedArray &vec)
+  {
+    data = vec_add(data,vec.data);
+    return *this;
+  }
+
+  /**
+   * Subtraction.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray &
+  operator -= (const VectorizedArray &vec)
+  {
+    data = vec_sub(data,vec.data);
+    return *this;
+  }
+
+  /**
+   * Multiplication.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray &
+  operator *= (const VectorizedArray &vec)
+  {
+    data = vec_mul(data,vec.data);
+    return *this;
+  }
+
+  /**
+   * Division.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray &
+  operator /= (const VectorizedArray &vec)
+  {
+    data = vec_div(data,vec.data);
+    return *this;
+  }
+
+  /**
+   * Load @p n_array_elements from memory into the calling class, starting at
+   * the given address.
+   */
+  DEAL_II_ALWAYS_INLINE void load (const double *ptr)
+  {
+    data = vec_vsx_ld (0, ptr);
+  }
+
+  /**
+   * Write the content of the calling class into memory in form of @p
+   * n_array_elements to the given address.
+   */
+  DEAL_II_ALWAYS_INLINE void store (double *ptr) const
+  {
+    vec_vsx_st (data, 0, ptr);
+  }
+
+  /** @copydoc VectorizedArray<Number>::streaming_store()
+   */
+  DEAL_II_ALWAYS_INLINE
+  void streaming_store (double *ptr) const
+  {
+    store(ptr);
+  }
+
+  /** @copydoc VectorizedArray<Number>::gather()
+   */
+  DEAL_II_ALWAYS_INLINE
+  void gather (const double       *base_ptr,
+               const unsigned int *offsets)
+  {
+    for (unsigned int i=0; i<2; ++i)
+      *(reinterpret_cast<double *>(&data)+i) = base_ptr[offsets[i]];
+  }
+
+  /** @copydoc VectorizedArray<Number>::scatter
+   */
+  DEAL_II_ALWAYS_INLINE
+  void scatter (const unsigned int *offsets,
+                double             *base_ptr) const
+  {
+    for (unsigned int i=0; i<2; ++i)
+      base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data)+i);
+  }
+
+  /**
+   * Actual data field. Since this class represents a POD data type, it
+   * remains public.
+   */
+  __vector double data;
+
+private:
+  /**
+   * Return the square root of this field. Not for use in user code. Use
+   * sqrt(x) instead.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray
+  get_sqrt () const
+  {
+    VectorizedArray res;
+    res.data = vec_sqrt(data);
+    return res;
+  }
+
+  /**
+   * Return the absolute value of this field. Not for use in user code. Use
+   * abs(x) instead.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray
+  get_abs () const
+  {
+    VectorizedArray res;
+    res.data = vec_abs (data);
+    return res;
+  }
+
+  /**
+   * Return the component-wise maximum of this field and another one. Not for
+   * use in user code. Use max(x,y) instead.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray
+  get_max (const VectorizedArray &other) const
+  {
+    VectorizedArray res;
+    res.data = vec_max (data, other.data);
+    return res;
+  }
+
+  /**
+   * Return the component-wise minimum of this field and another one. Not for
+   * use in user code. Use min(x,y) instead.
+   */
+  DEAL_II_ALWAYS_INLINE VectorizedArray
+  get_min (const VectorizedArray &other) const
+  {
+    VectorizedArray res;
+    res.data = vec_min (data, other.data);
+    return res;
+  }
+
+  /**
+   * Make a few functions friends.
+   */
+  template <typename Number2> friend VectorizedArray<Number2>
+  std::sqrt (const VectorizedArray<Number2> &);
+  template <typename Number2> friend VectorizedArray<Number2>
+  std::abs  (const VectorizedArray<Number2> &);
+  template <typename Number2> friend VectorizedArray<Number2>
+  std::max  (const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
+  template <typename Number2> friend VectorizedArray<Number2>
+  std::min  (const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
+};
+
+#endif // ALTIVEC && VSX
+
+
+
 /**
  * Relational operator == for VectorizedArray
  *
-- 
2.39.5