aligned_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 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 dealii_aligned_vector_h
#define dealii_aligned_vector_h

#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/parallel.h>

// boost::serialization::make_array used to be in array.hpp, but was
// moved to a different file in BOOST 1.64
#include <boost/version.hpp>
#if BOOST_VERSION >= 106400
#  include <boost/serialization/array_wrapper.hpp>
#else
#  include <boost/serialization/array.hpp>
#endif
#include <boost/serialization/split_member.hpp>

#include <cstring>
#include <type_traits>



DEAL_II_NAMESPACE_OPEN


template < class T >
class AlignedVector
{
public:
  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;

  AlignedVector ();

  AlignedVector (const size_type size,
                 const T        &init = T());

  ~AlignedVector ();

  AlignedVector (const AlignedVector<T> &vec);

  AlignedVector (AlignedVector<T> &&vec) noexcept;

  AlignedVector &
  operator = (const AlignedVector<T> &vec);

  AlignedVector &
  operator = (AlignedVector<T> &&vec) noexcept;

  void resize_fast (const size_type size);

  void resize (const size_type size_in);

  void resize (const size_type size_in,
               const T        &init);

  void reserve (const size_type size_alloc);

  void clear ();

  void push_back (const T in_data);

  reference back ();

  const_reference back () const;

  template <typename ForwardIterator>
  void insert_back (ForwardIterator begin,
                    ForwardIterator end);

  void fill ();

  void fill (const T &element);

  void swap (AlignedVector<T> &vec);

  bool empty () const;

  size_type size () const;

  size_type capacity () const;

  reference
  operator [] (const size_type index);

  const_reference operator [] (const size_type index) const;

  iterator begin ();

  iterator end ();

  const_iterator begin () const;

  const_iterator end () const;

  size_type memory_consumption () const;

  template <class Archive>
  void save (Archive &ar, const unsigned int version) const;

  template <class Archive>
  void load (Archive &ar, const unsigned int version);

  BOOST_SERIALIZATION_SPLIT_MEMBER()

private:

  T *_data;

  T *_end_data;

  T *_end_allocated;
};


// ------------------------------- inline functions --------------------------

namespace internal
{
  template <typename T>
  class AlignedVectorCopy : private parallel::ParallelForInteger
  {
    static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
  public:
    AlignedVectorCopy (const T *const source_begin,
                       const T *const source_end,
                       T *const destination)
      :
      source_ (source_begin),
      destination_ (destination)
    {
      Assert (source_end >= source_begin, ExcInternalError());
      Assert (source_end == source_begin || destination != nullptr, 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);
    }

    virtual void apply_to_subrange (const std::size_t begin,
                                    const std::size_t end) const
    {
      if (end == begin)
        return;

      // for classes trivial assignment can use memcpy. cast element to
      // (void*) to silence compiler warning for virtual classes (they will
      // never arrive here because they are non-trivial).

      if (std::is_trivial<T>::value == true)
        std::memcpy ((void *)(destination_+begin), (void *)(source_+begin),
                     (end-begin)*sizeof(T));
      else
        for (std::size_t i=begin; i<end; ++i)
          new (&destination_[i]) T(source_[i]);
    }

  private:
    const T *const source_;
    T *const destination_;
  };


  template <typename T>
  class AlignedVectorMove : private parallel::ParallelForInteger
  {
    static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
  public:
    AlignedVectorMove (T *const source_begin,
                       T *const source_end,
                       T *const destination)
      :
      source_ (source_begin),
      destination_ (destination)
    {
      Assert (source_end >= source_begin, ExcInternalError());
      Assert (source_end == source_begin || destination != nullptr, 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);
    }

    virtual void apply_to_subrange (const std::size_t begin,
                                    const std::size_t end) const
    {
      if (end == begin)
        return;

      // for classes trivial assignment can use memcpy. cast element to
      // (void*) to silence compiler warning for virtual classes (they will
      // never arrive here because they are non-trivial).

      if (std::is_trivial<T>::value == true)
        std::memcpy ((void *)(destination_+begin), (void *)(source_+begin),
                     (end-begin)*sizeof(T));
      else
        for (std::size_t i=begin; i<end; ++i)
          {
            // initialize memory (copy construct by placement new), and
            // destruct the source
            new (&destination_[i]) T(std::move(source_[i]));
            source_[i].~T();
          }
    }

  private:
    T *const source_;
    T *const destination_;
  };


  template <typename T, bool initialize_memory>
  class AlignedVectorSet : private parallel::ParallelForInteger
  {
    static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
  public:
    AlignedVectorSet (const std::size_t size,
                      const T &element,
                      T *const destination)
      :
      element_ (element),
      destination_ (destination),
      trivial_element (false)
    {
      if (size == 0)
        return;
      Assert (destination != nullptr, ExcInternalError());

      // do not use memcmp for long double because on some systems it does not
      // completely fill its memory and may lead to false positives in
      // e.g. valgrind
      if (std::is_trivial<T>::value == true &&
          std::is_same<T,long double>::value == false)
        {
          const unsigned char zero [sizeof(T)] = {};
          // cast element to (void*) to silence compiler warning for virtual
          // classes (they will never arrive here because they are
          // non-trivial).
          if (std::memcmp(zero, (void *)&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);
    }

    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. cast
      // element to (void*) to silence compiler warning for virtual
      // classes (they will never arrive here because they are
      // non-trivial).
      if (std::is_trivial<T>::value == true && trivial_element)
        std::memset ((void *)(destination_+begin), 0, (end-begin)*sizeof(T));
      else
        copy_construct_or_assign(begin, end,
                                 std::integral_constant<bool, initialize_memory>());
    }

  private:
    const T &element_;
    mutable T *destination_;
    bool trivial_element;

    // copy assignment operation
    void copy_construct_or_assign(const std::size_t begin,
                                  const std::size_t end,
                                  std::integral_constant<bool, false>) const
    {
      for (std::size_t i=begin; i<end; ++i)
        destination_[i] = element_;
    }

    // copy constructor (memory initialization)
    void copy_construct_or_assign(const std::size_t begin,
                                  const std::size_t end,
                                  std::integral_constant<bool, true>) const
    {
      for (std::size_t i=begin; i<end; ++i)
        new (&destination_[i]) T(element_);
    }
  };



  template <typename T, bool initialize_memory>
  class AlignedVectorDefaultInitialize : private parallel::ParallelForInteger
  {
    static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
  public:
    AlignedVectorDefaultInitialize (const std::size_t size,
                                    T *const destination)
      :
      destination_ (destination)
    {
      if (size == 0)
        return;
      Assert (destination != nullptr, ExcInternalError());

      if (size < minimum_parallel_grain_size)
        apply_to_subrange (0, size);
      else
        apply_parallel (0, size, minimum_parallel_grain_size);
    }

    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. cast
      // element to (void*) to silence compiler warning for virtual
      // classes (they will never arrive here because they are
      // non-trivial).
      if (std::is_trivial<T>::value == true)
        std::memset ((void *)(destination_+begin), 0, (end-begin)*sizeof(T));
      else
        default_construct_or_assign(begin, end,
                                    std::integral_constant<bool, initialize_memory>());
    }

  private:
    mutable T *destination_;

    // copy assignment operation
    void default_construct_or_assign(const std::size_t begin,
                                     const std::size_t end,
                                     std::integral_constant<bool, false>) const
    {
      for (std::size_t i=begin; i<end; ++i)
        destination_[i] = std::move(T());
    }

    // copy constructor (memory initialization)
    void default_construct_or_assign(const std::size_t begin,
                                     const std::size_t end,
                                     std::integral_constant<bool, true>) const
    {
      for (std::size_t i=begin; i<end; ++i)
        new (&destination_[i]) T;
    }
  };

} // end of namespace internal


#ifndef DOXYGEN


template < class T >
inline
AlignedVector<T>::AlignedVector ()
  :
  _data (nullptr),
  _end_data (nullptr),
  _end_allocated (nullptr)
{}



template < class T >
inline
AlignedVector<T>::AlignedVector (const size_type size,
                                 const T &init)
  :
  _data (nullptr),
  _end_data (nullptr),
  _end_allocated (nullptr)
{
  if (size > 0)
    resize (size, init);
}



template < class T >
inline
AlignedVector<T>::~AlignedVector ()
{
  clear();
}



template < class T >
inline
AlignedVector<T>::AlignedVector (const AlignedVector<T> &vec)
  :
  _data (nullptr),
  _end_data (nullptr),
  _end_allocated (nullptr)
{
  // copy the data from vec
  reserve (vec._end_data - vec._data);
  _end_data = _end_allocated;
  internal::AlignedVectorCopy<T> (vec._data, vec._end_data, _data);
}



template < class T >
inline
AlignedVector<T>::AlignedVector (AlignedVector<T> &&vec) noexcept
:
_data (vec._data),
      _end_data (vec._end_data),
      _end_allocated (vec._end_allocated)
{
  vec._data = nullptr;
  vec._end_data = nullptr;
  vec._end_allocated = nullptr;
}



template < class T >
inline
AlignedVector<T> &
AlignedVector<T>::operator = (const AlignedVector<T> &vec)
{
  resize(0);
  resize_fast (vec._end_data - vec._data);
  internal::AlignedVectorCopy<T> (vec._data, vec._end_data, _data);
  return *this;
}



template < class T >
inline
AlignedVector<T> &
AlignedVector<T>::operator = (AlignedVector<T> &&vec) noexcept
{
  clear();

  _data = vec._data;
  _end_data = vec._end_data;
  _end_allocated = vec._end_allocated;

  vec._data = nullptr;
  vec._end_data = nullptr;
  vec._end_allocated = nullptr;

  return *this;
}



template < class T >
inline
void
AlignedVector<T>::resize_fast (const size_type size_in)
{
  const size_type old_size = size();
  if (std::is_trivial<T>::value == false && size_in < old_size)
    {
      // call destructor on fields that are released. doing it backward
      // releases the elements in reverse order as compared to how they were
      // created
      while (_end_data != _data+size_in)
        (--_end_data)->~T();
    }
  reserve (size_in);
  _end_data = _data + size_in;

  // need to still set the values in case the class is non-trivial because
  // virtual classes etc. need to run their (default) constructor
  if (std::is_trivial<T>::value == false && size_in > old_size)
    ::internal::AlignedVectorDefaultInitialize<T,true> (size_in-old_size, _data+old_size);
}



template < class T >
inline
void
AlignedVector<T>::resize (const size_type size_in)
{
  const size_type old_size = size();
  if (std::is_trivial<T>::value == false && size_in < old_size)
    {
      // call destructor on fields that are released. doing it backward
      // releases the elements in reverse order as compared to how they were
      // created
      while (_end_data != _data+size_in)
        (--_end_data)->~T();
    }
  reserve (size_in);
  _end_data = _data + size_in;

  // finally set the desired init values
  if (size_in > old_size)
    ::internal::AlignedVectorDefaultInitialize<T,true> (size_in-old_size, _data+old_size);
}



template < class T >
inline
void
AlignedVector<T>::resize (const size_type size_in,
                          const T        &init)
{
  const size_type old_size = size();
  if (std::is_trivial<T>::value == false && size_in < old_size)
    {
      // call destructor on fields that are released. doing it backward
      // releases the elements in reverse order as compared to how they were
      // created
      while (_end_data != _data+size_in)
        (--_end_data)->~T();
    }
  reserve (size_in);
  _end_data = _data + size_in;

  // finally set the desired init values
  if (size_in > old_size)
    ::internal::AlignedVectorSet<T,true> (size_in-old_size, init, _data+old_size);
}



template < class T >
inline
void
AlignedVector<T>::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);

      // allocate and align along 64-byte boundaries (this is enough for all
      // levels of vectorization currently supported by deal.II)
      T *new_data;
      Utilities::System::posix_memalign ((void **)&new_data, 64, size_actual_allocate);

      // 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<T>(new_data, new_data + old_size,
                                                 _data);
          free(new_data);
        }
      else
        Assert(new_data == nullptr, ExcInternalError());
    }
  else if (size_alloc == 0)
    clear();
}



template < class T >
inline
void
AlignedVector<T>::clear ()
{
  if (_data != nullptr)
    {
      if (std::is_trivial<T>::value == false)
        while (_end_data != _data)
          (--_end_data)->~T();

      free(_data);
    }
  _data = nullptr;
  _end_data = nullptr;
  _end_allocated = nullptr;
}



template < class T >
inline
void
AlignedVector<T>::push_back (const T in_data)
{
  Assert (_end_data <= _end_allocated, ExcInternalError());
  if (_end_data == _end_allocated)
    reserve (std::max(2*capacity(),static_cast<size_type>(16)));
  if (std::is_trivial<T>::value == false)
    new (_end_data++) T(in_data);
  else
    *_end_data++ = in_data;
}



template < class T >
inline
typename AlignedVector<T>::reference
AlignedVector<T>::back ()
{
  AssertIndexRange (0, size());
  T *field = _end_data - 1;
  return *field;
}



template < class T >
inline
typename AlignedVector<T>::const_reference
AlignedVector<T>::back () const
{
  AssertIndexRange (0, size());
  const T *field = _end_data - 1;
  return *field;
}



template < class T >
template <typename ForwardIterator>
inline
void
AlignedVector<T>::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::is_trivial<T>::value == false)
        new (_end_data) T;
      *_end_data = *begin;
    }
}



template < class T >
inline
void
AlignedVector<T>::fill ()
{
  ::internal::AlignedVectorDefaultInitialize<T,false> (size(), _data);
}



template < class T >
inline
void
AlignedVector<T>::fill (const T &value)
{
  ::internal::AlignedVectorSet<T,false> (size(), value, _data);
}



template < class T >
inline
void
AlignedVector<T>::swap (AlignedVector<T> &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<T>::empty () const
{
  return _end_data == _data;
}



template < class T >
inline
typename AlignedVector<T>::size_type
AlignedVector<T>::size () const
{
  return _end_data - _data;
}



template < class T >
inline
typename AlignedVector<T>::size_type
AlignedVector<T>::capacity () const
{
  return _end_allocated - _data;
}



template < class T >
inline
typename AlignedVector<T>::reference
AlignedVector<T>::operator [] (const size_type index)
{
  AssertIndexRange (index, size());
  return _data[index];
}



template < class T >
inline
typename AlignedVector<T>::const_reference
AlignedVector<T>::operator [] (const size_type index) const
{
  AssertIndexRange (index, size());
  return _data[index];
}



template < class T >
inline
typename AlignedVector<T>::iterator
AlignedVector<T>::begin ()
{
  return _data;
}



template < class T >
inline
typename AlignedVector<T>::iterator
AlignedVector<T>::end ()
{
  return _end_data;
}



template < class T >
inline
typename AlignedVector<T>::const_iterator
AlignedVector<T>::begin () const
{
  return _data;
}



template < class T >
inline
typename AlignedVector<T>::const_iterator
AlignedVector<T>::end () const
{
  return _end_data;
}



template < class T >
template < class Archive >
inline
void
AlignedVector<T>::save (Archive &ar, const unsigned int) const
{
  size_type vec_size (size());
  ar &vec_size;
  if (vec_size > 0)
    ar &boost::serialization::make_array(_data, vec_size);
}



template < class T >
template < class Archive >
inline
void
AlignedVector<T>::load (Archive &ar, const unsigned int)
{
  size_type vec_size = 0;
  ar &vec_size ;

  if (vec_size > 0)
    {
      reserve(vec_size);
      ar &boost::serialization::make_array(_data, vec_size);
      _end_data = _data + vec_size;
    }
}



template < class T >
inline
typename AlignedVector<T>::size_type
AlignedVector<T>::memory_consumption () const
{
  size_type memory = sizeof(*this);
  for (const T *t = _data ; t != _end_data; ++t)
    memory += ::MemoryConsumption::memory_consumption(*t);
  memory += sizeof(T) * (_end_allocated-_end_data);
  return memory;
}


#endif // ifndef DOXYGEN


template < class T >
bool operator == (const AlignedVector<T> &lhs,
                  const AlignedVector<T> &rhs)
{
  if (lhs.size() != rhs.size())
    return false;
  for (typename AlignedVector<T>::const_iterator lit = lhs.begin(),
       rit = rhs.begin(); lit != lhs.end(); ++lit, ++rit)
    if (*lit != *rit)
      return false;
  return true;
}




template < class T >
bool operator != (const AlignedVector<T> &lhs,
                  const AlignedVector<T> &rhs)
{
  return !(operator==(lhs, rhs));
}


DEAL_II_NAMESPACE_CLOSE

#endif