array_view.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2019 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_array_view_h
#define dealii_array_view_h

#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_space.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/table.h>
#include <deal.II/base/tensor.h>

#include <deal.II/lac/lapack_full_matrix.h>

#include <type_traits>
#include <vector>

DEAL_II_NAMESPACE_OPEN


template <typename ElementType, typename MemorySpaceType = MemorySpace::Host>
class ArrayView
{
public:
  using value_type = ElementType;

  using iterator = value_type *;

  using const_iterator = const ElementType *;

  ArrayView();

  ArrayView(value_type *starting_element, const std::size_t n_elements);

  ArrayView(const ArrayView<typename std::remove_cv<value_type>::type,
                            MemorySpaceType> &view);

  ArrayView(
    const std::vector<typename std::remove_cv<value_type>::type> &vector);

  ArrayView(std::vector<typename std::remove_cv<value_type>::type> &vector);

  void
  reinit(value_type *starting_element, const std::size_t n_elements);

  bool
  operator==(
    const ArrayView<const value_type, MemorySpaceType> &other_view) const;

  bool
  operator==(const ArrayView<typename std::remove_cv<value_type>::type,
                             MemorySpaceType> &other_view) const;

  bool
  operator!=(
    const ArrayView<const value_type, MemorySpaceType> &other_view) const;

  bool
  operator!=(const ArrayView<typename std::remove_cv<value_type>::type,
                             MemorySpaceType> &other_view) const;

  std::size_t
  size() const;

  value_type *
  data() const noexcept;

  iterator
  begin() const;

  iterator
  end() const;

  const_iterator
  cbegin() const;

  const_iterator
  cend() const;

  value_type &operator[](const std::size_t i) const;

private:
  value_type *starting_element;

  std::size_t n_elements;

  friend class ArrayView<const ElementType, MemorySpaceType>;
};



//---------------------------------------------------------------------------


namespace internal
{
  namespace ArrayViewHelper
  {
    template <typename MemorySpaceType>
    inline bool
    is_in_correct_memory_space(const void *const ptr)
    {
#ifndef DEAL_II_COMPILER_CUDA_AWARE
      (void)ptr;
      static_assert(std::is_same<MemorySpaceType, MemorySpace::Host>::value,
                    "If the compiler doesn't understand CUDA code, "
                    "the only possible memory space is 'MemorySpace::Host'!");
      return true;
#else
      cudaPointerAttributes attributes;
      const cudaError_t cuda_error = cudaPointerGetAttributes(&attributes, ptr);
      if (cuda_error != cudaErrorInvalidValue)
        {
          AssertCuda(cuda_error);
          if (std::is_same<MemorySpaceType, MemorySpace::Host>::value)
            return attributes.memoryType == cudaMemoryTypeHost;
          else
            return attributes.memoryType == cudaMemoryTypeDevice;
        }
      else
        {
          // ignore and reset the error since host pointers produce an error
          cudaGetLastError();
          return std::is_same<MemorySpaceType, MemorySpace::Host>::value;
        }
#endif
    }
  } // namespace ArrayViewHelper
} // namespace internal



template <typename ElementType, typename MemorySpaceType>
inline ArrayView<ElementType, MemorySpaceType>::ArrayView()
  : starting_element(nullptr)
  , n_elements(0)
{}



template <typename ElementType, typename MemorySpaceType>
inline ArrayView<ElementType, MemorySpaceType>::ArrayView(
  value_type *      starting_element,
  const std::size_t n_elements)
  : starting_element(starting_element)
  , n_elements(n_elements)
{
  Assert(
    n_elements == 0 ||
      internal::ArrayViewHelper::is_in_correct_memory_space<MemorySpaceType>(
        starting_element),
    ExcMessage("The memory space indicated by the template parameter "
               "and the one derived from the pointer value do not match!"));
}



template <typename ElementType, typename MemorySpaceType>
inline void
ArrayView<ElementType, MemorySpaceType>::reinit(value_type *starting_element,
                                                const std::size_t n_elements)
{
  *this = ArrayView(starting_element, n_elements);
}



template <typename ElementType, typename MemorySpaceType>
inline ArrayView<ElementType, MemorySpaceType>::ArrayView(
  const ArrayView<typename std::remove_cv<value_type>::type, MemorySpaceType>
    &view)
  : starting_element(view.starting_element)
  , n_elements(view.n_elements)
{}



template <typename ElementType, typename MemorySpaceType>
inline ArrayView<ElementType, MemorySpaceType>::ArrayView(
  const std::vector<typename std::remove_cv<value_type>::type> &vector)
  : // use delegating constructor
  ArrayView(vector.data(), vector.size())
{
  // the following static_assert is not strictly necessary because,
  // if we got a const std::vector reference argument but ElementType
  // is not itself const, then the call to the forwarding constructor
  // above will already have failed: vector.data() will have returned
  // a const pointer, but we need a non-const pointer.
  //
  // nevertheless, leave the static_assert in since it provides a
  // more descriptive error message that will simply come after the first
  // error produced above
  static_assert(std::is_const<value_type>::value == true,
                "This constructor may only be called if the ArrayView "
                "object has a const value_type. In other words, you can "
                "only create an ArrayView to const values from a const "
                "std::vector.");
}



template <typename ElementType, typename MemorySpaceType>
inline ArrayView<ElementType, MemorySpaceType>::ArrayView(
  std::vector<typename std::remove_cv<value_type>::type> &vector)
  : // use delegating constructor
  ArrayView(vector.data(), vector.size())
{}



template <typename ElementType, typename MemorySpaceType>
inline bool
ArrayView<ElementType, MemorySpaceType>::
operator==(const ArrayView<const value_type, MemorySpaceType> &other_view) const
{
  return (other_view.data() == starting_element) &&
         (other_view.size() == n_elements);
}



template <typename ElementType, typename MemorySpaceType>
inline bool
ArrayView<ElementType, MemorySpaceType>::
operator==(const ArrayView<typename std::remove_cv<value_type>::type,
                           MemorySpaceType> &other_view) const
{
  return (other_view.data() == starting_element) &&
         (other_view.size() == n_elements);
}



template <typename ElementType, typename MemorySpaceType>
inline bool
ArrayView<ElementType, MemorySpaceType>::
operator!=(const ArrayView<const value_type, MemorySpaceType> &other_view) const
{
  return !(*this == other_view);
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::value_type *
ArrayView<ElementType, MemorySpaceType>::data() const noexcept
{
  if (n_elements == 0)
    return nullptr;
  else
    return starting_element;
}



template <typename ElementType, typename MemorySpaceType>
inline bool
ArrayView<ElementType, MemorySpaceType>::
operator!=(const ArrayView<typename std::remove_cv<value_type>::type,
                           MemorySpaceType> &other_view) const
{
  return !(*this == other_view);
}



template <typename ElementType, typename MemorySpaceType>
inline std::size_t
ArrayView<ElementType, MemorySpaceType>::size() const
{
  return n_elements;
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::iterator
ArrayView<ElementType, MemorySpaceType>::begin() const
{
  return starting_element;
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::iterator
ArrayView<ElementType, MemorySpaceType>::end() const
{
  return starting_element + n_elements;
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::const_iterator
ArrayView<ElementType, MemorySpaceType>::cbegin() const
{
  return starting_element;
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::const_iterator
ArrayView<ElementType, MemorySpaceType>::cend() const
{
  return starting_element + n_elements;
}



template <typename ElementType, typename MemorySpaceType>
inline typename ArrayView<ElementType, MemorySpaceType>::value_type &
  ArrayView<ElementType, MemorySpaceType>::operator[](const std::size_t i) const
{
  Assert(i < n_elements, ExcIndexRange(i, 0, n_elements));
  Assert(
    (std::is_same<MemorySpaceType, MemorySpace::Host>::value),
    ExcMessage(
      "Accessing elements is only allowed if the data is stored in CPU memory!"));

  return *(starting_element + i);
}



#ifndef DOXYGEN
namespace internal
{
  namespace ArrayViewHelper
  {
    template <class Iterator>
    bool
    is_contiguous(const Iterator &first, const Iterator &last)
    {
      const auto n = std::distance(first, last);
      for (typename std::decay<decltype(n)>::type i = 0; i < n; ++i)
        if (std::addressof(*(std::next(first, i))) !=
            std::next(std::addressof(*first), i))
          return false;
      return true;
    }


    template <class T>
    constexpr bool
    is_contiguous(T *, T *)
    {
      return true;
    }
  } // namespace ArrayViewHelper
} // namespace internal
#endif



template <typename Iterator, typename MemorySpaceType = MemorySpace::Host>
ArrayView<typename std::remove_reference<
            typename std::iterator_traits<Iterator>::reference>::type,
          MemorySpaceType>
make_array_view(const Iterator begin, const Iterator end)
{
  static_assert(
    std::is_same<typename std::iterator_traits<Iterator>::iterator_category,
                 typename std::random_access_iterator_tag>::value,
    "The provided iterator should be a random access iterator.");
  Assert(begin <= end,
         ExcMessage(
           "The beginning of the array view should be before the end."));
  Assert(internal::ArrayViewHelper::is_contiguous(begin, end),
         ExcMessage("The provided range isn't contiguous in memory!"));
  // the reference type, not the value type, knows the constness of the iterator
  return ArrayView<typename std::remove_reference<
                     typename std::iterator_traits<Iterator>::reference>::type,
                   MemorySpaceType>(std::addressof(*begin), end - begin);
}



template <typename ElementType, typename MemorySpaceType = MemorySpace::Host>
ArrayView<ElementType, MemorySpaceType>
make_array_view(ElementType *const begin, ElementType *const end)
{
  Assert(begin <= end,
         ExcMessage(
           "The beginning of the array view should be before the end."));
  return ArrayView<ElementType, MemorySpaceType>(begin, end - begin);
}



template <typename Number, typename MemorySpaceType>
inline ArrayView<const Number, MemorySpaceType>
make_array_view(const ArrayView<Number, MemorySpaceType> &array_view)
{
  return make_array_view(array_view.cbegin(), array_view.cend());
}



template <typename Number, typename MemorySpaceType>
inline ArrayView<Number, MemorySpaceType>
make_array_view(ArrayView<Number, MemorySpaceType> &array_view)
{
  return make_array_view(array_view.begin(), array_view.end());
}



template <int rank, int dim, typename Number>
inline ArrayView<const Number>
make_array_view(const Tensor<rank, dim, Number> &tensor)
{
  return make_array_view(tensor.begin_raw(), tensor.end_raw());
}



template <int rank, int dim, typename Number>
inline ArrayView<Number>
make_array_view(Tensor<rank, dim, Number> &tensor)
{
  return make_array_view(tensor.begin_raw(), tensor.end_raw());
}



template <int rank, int dim, typename Number>
inline ArrayView<const Number>
make_array_view(const SymmetricTensor<rank, dim, Number> &tensor)
{
  return make_array_view(tensor.begin_raw(), tensor.end_raw());
}



template <int rank, int dim, typename Number>
inline ArrayView<Number>
make_array_view(SymmetricTensor<rank, dim, Number> &tensor)
{
  return make_array_view(tensor.begin_raw(), tensor.end_raw());
}



template <typename ElementType, int N>
inline ArrayView<ElementType> make_array_view(ElementType (&array)[N])
{
  return ArrayView<ElementType>(array, N);
}



template <typename ElementType>
inline ArrayView<ElementType>
make_array_view(Vector<ElementType> &vector)
{
  return ArrayView<ElementType>(vector.begin(), vector.size());
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const Vector<ElementType> &vector)
{
  return ArrayView<const ElementType>(vector.begin(), vector.size());
}



template <typename ElementType>
inline ArrayView<ElementType>
make_array_view(std::vector<ElementType> &vector)
{
  return ArrayView<ElementType>(vector.data(), vector.size());
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const std::vector<ElementType> &vector)
{
  return ArrayView<const ElementType>(vector.data(), vector.size());
}



template <typename ElementType>
inline ArrayView<ElementType>
make_array_view(std::vector<ElementType> &vector,
                const std::size_t         starting_index,
                const std::size_t         size_of_view)
{
  Assert(starting_index + size_of_view <= vector.size(),
         ExcMessage("The starting index and size of the view you want to "
                    "create would lead to a view that extends beyond the end "
                    "of the given vector."));
  return ArrayView<ElementType>(&vector[starting_index], size_of_view);
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const std::vector<ElementType> &vector,
                const std::size_t               starting_index,
                const std::size_t               size_of_view)
{
  Assert(starting_index + size_of_view <= vector.size(),
         ExcMessage("The starting index and size of the view you want to "
                    "create would lead to a view that extends beyond the end "
                    "of the given vector."));
  return ArrayView<const ElementType>(&vector[starting_index], size_of_view);
}



template <typename ElementType>
inline ArrayView<ElementType>
  make_array_view(Table<2, ElementType> &                         table,
                  const typename Table<2, ElementType>::size_type row)
{
  AssertIndexRange(row, table.size()[0]);
  return ArrayView<ElementType>(&table[row][0], table.size()[1]);
}



template <typename ElementType>
inline ArrayView<ElementType> make_array_view(Table<2, ElementType> &table)
{
  return ArrayView<ElementType>(&table[0][0], table.n_elements());
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const Table<2, ElementType> &table)
{
  return ArrayView<const ElementType>(&table[0][0], table.n_elements());
}


template <typename ElementType>
inline ArrayView<ElementType>
make_array_view(LAPACKFullMatrix<ElementType> &matrix)
{
  return ArrayView<ElementType>(&matrix(0, 0), matrix.n_elements());
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const LAPACKFullMatrix<ElementType> &matrix)
{
  return ArrayView<const ElementType>(&matrix(0, 0), matrix.n_elements());
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const Table<2, ElementType> &                   table,
                const typename Table<2, ElementType>::size_type row)
{
  AssertIndexRange(row, table.size()[0]);
  return ArrayView<const ElementType>(&table[row][0], table.size()[1]);
}



template <typename ElementType>
inline ArrayView<ElementType> make_array_view(
  Table<2, ElementType> &                         table,
  const typename Table<2, ElementType>::size_type row,
  const typename Table<2, ElementType>::size_type starting_column,
  const std::size_t                               size_of_view)
{
  AssertIndexRange(row, table.size()[0]);
  AssertIndexRange(starting_column, table.size()[1]);
  Assert(starting_column + size_of_view <= table.size()[1],
         ExcMessage("The starting index and size of the view you want to "
                    "create would lead to a view that extends beyond the end "
                    "of a column of the given table."));
  return ArrayView<ElementType>(&table[row][starting_column], size_of_view);
}



template <typename ElementType>
inline ArrayView<const ElementType>
make_array_view(const Table<2, ElementType> &                   table,
                const typename Table<2, ElementType>::size_type row,
                const typename Table<2, ElementType>::size_type starting_column,
                const std::size_t                               size_of_view)
{
  AssertIndexRange(row, table.size()[0]);
  AssertIndexRange(starting_column, table.size()[1]);
  Assert(starting_column + size_of_view <= table.size()[1],
         ExcMessage("The starting index and size of the view you want to "
                    "create would lead to a view that extends beyond the end "
                    "of a column of the given table."));
  return ArrayView<const ElementType>(&table[row][starting_column],
                                      size_of_view);
}



/*
 * Create a view that doesn't allow the container it points to to be modified.
 * This is useful if the object passed in is not `const` already and a function
 * requires a view to constant memory in its signature.
 *
 * This function returns an object of type `ArrayView<const T>` where `T` is the
 * element type of the container.
 *
 * @relatesalso ArrayView
 */
template <typename Container>
inline auto
make_const_array_view(const Container &container)
  -> decltype(make_array_view(container))
{
  return make_array_view(container);
}


DEAL_II_NAMESPACE_CLOSE

#endif