From 989fef6e50ff05510809917a8be61b62a0b99ff0 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 18 Feb 2014 21:37:44 +0000 Subject: [PATCH] Put AlignedVector into its own file git-svn-id: https://svn.dealii.org/trunk@32508 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/base/aligned_vector.h | 776 +++++++++++++++ deal.II/include/deal.II/base/vectorization.h | 913 +++++++----------- deal.II/source/distributed/tria.cc | 2 +- 3 files changed, 1132 insertions(+), 559 deletions(-) create mode 100644 deal.II/include/deal.II/base/aligned_vector.h diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h new file mode 100644 index 0000000000..c013786d69 --- /dev/null +++ b/deal.II/include/deal.II/base/aligned_vector.h @@ -0,0 +1,776 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2011 - 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. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__aligned_vector_h +#define __deal2__aligned_vector_h + +#include +#include +#include +#include +#include +#include + + +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 +#include +#endif + + + +DEAL_II_NAMESPACE_OPEN + + +/** + * This is a replacement class for std::vector to be used in combination with + * VectorizedArray and derived data types. It allocates memory aligned to + * addresses of a vectorized data type (in order to avoid segmentation faults + * when a variable of type VectorizedArray which the compiler assumes to be + * aligned to certain memory addresses does not actually follow these + * rules). This could also be achieved by proving std::vector with a + * user-defined allocator. On the other hand, writing an own small vector + * class lets us implement parallel copy and move operations with TBB, insert + * deal.II-style assertions, and cut some unnecessary functionality. Note that + * this vector is a bit more memory-consuming than std::vector because of + * alignment, so it is recommended to only use this vector on long vectors. + * + * @p author Katharina Kormann, Martin Kronbichler, 2011 + */ +template < class T > +class AlignedVector +{ +public: + /** + * Declare standard types used in all containers. These types parallel those + * in the C++ standard libraries vector<...> class. + */ + typedef T value_type; + typedef value_type *pointer; + typedef const value_type *const_pointer; + typedef value_type *iterator; + typedef const value_type *const_iterator; + typedef value_type &reference; + typedef const value_type &const_reference; + typedef std::size_t size_type; + + /** + * Empty constructor. Sets the vector size to zero. + */ + AlignedVector (); + + /** + * Sets the vector size to the given size and initializes all elements with + * T(). + */ + AlignedVector (const size_type size); + + /** + * Destructor. + */ + ~AlignedVector (); + + /** + * Copy constructor. + */ + AlignedVector (const AlignedVector &vec); + + /** + * Assignment to the input vector @p vec. + */ + AlignedVector & + operator = (const AlignedVector &vec); + + /** + * Change the size of the vector. It keeps old elements previously available + * but does not initialize the newly allocated memory, leaving it in an + * undefined state. + */ + void resize_fast (const size_type size); + + /** + * Change the size of the vector. It keeps old elements previously + * available, and initializes each element with the specified data. If the + * new vector size is shorter than the old one, the memory is not released + * unless the new size is zero. + */ + void resize (const size_type size_in, + const T &init = T()); + + /** + * Reserve memory space for @p size elements. If the argument @p size is set + * to zero, all previously allocated memory is released. + * + * In order to avoid too frequent reallocation (which involves copy of the + * data), this function doubles the amount of memory occupied when the given + * size is larger than the previously allocated size. + */ + void reserve (const size_type size_alloc); + + /** + * Releases all previously allocated memory and leaves the vector in a state + * equivalent to the state after the default constructor has been called. + */ + void clear (); + + /** + * Inserts an element at the end of the vector, increasing the vector size + * by one. Note that the allocated size will double whenever the previous + * space is not enough to hold the new element. + */ + void push_back (const T in_data); + + /** + * Returns the last element of the vector (read and write access). + */ + reference back (); + + /** + * Returns the last element of the vector (read-only access). + */ + const_reference back () const; + + /** + * Inserts several elements at the end of the vector given by a range of + * elements. + */ + template + void insert_back (ForwardIterator begin, + ForwardIterator end); + + /** + * Swaps the given vector with the calling vector. + */ + void swap (AlignedVector &vec); + + /** + * Returns whether the vector is empty, i.e., its size is zero. + */ + bool empty () const; + + /** + * Returns the size of the vector. + */ + size_type size () const; + + /** + * Returns the capacity of the vector, i.e., the size this vector can hold + * without reallocation. Note that capacity() >= size(). + */ + size_type capacity () const; + + /** + * Read-write access to entry @p index in the vector. + */ + reference + operator [] (const size_type index); + + /** + * Read-only access to entry @p index in the vector. + */ + const_reference operator [] (const size_type index) const; + + /** + * Returns a read and write pointer to the beginning of the data array. + */ + iterator begin (); + + /** + * Returns a read and write pointer to the end of the data array. + */ + iterator end (); + + /** + * Returns a read-only pointer to the beginning of the data array. + */ + const_iterator begin () const; + + /** + * Returns a read-only pointer to the end of the data array. + */ + const_iterator end () const; + + /** + * Returns the memory consumption of the allocated memory in this class. If + * the underlying type @p T allocates memory by itself, this memory is not + * counted. + */ + size_type memory_consumption () const; + +private: + + /** + * Pointer to actual class data. + */ + T *_data; + + /** + * Pointer to the end of valid data fields. + */ + T *_end_data; + + /** + * Pointer to the end of the allocated memory. + */ + T *_end_allocated; +}; + + +// ------------------------------- inline functions -------------------------- + +/** + * This namespace defines the copy and set functions used in + * AlignedVector. These functions operate in parallel when there are enough + * elements in the vector. + */ +namespace internal +{ + /** + * Move and class that actually issues the copy commands in + * AlignedVector. This class is based on the specialized for loop base class + * ParallelForLoop in parallel.h whose purpose is the following: When + * calling a parallel for loop on AlignedVector with apply_to_subranges, it + * generates different code for every different argument we might choose (as + * it is templated). This gives a lot of code (e.g. it triples the memory + * required for compiling the file matrix_free.cc and the final object size + * is several times larger) which is completely useless. Therefore, this + * class channels all copy commands through one call to apply_to_subrange + * for all possible types, which makes the copy operation much cleaner + * (thanks to a virtual function, whose cost is negligible in this context). + * + * @relates AlignedVector + */ + template + class AlignedVectorMove : private parallel::ParallelForInteger + { + static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1; + public: + /** + * Constructor. Issues a parallel call if + * there are sufficiently many elements, + * otherwise work in serial. Copies the data + * from source to destination and then calls + * destructor on the source. If the optional + * argument is set to true, the source is left + * untouched instead. + */ + AlignedVectorMove (T *source_begin, + T *source_end, + T *destination, + bool copy_only = false) + : + source_ (source_begin), + destination_ (destination), + copy_only_ (copy_only) + { + Assert (source_end >= source_begin, ExcInternalError()); + const std::size_t size = source_end - source_begin; + if (size < minimum_parallel_grain_size) + apply_to_subrange (0, size); + else + apply_parallel (0, size, minimum_parallel_grain_size); + } + + /** + * This method moves elements from the source + * to the destination given in the constructor + * on a subrange given by two integers. + */ + virtual void apply_to_subrange (const std::size_t begin, + const std::size_t end) const + { + // for classes trivial assignment can use + // memcpy + if (std_cxx1x::is_trivial::value == true) + std::memcpy (destination_+begin, source_+begin, (end-begin)*sizeof(T)); + else if (copy_only_ == false) + for (std::size_t i=begin; i + class AlignedVectorSet : private parallel::ParallelForInteger + { + static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1; + public: + /** + * Constructor. Issues a parallel call if + * there are sufficiently many elements, + * otherwise work in serial. + */ + AlignedVectorSet (const std::size_t size, + const T &element, + T *destination) + : + element_ (element), + destination_ (destination), + trivial_element (false) + { + if (size == 0) + return; + + if (std_cxx1x::is_trivial::value == true) + { + const unsigned char zero [sizeof(T)] = {}; + if (std::memcmp(zero, &element, sizeof(T)) == 0) + trivial_element = true; + } + if (size < minimum_parallel_grain_size) + apply_to_subrange (0, size); + else + apply_parallel (0, size, minimum_parallel_grain_size); + } + + private: + + /** + * This sets elements on a subrange given by + * two integers. + */ + virtual void apply_to_subrange (const std::size_t begin, + const std::size_t end) const + { + // for classes with trivial assignment of zero + // can use memset + if (std_cxx1x::is_trivial::value == true && trivial_element) + std::memset (destination_+begin, 0, (end-begin)*sizeof(T)); + else + for (std::size_t i=begin; i +inline +AlignedVector::AlignedVector () + : + _data (0), + _end_data (0), + _end_allocated (0) +{} + + + +template < class T > +inline +AlignedVector::AlignedVector (const size_type size) + : + _data (0), + _end_data (0), + _end_allocated (0) +{ + if (size > 0) + resize (size); +} + + + +template < class T > +inline +AlignedVector::~AlignedVector () +{ + clear(); +} + + + +template < class T > +inline +AlignedVector::AlignedVector (const AlignedVector &vec) + : + _data (0), + _end_data (0), + _end_allocated (0) +{ + // do not invalidate old data + resize_fast (vec._end_data - vec._data); + internal::AlignedVectorMove (vec._data, vec._end_data, _data, true); +} + + + +template < class T > +inline +AlignedVector& +AlignedVector::operator = (const AlignedVector &vec) +{ + clear(); + resize_fast (vec._end_data - vec._data); + internal::AlignedVectorMove (vec._data, vec._end_data, _data, true); + return *this; +} + + + +template < class T > +inline +void +AlignedVector::resize_fast (const size_type size) +{ + reserve (size); + _end_data = _data + size; +} + + +template < class T > +inline +void +AlignedVector::resize (const size_type size_in, + const T &init) +{ + const size_type old_size = size(); + if (std_cxx1x::is_trivial::value == false && size_in < old_size) + { + // call destructor on fields that are released + while (_end_data != _data+size_in) + (--_end_data)->~T(); + } + + resize_fast (size_in); + // now _size is set correctly, need to set the + // values + if (size_in > old_size) + internal::AlignedVectorSet (size_in-old_size, init, + _data+old_size); +} + + + +template < class T > +inline +void +AlignedVector::reserve (const size_type size_alloc) +{ + const size_type old_size = _end_data - _data; + const size_type allocated_size = _end_allocated - _data; + if (size_alloc > allocated_size) + { + // if we continuously increase the size of the vector, we might be + // reallocating a lot of times. therefore, try to increase the size more + // aggressively + size_type new_size = size_alloc; + if (size_alloc < (2 * allocated_size)) + new_size = 2 * allocated_size; + + const size_type size_actual_allocate = new_size * sizeof(T); + +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 + + // allocate and align along boundaries of the size of + // VectorizedArray, which is 16 bytes for SSE and 32 bytes for + // AVX + T *new_data = static_cast(_mm_malloc (size_actual_allocate, + sizeof(VectorizedArray))); +#else + T *new_data = static_cast(malloc (size_actual_allocate)); +#endif + if (new_data == 0) + throw std::bad_alloc(); + + // copy data in case there was some content before and release the old + // memory with the function corresponding to the one used for allocating + std::swap (_data, new_data); + _end_data = _data + old_size; + _end_allocated = _data + new_size; + if (_end_data != _data) + { + internal::AlignedVectorMove(new_data, new_data + old_size, + _data); +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 + _mm_free(new_data); +#else + free(new_data); +#endif + } + } + else if (size_alloc == 0) + clear(); +} + + + +template < class T > +inline +void +AlignedVector::clear () +{ + if (_data != 0) + { + if (std_cxx1x::is_trivial::value == false) + while (_end_data != _data) + (--_end_data)->~T(); + +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 + _mm_free(_data); +#else + free(_data); +#endif + } + _data = 0; + _end_data = 0; + _end_allocated = 0; +} + + + +template < class T > +inline +void +AlignedVector::push_back (const T in_data) +{ + Assert (_end_data <= _end_allocated, ExcInternalError()); + if (_end_data == _end_allocated) + reserve (std::max(2*capacity(),static_cast(16))); + if (std_cxx1x::is_trivial::value == false) + new (_end_data) T; + *_end_data++ = in_data; +} + + + +template < class T > +inline +typename AlignedVector::reference +AlignedVector::back () +{ + AssertIndexRange (0, size()); + T *field = _end_data - 1; + return *field; +} + + + +template < class T > +inline +typename AlignedVector::const_reference +AlignedVector::back () const +{ + AssertIndexRange (0, size()); + const T *field = _end_data - 1; + return *field; +} + + + +template < class T > +template +inline +void +AlignedVector::insert_back (ForwardIterator begin, + ForwardIterator end) +{ + const unsigned int old_size = size(); + reserve (old_size + (end-begin)); + for ( ; begin != end; ++begin, ++_end_data) + { + if (std_cxx1x::is_trivial::value == false) + new (_end_data) T; + *_end_data = *begin; + } +} + + + +template < class T > +inline +void +AlignedVector::swap (AlignedVector &vec) +{ + std::swap (_data, vec._data); + std::swap (_end_data, vec._end_data); + std::swap (_end_allocated, vec._end_allocated); +} + + + +template < class T > +inline +bool +AlignedVector::empty () const +{ + return _end_data == _data; +} + + + +template < class T > +inline +typename AlignedVector::size_type +AlignedVector::size () const +{ + return _end_data - _data; +} + + + +template < class T > +inline +typename AlignedVector::size_type +AlignedVector::capacity () const +{ + return _end_allocated - _data; +} + + + +template < class T > +inline +typename AlignedVector::reference +AlignedVector::operator [] (const size_type index) +{ + AssertIndexRange (index, size()); + return _data[index]; +} + + + +template < class T > +inline +typename AlignedVector::const_reference +AlignedVector::operator [] (const size_type index) const +{ + AssertIndexRange (index, size()); + return _data[index]; +} + + + +template < class T > +inline +typename AlignedVector::iterator +AlignedVector::begin () +{ + return _data; +} + + + +template < class T > +inline +typename AlignedVector::iterator +AlignedVector::end () +{ + return _end_data; +} + + + +template < class T > +inline +typename AlignedVector::const_iterator +AlignedVector::begin () const +{ + return _data; +} + + + +template < class T > +inline +typename AlignedVector::const_iterator +AlignedVector::end () const +{ + return _end_data; +} + + + +template < class T > +inline +typename AlignedVector::size_type +AlignedVector::memory_consumption () const +{ + size_type memory = sizeof(this); + memory += sizeof(T) * capacity(); + return memory; +} + + +#endif // ifndef DOXYGEN + + +/** + * Relational operator == for AlignedVector + * + * @relates AlignedVector + */ +template < class T > +bool operator == (const AlignedVector &lhs, + const AlignedVector &rhs) +{ + if (lhs.size() != rhs.size()) + return false; + for (typename AlignedVector::const_iterator lit = lhs.begin(), + rit = rhs.begin(); lit != lhs.end(); ++lit, ++rit) + if (*lit != *rit) + return false; + return true; +} + + + + +/** + * Relational operator != for AlignedVector + * + * @relates AlignedVector + */ +template < class T > +bool operator != (const AlignedVector &lhs, + const AlignedVector &rhs) +{ + return !(operator==(lhs, rhs)); +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/include/deal.II/base/vectorization.h b/deal.II/include/deal.II/base/vectorization.h index ec48938935..cdb8017464 100644 --- a/deal.II/include/deal.II/base/vectorization.h +++ b/deal.II/include/deal.II/base/vectorization.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -19,17 +19,15 @@ #define __deal2__vectorization_h #include -#include #include -#include -#include #include -#include // Note: // The flag DEAL_II_COMPILER_VECTORIZATION_LEVEL is essentially constructed // according to the following scheme +// #ifdef __AVX512F__ +// #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 3 // #ifdef __AVX__ // #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 2 // #elif defined (__SSE2__) @@ -40,12 +38,10 @@ // In addition to checking the flags __AVX__ and __SSE2__, a configure test // ensures that these feature are not only present but also working properly. -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL == 2 // AVX +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 // AVX, AVX-512 #include -#include #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL == 1 // SSE2 #include -#include #endif @@ -71,10 +67,359 @@ namespace std DEAL_II_NAMESPACE_OPEN -// for safety, also check that __AVX__ is defined in case the user manually +// for safety, also check that __AVX512F__ is defined in case the user manually // set some conflicting compile flags which prevent compilation -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL == 2 && defined(__AVX__) +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL == 3 && defined(__AVX512F__) + +/** + * Specialization of VectorizedArray class for double and AVX-512. + */ +template <> +class VectorizedArray +{ +public: + /** + * This gives the number of vectors collected in this class. + */ + static const unsigned int n_array_elements = 8; + + /** + * This function can be used to set all data fields to a given scalar. + */ + VectorizedArray & + operator = (const double x) + { + data = _mm256_set_pd(x, x, x, x, x, x, x, x); + return *this; + } + + /** + * Access operator. + */ + double & + operator [] (const unsigned int comp) + { + AssertIndexRange (comp, 8); + return *(reinterpret_cast(&data)+comp); + } + + /** + * Constant access operator. + */ + const double & + operator [] (const unsigned int comp) const + { + AssertIndexRange (comp, 8); + return *(reinterpret_cast(&data)+comp); + } + + /** + * Addition. + */ + VectorizedArray & + operator += (const VectorizedArray &vec) + { + // if the compiler supports vector arithmetics, we can simply use += + // operator on the given data type. Otherwise, we need to use the built-in + // intrinsic command for __m512d +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data += vec.data; +#else + data = _mm512_add_pd(data,vec.data); +#endif + return *this; + } + + /** + * Subtraction. + */ + VectorizedArray & + operator -= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data -= vec.data; +#else + data = _mm512_sub_pd(data,vec.data); +#endif + return *this; + } + /** + * Multiplication. + */ + VectorizedArray & + operator *= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data *= vec.data; +#else + data = _mm512_mul_pd(data,vec.data); +#endif + return *this; + } + + /** + * Division. + */ + VectorizedArray & + operator /= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data /= vec.data; +#else + data = _mm512_div_pd(data,vec.data); +#endif + return *this; + } + + /** + * Actual data field. Since this class represents a POD data type, it + * remains public. + */ + __m512d data; + +private: + /** + * Returns the square root of this field. Not for use in user code. Use + * sqrt(x) instead. + */ + VectorizedArray + get_sqrt () const + { + VectorizedArray res; + res.data = _mm512_sqrt_pd(data); + return res; + } + + /** + * Returns the absolute value of this field. Not for use in user code. Use + * abs(x) instead. + */ + VectorizedArray + get_abs () const + { + // to compute the absolute value, perform bitwise andnot with -0. This + // will leave all value and exponent bits unchanged but force the sign + // value to +. + __m256d mask = _mm256_set_pd (-0., -0., -0., -0., -0., -0., -0., -0.); + VectorizedArray res; + res.data = _mm256_andnot_pd(mask, data); + return res; + } + + /** + * Returns the component-wise maximum of this field and another one. Not for + * use in user code. Use max(x,y) instead. + */ + VectorizedArray + get_max (const VectorizedArray &other) const + { + VectorizedArray res; + res.data = _mm512_max_pd (data, other.data); + return res; + } + + /** + * Returns the component-wise minimum of this field and another one. Not for + * use in user code. Use min(x,y) instead. + */ + VectorizedArray + get_min (const VectorizedArray &other) const + { + VectorizedArray res; + res.data = _mm512_min_pd (data, other.data); + return res; + } + + /** + * Make a few functions friends. + */ + template friend VectorizedArray + std::sqrt (const VectorizedArray &); + template friend VectorizedArray + std::abs (const VectorizedArray &); + template friend VectorizedArray + std::max (const VectorizedArray &, const VectorizedArray &); + template friend VectorizedArray + std::min (const VectorizedArray &, const VectorizedArray &); +}; + + + +/** + * Specialization for float and AVX. + */ +template<> +class VectorizedArray +{ +public: + /** + * This gives the number of vectors collected in this class. + */ + static const unsigned int n_array_elements = 16; + + /** + * This function can be used to set all data fields to a given scalar. + */ + VectorizedArray & + operator = (const float x) + { + data = _mm256_set_ps(x, x, x, x, x, x, x, x, x, x, x, x, x, x, x, x); + return *this; + } + + /** + * Access operator. + */ + float & + operator [] (const unsigned int comp) + { + AssertIndexRange (comp, 16); + return *(reinterpret_cast(&data)+comp); + } + + /** + * Constant access operator. + */ + const float & + operator [] (const unsigned int comp) const + { + AssertIndexRange (comp, 16); + return *(reinterpret_cast(&data)+comp); + } + + /** + * Addition. + */ + VectorizedArray & + operator += (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data += vec.data; +#else + data = _mm512_add_ps(data,vec.data); +#endif + return *this; + } + + /** + * Subtraction. + */ + VectorizedArray & + operator -= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data -= vec.data; +#else + data = _mm512_sub_ps(data,vec.data); +#endif + return *this; + } + /** + * Multiplication. + */ + VectorizedArray & + operator *= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data *= vec.data; +#else + data = _mm512_mul_ps(data,vec.data); +#endif + return *this; + } + + /** + * Division. + */ + VectorizedArray & + operator /= (const VectorizedArray &vec) + { +#ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS + data /= vec.data; +#else + data = _mm512_div_ps(data,vec.data); +#endif + return *this; + } + + /** + * Actual data field. Since this class represents a POD data type, it + * remains public. + */ + __m512 data; + +private: + + /** + * Returns the square root of this field. Not for use in user code. Use + * sqrt(x) instead. + */ + VectorizedArray + get_sqrt () const + { + VectorizedArray res; + res.data = _mm512_sqrt_ps(data); + return res; + } + + /** + * Returns the absolute value of this field. Not for use in user code. Use + * abs(x) instead. + */ + VectorizedArray + get_abs () const + { + // to compute the absolute value, perform + // bitwise andnot with -0. This will leave all + // value and exponent bits unchanged but force + // the sign value to +. + __m256 mask = _mm512_setzero_ps (-0.f, -0.f, -0.f, -0.f, -0.f, -0.f, -0.f, -0.f, + -0.f, -0.f, -0.f, -0.f, -0.f, -0.f, -0.f, -0.f); + VectorizedArray res; + res.data = _mm512_andnot_ps(mask, data); + return res; + } + + /** + * Returns the component-wise maximum of this field and another one. Not for + * use in user code. Use max(x,y) instead. + */ + VectorizedArray + get_max (const VectorizedArray &other) const + { + VectorizedArray res; + res.data = _mm512_max_ps (data, other.data); + return res; + } + + /** + * Returns the component-wise minimum of this field and another one. Not for + * use in user code. Use min(x,y) instead. + */ + VectorizedArray + get_min (const VectorizedArray &other) const + { + VectorizedArray res; + res.data = _mm512_min_ps (data, other.data); + return res; + } + + /** + * Make a few functions friends. + */ + template friend VectorizedArray + std::sqrt (const VectorizedArray &); + template friend VectorizedArray + std::abs (const VectorizedArray &); + template friend VectorizedArray + std::max (const VectorizedArray &, const VectorizedArray &); + template friend VectorizedArray + std::min (const VectorizedArray &, const VectorizedArray &); +}; + + +#elif DEAL_II_COMPILER_VECTORIZATION_LEVEL == 2 && defined(__AVX__) /** * Specialization of VectorizedArray class for double and AVX. @@ -1392,554 +1737,6 @@ operator - (const VectorizedArray &u) -/** - * This namespace defines the copy and set functions used in - * AlignedVector. These functions operate in parallel when there are enough - * elements in the vector. - */ -namespace internal -{ - /** - * Move and class that actually issues the copy commands in - * AlignedVector. This class is based on the specialized for loop base class - * ParallelForLoop in parallel.h whose purpose is the following: When - * calling a parallel for loop on AlignedVector with apply_to_subranges, it - * generates different code for every different argument we might choose (as - * it is templated). This gives a lot of code (e.g. it triples the memory - * required for compiling the file matrix_free.cc and the final object size - * is several times larger) which is completely useless. Therefore, this - * class channels all copy commands through one call to apply_to_subrange - * for all possible types, which makes the copy operation much cleaner - * (thanks to a virtual function, whose cost is negligible in this context). - * - * @relates AlignedVector - */ - template - class AlignedVectorMove : private parallel::ParallelForInteger - { - static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1; - public: - /** - * Constructor. Issues a parallel call if - * there are sufficiently many elements, - * otherwise work in serial. Copies the data - * from source to destination and then calls - * destructor on the source. If the optional - * argument is set to true, the source is left - * untouched instead. - */ - AlignedVectorMove (T *source_begin, - T *source_end, - T *destination, - bool copy_only = false) - : - source_ (source_begin), - destination_ (destination), - copy_only_ (copy_only) - { - Assert (source_end >= source_begin, ExcInternalError()); - const std::size_t size = source_end - source_begin; - if (size < minimum_parallel_grain_size) - apply_to_subrange (0, size); - else - apply_parallel (0, size, minimum_parallel_grain_size); - } - - /** - * This method moves elements from the source - * to the destination given in the constructor - * on a subrange given by two integers. - */ - virtual void apply_to_subrange (const std::size_t begin, - const std::size_t end) const - { - // for classes trivial assignment can use - // memcpy - if (std_cxx1x::is_trivial::value == true) - std::memcpy (destination_+begin, source_+begin, (end-begin)*sizeof(T)); - else if (copy_only_ == false) - for (std::size_t i=begin; i - class AlignedVectorSet : private parallel::ParallelForInteger - { - static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1; - public: - /** - * Constructor. Issues a parallel call if - * there are sufficiently many elements, - * otherwise work in serial. - */ - AlignedVectorSet (const std::size_t size, - const T &element, - T *destination) - : - element_ (element), - destination_ (destination), - trivial_element (false) - { - if (size == 0) - return; - - if (std_cxx1x::is_trivial::value == true) - { - const unsigned char zero [sizeof(T)] = {}; - if (std::memcmp(zero, &element, sizeof(T)) == 0) - trivial_element = true; - } - if (size < minimum_parallel_grain_size) - apply_to_subrange (0, size); - else - apply_parallel (0, size, minimum_parallel_grain_size); - } - - private: - - /** - * This sets elements on a subrange given by - * two integers. - */ - virtual void apply_to_subrange (const std::size_t begin, - const std::size_t end) const - { - // for classes with trivial assignment of zero - // can use memset - if (std_cxx1x::is_trivial::value == true && trivial_element) - std::memset (destination_+begin, 0, (end-begin)*sizeof(T)); - else - for (std::size_t i=begin; i -class AlignedVector -{ -public: - /** - * Declare standard types used in all - * containers. These types parallel those - * in the C++ standard libraries - * vector<...> class. - */ - typedef T value_type; - typedef value_type *pointer; - typedef const value_type *const_pointer; - typedef value_type *iterator; - typedef const value_type *const_iterator; - typedef value_type &reference; - typedef const value_type &const_reference; - typedef std::size_t size_type; - - /** - * Empty constructor. Sets the vector size to - * zero. - */ - AlignedVector () - : - _data (0), - _end_data (0), - _end_allocated (0) - {}; - - /** - * Sets the vector size to the given size and - * initializes all elements with T(). - */ - AlignedVector (const size_type size) - : - _data (0), - _end_data (0), - _end_allocated (0) - { - if (size > 0) - resize (size); - } - - /** - * Destructor. - */ - ~AlignedVector () - { - clear(); - } - - /** - * Copy constructor. - */ - AlignedVector (const AlignedVector &vec) - : - _data (0), - _end_data (0), - _end_allocated (0) - { - // do not invalidate old data - resize_fast (vec._end_data - vec._data); - internal::AlignedVectorMove (vec._data, vec._end_data, _data, true); - } - - /** - * Assignment to the input vector @p vec. - */ - AlignedVector & - operator = (const AlignedVector &vec) - { - clear(); - resize_fast (vec._end_data - vec._data); - internal::AlignedVectorMove (vec._data, vec._end_data, _data, true); - return *this; - } - - /** - * Change the size of the vector. It keeps old - * elements previously available but does not - * initialize the newly allocated memory, - * leaving it in an undefined state. - */ - void resize_fast (const size_type size) - { - reserve (size); - _end_data = _data + size; - } - - /** - * Change the size of the vector. It keeps old - * elements previously available, and - * initializes each element with the specified - * data. If the new vector size is shorter - * than the old one, the memory is not - * released unless the new size is zero. - */ - void resize (const size_type size_in, - const T &init = T()) - { - const size_type old_size = size(); - if (std_cxx1x::is_trivial::value == false && size_in < old_size) - { - // call destructor on fields that are released - while (_end_data != _data+size_in) - (--_end_data)->~T(); - } - - resize_fast (size_in); - // now _size is set correctly, need to set the - // values - if (size_in > old_size) - internal::AlignedVectorSet (size_in-old_size, init, - _data+old_size); - } - - /** - * Reserve memory space for @p size - * elements. If the argument @p size is set to - * zero, all previously allocated memory is - * released. - * - * In order to avoid too frequent reallocation - * (which involves copy of the data), this - * function doubles the amount of memory - * occupied when the given size is larger than - * the previously allocated size. - */ - void reserve (const size_type size_alloc) - { - const size_type old_size = _end_data - _data; - const size_type allocated_size = _end_allocated - _data; - if (size_alloc > allocated_size) - { - // if we continuously increase the size of the - // vector, we might be reallocating a lot of - // times. therefore, try to increase the size - // more aggressively - size_type new_size = size_alloc; - if (size_alloc < (2 * allocated_size)) - new_size = 2 * allocated_size; - - const size_type size_actual_allocate = new_size * sizeof(T); - -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - - // allocate and align along boundaries of the - // size of VectorizedArray, which is - // 16 bytes for SSE and 32 bytes for AVX - T *new_data = static_cast(_mm_malloc (size_actual_allocate, - sizeof(VectorizedArray))); -#else - T *new_data = static_cast(malloc (size_actual_allocate)); -#endif - if (new_data == 0) - throw std::bad_alloc(); - - // copy data in case there was some content - // before and release the old memory with the - // function corresponding to the one used for - // allocating - std::swap (_data, new_data); - _end_data = _data + old_size; - _end_allocated = _data + new_size; - if (_end_data != _data) - { - internal::AlignedVectorMove(new_data, new_data + old_size, - _data); -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - _mm_free(new_data); -#else - free(new_data); -#endif - } - } - else if (size_alloc == 0) - clear(); - } - - /** - * Releases all previously allocated memory - * and leaves the vector in a state equivalent - * to the state after the default constructor - * has been called. - */ - void clear () - { - if (_data != 0) - { - if (std_cxx1x::is_trivial::value == false) - while (_end_data != _data) - (--_end_data)->~T(); - -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - _mm_free(_data); -#else - free(_data); -#endif - } - _data = 0; - _end_data = 0; - _end_allocated = 0; - }; - - /** - * Inserts an element at the end of the - * vector, increasing the vector size by - * one. Note that the allocated size will - * double whenever the previous space is not - * enough to hold the new element. - */ - void push_back (const T in_data) - { - Assert (_end_data <= _end_allocated, ExcInternalError()); - if (_end_data == _end_allocated) - reserve (std::max(2*capacity(),static_cast(16))); - if (std_cxx1x::is_trivial::value == false) - new (_end_data) T; - *_end_data++ = in_data; - } - - /** - * Returns the last element of the vector - * (read and write access). - */ - reference back () - { - AssertIndexRange (0, size()); - T *field = _end_data - 1; - return *field; - } - - /** - * Returns the last element of the vector - * (read-only access). - */ - const_reference back () const - { - AssertIndexRange (0, size()); - const T *field = _end_data - 1; - return *field; - } - - /** - * Inserts several elements at the end of the - * vector given by a range of elements. - */ - template - void insert_back (ForwardIterator begin, - ForwardIterator end) - { - const unsigned int old_size = size(); - reserve (old_size + (end-begin)); - for ( ; begin != end; ++begin, ++_end_data) - { - if (std_cxx1x::is_trivial::value == false) - new (_end_data) T; - *_end_data = *begin; - } - } - - /** - * Swaps the given vector with the calling - * vector. - */ - void swap (AlignedVector &vec) - { - std::swap (_data, vec._data); - std::swap (_end_data, vec._end_data); - std::swap (_end_allocated, vec._end_allocated); - } - - /** - * Returns the size of the vector. - */ - size_type size () const - { - return _end_data - _data; - } - - /** - * Returns the capacity of the vector, i.e., - * the size this vector can hold without - * reallocation. Note that capacity() >= - * size(). - */ - size_type capacity () const - { - return _end_allocated - _data; - } - - /** - * Read-write access to entry @p index in the - * vector. - */ - reference - operator [] (const size_type index) - { - AssertIndexRange (index, size()); - return _data[index]; - }; - - /** - * Read-only access to entry @p index in the - * vector. - */ - const_reference operator [] (const size_type index) const - { - AssertIndexRange (index, size()); - return _data[index]; - }; - - /** - * Returns a read and write pointer to the - * beginning of the data array. - */ - iterator begin () - { - return _data; - } - - /** - * Returns a read and write pointer to the - * end of the data array. - */ - iterator end () - { - return _end_data; - } - - /** - * Returns a read-only pointer to the - * beginning of the data array. - */ - const_iterator begin () const - { - return _data; - } - - /** - * Returns a read-only pointer to the - * end of the data array. - */ - const_iterator end () const - { - return _end_data; - } - - /** - * Returns the memory consumption of the - * allocated memory in this class. If the - * underlying type @p T allocates memory by - * itself, this memory is not counted. - */ - size_type memory_consumption () const - { - size_type memory = sizeof(this); - memory += sizeof(T) * capacity(); - return memory; - } - -private: - - /** - * Pointer to actual class data. - */ - T *_data; - - /** - * Pointer to the end of valid data fields. - */ - T *_end_data; - - /** - * Pointer to the end of the allocated memory. - */ - T *_end_allocated; -}; - - DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index c9f3b4434e..c50275614a 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -3897,7 +3897,7 @@ namespace parallel for (i = 0; i < l; i++) { typename dealii::Triangulation::cell_iterator cell (triangulation, i, dealii_index); - child_id = internal::p4est::functions::quadrant_ancestor_id (&quad, i + 1); + child_id = dealii::internal::p4est::functions::quadrant_ancestor_id (&quad, i + 1); Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!")); dealii_index = cell->child_index(child_id); } -- 2.39.5