From: kronbichler Date: Mon, 5 May 2014 07:48:39 +0000 (+0000) Subject: Add methods VectorizedArray::load and VectorizedArray store for vectorized read and... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cfd117c691e43002b95ac6cbca81462cc5916698;p=dealii-svn.git Add methods VectorizedArray::load and VectorizedArray store for vectorized read and write of unaligned memory addresses. This also fixes the step-48 test case problems git-svn-id: https://svn.dealii.org/trunk@32886 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/base/vectorization.h b/deal.II/include/deal.II/base/vectorization.h index a8e14d9dac..ca843ed5dd 100644 --- a/deal.II/include/deal.II/base/vectorization.h +++ b/deal.II/include/deal.II/base/vectorization.h @@ -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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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*. + */ + 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 sin (const ::dealii::VectorizedArray &x) { - ::dealii::VectorizedArray 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::n_array_elements]; for (unsigned int i=0; i::n_array_elements; ++i) - sin_val[i] = std::sin(x[i]); - return sin_val; + values[i] = std::sin(x[i]); + ::dealii::VectorizedArray 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 {tan(x[0]), tan(x[1]), ..., - * tan(x[n_array_elements-1])}. + * Computes the cosine of a vectorized data field. The result is returned as + * vectorized array in the form {cos(x[0]), cos(x[1]), ..., + * cos(x[n_array_elements-1])}. * * @relates VectorizedArray */ template inline ::dealii::VectorizedArray - tan (const ::dealii::VectorizedArray &x) + cos (const ::dealii::VectorizedArray &x) { - ::dealii::VectorizedArray tan_val; + Number values[::dealii::VectorizedArray::n_array_elements]; for (unsigned int i=0; i::n_array_elements; ++i) - tan_val[i] = std::tan(x[i]); - return tan_val; + values[i] = std::cos(x[i]); + ::dealii::VectorizedArray 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 {cos(x[0]), cos(x[1]), ..., - * cos(x[n_array_elements-1])}. + * Computes the tangent of a vectorized data field. The result is returned as + * vectorized array in the form {tan(x[0]), tan(x[1]), ..., + * tan(x[n_array_elements-1])}. * * @relates VectorizedArray */ template inline ::dealii::VectorizedArray - cos (const ::dealii::VectorizedArray &x) + tan (const ::dealii::VectorizedArray &x) { - ::dealii::VectorizedArray cos_val; + Number values[::dealii::VectorizedArray::n_array_elements]; for (unsigned int i=0; i::n_array_elements; ++i) - cos_val[i] = std::cos(x[i]); - return cos_val; + values[i] = std::tan(x[i]); + ::dealii::VectorizedArray 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 {exp(x[0]), exp(x[1]), ..., @@ -1819,13 +1980,16 @@ namespace std ::dealii::VectorizedArray exp (const ::dealii::VectorizedArray &x) { - ::dealii::VectorizedArray exp_val; + Number values[::dealii::VectorizedArray::n_array_elements]; for (unsigned int i=0; i::n_array_elements; ++i) - exp_val[i] = std::exp(x[i]); - return exp_val; + values[i] = std::exp(x[i]); + ::dealii::VectorizedArray 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 {log(x[0]), log(x[1]), ..., @@ -1838,10 +2002,12 @@ namespace std ::dealii::VectorizedArray log (const ::dealii::VectorizedArray &x) { - ::dealii::VectorizedArray log_val; + Number values[::dealii::VectorizedArray::n_array_elements]; for (unsigned int i=0; i::n_array_elements; ++i) - log_val[i] = std::log(x[i]); - return log_val; + values[i] = std::log(x[i]); + ::dealii::VectorizedArray out; + out.load(&values[0]); + return out; } diff --git a/tests/base/vectorization_01.cc b/tests/base/vectorization_01.cc index 3d0fa338a1..b0ba4f671f 100644 --- a/tests/base/vectorization_01.cc +++ b/tests/base/vectorization_01.cc @@ -34,7 +34,7 @@ void test () a = Number(2.); b = Number(-1.); for (unsigned int i=0; i 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 (); + test (); deallog.pop(); } diff --git a/tests/base/vectorization_03.cc b/tests/base/vectorization_03.cc index 5b9ac5a144..82936aae7c 100644 --- a/tests/base/vectorization_03.cc +++ b/tests/base/vectorization_03.cc @@ -22,19 +22,70 @@ #include "../tests.h" #include +#include + + +struct Evaluation +{ + VectorizedArray get_value(const unsigned int index) const + { + return values[index]; + } + + void submit_value(const VectorizedArray 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 cartesian_weight; + VectorizedArray jac_weight[1]; + VectorizedArray general_weight[1]; + VectorizedArray values[1]; +}; + + +void initialize(Evaluation &eval) +{ + eval.is_cartesian = true; + eval.cartesian_weight = static_cast(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(rand())/RAND_MAX; +} void test() { - VectorizedArray array[3]; - for (unsigned int i=0; i<3; ++i) - for (unsigned int v=0; v::n_array_elements; ++v) - array[i][v] = static_cast(rand())/RAND_MAX; + Evaluation current, old; + initialize(current); + initialize(old); + VectorizedArray weight; + weight = static_cast(rand())/RAND_MAX; + + VectorizedArray vec; + for (unsigned int v=0; v::n_array_elements; ++v) + vec[v] = static_cast(rand())/RAND_MAX; + + current.values[0] = vec; + old.values[0] = vec * 1.112 - std::max(2.*vec - 1., VectorizedArray()); + + Vector vector(200); + vector = 1.2; + + VectorizedArray cur = current.get_value(0); + VectorizedArray ol = old.get_value(0); + current.submit_value(2. * cur - ol - weight * std::sin(cur), 0); - const VectorizedArray 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::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 index 0000000000..bf5b059307 --- /dev/null +++ b/tests/base/vectorization_04.cc @@ -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 +#include +#include + +#include +#include + +template +void test () +{ + const unsigned int n_vectors = VectorizedArray::n_array_elements; + std::vector values(n_vectors*5); + for (unsigned int i=0; i > 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 stored(n_vectors*5); + for (unsigned int shift=0; shift < n_vectors; ++shift) + { + for (unsigned int i=0; i<4; ++i) + { + VectorizedArray 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 (); + deallog.pop(); + deallog.push("float"); + test (); + 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 (); + deallog.pop(); + + deallog.push("unsigned int"); + test (); + deallog.pop(); +} diff --git a/tests/base/vectorization_04.output b/tests/base/vectorization_04.output new file mode 100644 index 0000000000..057015edea --- /dev/null +++ b/tests/base/vectorization_04.output @@ -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 diff --git a/tests/matrix_free/step-48.cc b/tests/matrix_free/step-48.cc index 243050d7a4..b920f8b52c 100644 --- a/tests/matrix_free/step-48.cc +++ b/tests/matrix_free/step-48.cc @@ -141,9 +141,6 @@ namespace Step48 const VectorizedArray current_value = current.get_value(q); const VectorizedArray old_value = old.get_value(q); - // the first 'dummy' is needed to work around a compiler bug in - // the tester - const VectorizedArray 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 *