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.
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.
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
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
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
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
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
::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]), ...,
::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]), ...,
::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;
}
#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;
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+}
projects / dealii.git / commitdiff