]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add methods VectorizedArray::load and VectorizedArray store for vectorized read and...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 5 May 2014 07:48:39 +0000 (07:48 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 5 May 2014 07:48:39 +0000 (07:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@32886 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/base/vectorization.h
tests/base/vectorization_01.cc
tests/base/vectorization_03.cc
tests/base/vectorization_04.cc [new file with mode: 0644]
tests/base/vectorization_04.output [new file with mode: 0644]
tests/matrix_free/step-48.cc

index a8e14d9dac4534aaa17e24bd80a4ddd5fffe77df..ca843ed5dd0825e880cc6e5b80a9f9a6b3c6c508 100644 (file)
@@ -172,6 +172,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 64 bytes, as opposed
+   * to casting a double address to VectorizedArray<double>*.
+   */
+  void load (const double* ptr)
+  {
+    data = _mm512_loadu_pd (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 64 bytes, as opposed to casting a double address to
+   * VectorizedArray<double>*.
+   */
+  void store (double* ptr) const
+  {
+    _mm512_storeu_pd (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class represents a POD data type, it
    * remains public.
@@ -343,6 +364,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 64 bytes, as opposed
+   * to casting a float address to VectorizedArray<float>*.
+   */
+  void load (const float* ptr)
+  {
+    data = _mm512_loadu_ps (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 64 bytes, as opposed to casting a float address to
+   * VectorizedArray<float>*.
+   */
+  void store (float* ptr) const
+  {
+    _mm512_storeu_ps (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class represents a POD data type, it
    * remains public.
@@ -526,6 +568,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 32 bytes, as opposed
+   * to casting a double address to VectorizedArray<double>*.
+   */
+  void load (const double* ptr)
+  {
+    data = _mm256_loadu_pd (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 32 bytes, as opposed to casting a double address to
+   * VectorizedArray<double>*.
+   */
+  void store (double* ptr) const
+  {
+    _mm256_storeu_pd (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class
    * represents a POD data type, it remains
@@ -704,6 +767,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 32 bytes, as opposed
+   * to casting a float address to VectorizedArray<float>*.
+   */
+  void load (const float* ptr)
+  {
+    data = _mm256_loadu_ps (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 32 bytes, as opposed to casting a float address to
+   * VectorizedArray<float>*.
+   */
+  void store (float* ptr) const
+  {
+    _mm256_storeu_ps (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class
    * represents a POD data type, it remains
@@ -890,6 +974,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 16 bytes, as opposed
+   * to casting a double address to VectorizedArray<double>*.
+   */
+  void load (const double* ptr)
+  {
+    data = _mm_loadu_pd (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 16 bytes, as opposed to casting a double address to
+   * VectorizedArray<double>*.
+   */
+  void store (double* ptr) const
+  {
+    _mm_storeu_pd (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class
    * represents a POD data type, it remains
@@ -1069,6 +1174,27 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by 16 bytes, as opposed
+   * to casting a float address to VectorizedArray<float>*.
+   */
+  void load (const float* ptr)
+  {
+    data = _mm_loadu_ps (ptr);
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * 16 bytes, as opposed to casting a float address to
+   * VectorizedArray<float>*.
+   */
+  void store (float* ptr) const
+  {
+    _mm_storeu_ps (ptr, data);
+  }
+
   /**
    * Actual data field. Since this class
    * represents a POD data type, it remains
@@ -1288,6 +1414,28 @@ public:
     return *this;
   }
 
+  /**
+   * Loads @p n_array_elements from memory into the calling class, starting at
+   * the given address. The memory need not be aligned by the amount of bytes
+   * in the vectorized array, as opposed to casting a double address to
+   * VectorizedArray<double>*.
+   */
+  void load (const Number* ptr)
+  {
+    data = *ptr;
+  }
+
+  /**
+   * Writes the content of the calling class into memory in form of @p
+   * n_array_elements to the given address. The memory need not be aligned by
+   * the amount of bytes in the vectorized array, as opposed to casting a
+   * double address to VectorizedArray<double>*.
+   */
+  void store (Number* ptr) const
+  {
+    *ptr = data;
+  }
+
   /**
    * Actual data field. Since this class
    * represents a POD data type, it is declared
@@ -1761,52 +1909,65 @@ namespace std
   ::dealii::VectorizedArray<Number>
   sin (const ::dealii::VectorizedArray<Number> &x)
   {
-    ::dealii::VectorizedArray<Number> sin_val;
+    // put values in an array and later read in that array with an unaligned
+    // read. This should safe some instructions as compared to directly
+    // setting the individual elements and also circumvents a compiler
+    // optimization bug in gcc-4.6 (see also deal.II developers list from
+    // April 2014, topic "matrix_free/step-48 Test").
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
     for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
-      sin_val[i] = std::sin(x[i]);
-    return sin_val;
+      values[i] = std::sin(x[i]);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
   }
 
 
 
   /**
-   * Computes the tangent of a vectorized data field. The result is returned as
-   * vectorized array in the form <tt>{tan(x[0]), tan(x[1]), ...,
-   * tan(x[n_array_elements-1])}</tt>.
+   * Computes the cosine of a vectorized data field. The result is returned as
+   * vectorized array in the form <tt>{cos(x[0]), cos(x[1]), ...,
+   * cos(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
    */
   template <typename Number>
   inline
   ::dealii::VectorizedArray<Number>
-  tan (const ::dealii::VectorizedArray<Number> &x)
+  cos (const ::dealii::VectorizedArray<Number> &x)
   {
-    ::dealii::VectorizedArray<Number> tan_val;
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
     for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
-      tan_val[i] = std::tan(x[i]);
-    return tan_val;
+      values[i] = std::cos(x[i]);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
   }
 
 
+
   /**
-   * Computes the cosine of a vectorized data field. The result is returned as
-   * vectorized array in the form <tt>{cos(x[0]), cos(x[1]), ...,
-   * cos(x[n_array_elements-1])}</tt>.
+   * Computes the tangent of a vectorized data field. The result is returned as
+   * vectorized array in the form <tt>{tan(x[0]), tan(x[1]), ...,
+   * tan(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
    */
   template <typename Number>
   inline
   ::dealii::VectorizedArray<Number>
-  cos (const ::dealii::VectorizedArray<Number> &x)
+  tan (const ::dealii::VectorizedArray<Number> &x)
   {
-    ::dealii::VectorizedArray<Number> cos_val;
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
     for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
-      cos_val[i] = std::cos(x[i]);
-    return cos_val;
+      values[i] = std::tan(x[i]);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
   }
 
 
+
   /**
    * Computes the exponential of a vectorized data field. The result is returned
    * as vectorized array in the form <tt>{exp(x[0]), exp(x[1]), ...,
@@ -1819,13 +1980,16 @@ namespace std
   ::dealii::VectorizedArray<Number>
   exp (const ::dealii::VectorizedArray<Number> &x)
   {
-    ::dealii::VectorizedArray<Number> exp_val;
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
     for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
-      exp_val[i] = std::exp(x[i]);
-    return exp_val;
+      values[i] = std::exp(x[i]);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
   }
 
 
+
   /**
    * Computes the natural logarithm of a vectorized data field. The result is
    * returned as vectorized array in the form <tt>{log(x[0]), log(x[1]), ...,
@@ -1838,10 +2002,12 @@ namespace std
   ::dealii::VectorizedArray<Number>
   log (const ::dealii::VectorizedArray<Number> &x)
   {
-    ::dealii::VectorizedArray<Number> log_val;
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
     for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
-      log_val[i] = std::log(x[i]);
-    return log_val;
+      values[i] = std::log(x[i]);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
   }
 
 
index 3d0fa338a1c0b7376d7446cb89cb900416950620..b0ba4f671fe6926ccff5135edfbbc45679aadab4 100644 (file)
@@ -34,7 +34,7 @@ void test ()
   a = Number(2.);
   b = Number(-1.);
   for (unsigned int i=0; i<n_vectors; ++i)
-    c[i] = i;
+    c[i] = Number(i);
 
   deallog << "Addition: ";
   VectorizedArray<Number> d = a + b;
@@ -57,17 +57,17 @@ void test ()
     AssertThrow (e[i] == c[i], ExcInternalError());
   deallog << "OK" << std::endl
           << "Multiplication scalar: ";
-  a = 2. * a;
+  a = Number(2.) * a;
   for (unsigned int i=0; i<n_vectors; ++i)
     AssertThrow (a[i] == 4., ExcInternalError());
   deallog << "OK" << std::endl
           << "Division scalar left: ";
-  d = 1. / a;
+  d = Number(1.) / a;
   for (unsigned int i=0; i<n_vectors; ++i)
     AssertThrow (d[i] == 0.25, ExcInternalError());
   deallog << "OK" << std::endl
           << "Division scalar right: ";
-  e = d / 0.25;
+  e = d / Number(0.25);
   for (unsigned int i=0; i<n_vectors; ++i)
     AssertThrow (e[i] == 1, ExcInternalError());
   deallog << "OK" << std::endl
@@ -103,7 +103,7 @@ void test ()
     AssertThrow (d[i] == std::max(a[i], c[i]), ExcInternalError());
   deallog << "OK" << std::endl
           << "Minimum value: ";
-  d = std::min(0.5 * a + b, c);
+  d = std::min(Number(0.5) * a + b, c);
   for (unsigned int i=0; i<n_vectors; ++i)
     AssertThrow (d[i] == std::min(Number(0.5 * a[i] + b[i]), c[i]),
                  ExcInternalError());
@@ -167,6 +167,6 @@ int main()
   // path of VectorizedArray is taken no matter
   // what was done for double or float
   deallog.push("long double");
-  test<float> ();
+  test<long double> ();
   deallog.pop();
 }
index 5b9ac5a1446f1554acaab402c8a847466eb3cee1..82936aae7cb8dd47903a66177ad4fe05e67d3209 100644 (file)
 #include "../tests.h"
 
 #include <deal.II/base/vectorization.h>
+#include <deal.II/lac/vector.h>
+
+
+struct Evaluation
+{
+  VectorizedArray<double> get_value(const unsigned int index) const
+  {
+    return values[index];
+  }
+
+  void submit_value(const VectorizedArray<double> val,
+                    const unsigned int index)
+  {
+    if (is_cartesian)
+      values[index] = val * cartesian_weight * jac_weight[index];
+    else
+      values[index] = val * general_weight[index];
+  }
+
+  bool is_cartesian;
+  VectorizedArray<double> cartesian_weight;
+  VectorizedArray<double> jac_weight[1];
+  VectorizedArray<double> general_weight[1];
+  VectorizedArray<double> values[1];
+};
+
+
+void initialize(Evaluation &eval)
+{
+  eval.is_cartesian = true;
+  eval.cartesian_weight = static_cast<double>(rand())/RAND_MAX;
+  for (unsigned int i=0; i<4; ++i)
+    eval.cartesian_weight = std::max(eval.cartesian_weight, eval.cartesian_weight * eval.cartesian_weight);
+  eval.general_weight[0] = 0.2313342 * eval.cartesian_weight;
+  eval.jac_weight[0] = static_cast<double>(rand())/RAND_MAX;
+}
 
 
 void test()
 {
-  VectorizedArray<double> array[3];
-  for (unsigned int i=0; i<3; ++i)
-    for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
-      array[i][v] = static_cast<double>(rand())/RAND_MAX;
+  Evaluation current, old;
+  initialize(current);
+  initialize(old);
+  VectorizedArray<double> weight;
+  weight = static_cast<double>(rand())/RAND_MAX;
+
+  VectorizedArray<double> vec;
+  for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+    vec[v] = static_cast<double>(rand())/RAND_MAX;
+
+  current.values[0] = vec;
+  old.values[0] = vec * 1.112 - std::max(2.*vec - 1., VectorizedArray<double>());
+
+  Vector<double> vector(200);
+  vector = 1.2;
+
+  VectorizedArray<double> cur = current.get_value(0);
+  VectorizedArray<double> ol =  old.get_value(0);
+  current.submit_value(2. * cur - ol - weight * std::sin(cur), 0);
 
-  const VectorizedArray<double> tmp = 2. * array[0] - array[1] - array[2] * std::sin(array[0]);
+  vector *= 2.*current.get_value(0)[0];
 
   for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
-    deallog << tmp[v]-(2.*array[0][v]-array[1][v]-array[2][v]*std::sin(array[0][v])) << " ";
+    deallog << current.get_value(0)[v]/(current.cartesian_weight[v]*current.jac_weight[0][0])-(2.*vec[v]-ol[v]-weight[v]*std::sin(vec[v]) + std::cos(ol[v])) << " ";
   deallog << std::endl;
 }
 
diff --git a/tests/base/vectorization_04.cc b/tests/base/vectorization_04.cc
new file mode 100644 (file)
index 0000000..bf5b059
--- /dev/null
@@ -0,0 +1,91 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test for VectorizedArray::load and VectorizedArray::store
+
+#include "../tests.h"
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/vectorization.h>
+#include <deal.II/base/aligned_vector.h>
+
+template <typename Number>
+void test ()
+{
+  const unsigned int n_vectors = VectorizedArray<Number>::n_array_elements;
+  std::vector<Number> values(n_vectors*5);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = i;
+  AlignedVector<VectorizedArray<Number> > copied(4);
+
+  // test load operation for all possible values of alignment
+  for (unsigned int shift=0; shift < n_vectors; ++shift)
+    {
+      for (unsigned int i=0; i<4; ++i)
+        copied[i].load(&values[i*n_vectors+shift]);
+      for (unsigned int i=0; i<4; ++i)
+        for (unsigned int v=0; v<n_vectors; ++v)
+          AssertThrow(copied[i][v] == values[i*n_vectors+v+shift],
+                      ExcInternalError());
+    }
+  deallog << "load OK" << std::endl;
+
+  // test store operation
+  std::vector<Number> stored(n_vectors*5);
+  for (unsigned int shift=0; shift < n_vectors; ++shift)
+    {
+      for (unsigned int i=0; i<4; ++i)
+        {
+          VectorizedArray<Number> tmp;
+          tmp.load(&values[i*n_vectors]);
+          tmp.store(&stored[i*n_vectors+shift]);
+        }
+      for (unsigned int i=0; i<4*n_vectors; ++i)
+        AssertThrow(stored[i+shift] == i, ExcInternalError());
+    }
+  deallog << "store OK" << std::endl;
+}
+
+
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog.push("double");
+  test<double> ();
+  deallog.pop();
+  deallog.push("float");
+  test<float> ();
+  deallog.pop();
+
+  // test long double and unsigned int: in these cases, the default path of
+  // VectorizedArray is taken no matter what was done for double or float
+  deallog.push("long double");
+  test<long double> ();
+  deallog.pop();
+
+  deallog.push("unsigned int");
+  test<unsigned int> ();
+  deallog.pop();
+}
diff --git a/tests/base/vectorization_04.output b/tests/base/vectorization_04.output
new file mode 100644 (file)
index 0000000..057015e
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL:double::load OK
+DEAL:double::store OK
+DEAL:float::load OK
+DEAL:float::store OK
+DEAL:long double::load OK
+DEAL:long double::store OK
+DEAL:unsigned int::load OK
+DEAL:unsigned int::store OK
index 243050d7a4584e6640493bef3c9cf1f3abdc1c2c..b920f8b52c49d19621ff524a07dfc629805d532f 100644 (file)
@@ -141,9 +141,6 @@ namespace Step48
             const VectorizedArray<double> current_value = current.get_value(q);
             const VectorizedArray<double> old_value     = old.get_value(q);
 
-            // the first 'dummy' is needed to work around a compiler bug in
-            // the tester
-            const VectorizedArray<double> dummy = -std::sin(current_value);
             current.submit_value (2.*current_value - old_value -
                                   delta_t_sqr * std::sin(current_value),q);
             current.submit_gradient (- delta_t_sqr *

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.