vector.h

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

#ifndef dealii_vector_h
#define dealii_vector_h


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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/differentiation/ad/ad_number_traits.h>

#include <deal.II/lac/vector_operation.h>
#include <deal.II/lac/vector_type_traits.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 <cstdio>
#include <cstring>
#include <iostream>
#include <vector>

DEAL_II_NAMESPACE_OPEN


#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  class VectorBase;
}
#endif

#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace TrilinosWrappers
#endif

template <typename number>
class LAPACKFullMatrix;

template <typename>
class BlockVector;

template <typename>
class VectorView;

namespace parallel
{
  namespace internal
  {
    class TBBPartitioner;
  }
} // namespace parallel



template <typename Number>
class Vector : public Subscriptor
{
public:
  // The assertion in vector.templates.h for whether or not a number is
  // finite is not compatible for AD number types.
  static_assert(
    !Differentiation::AD::is_ad_number<Number>::value,
    "The Vector class does not support auto-differentiable numbers.");

  using value_type      = Number;
  using pointer         = value_type *;
  using const_pointer   = const value_type *;
  using iterator        = value_type *;
  using const_iterator  = const value_type *;
  using reference       = value_type &;
  using const_reference = const value_type &;
  using size_type       = types::global_dof_index;

  using real_type = typename numbers::NumberTraits<Number>::real_type;

public:
  Vector();

  Vector(const Vector<Number> &v);

  Vector(Vector<Number> &&v) noexcept;

  template <typename OtherNumber>
  explicit Vector(const Vector<OtherNumber> &v);

#ifdef DEAL_II_WITH_PETSC

  explicit Vector(const PETScWrappers::VectorBase &v);
#endif

#ifdef DEAL_II_WITH_TRILINOS

  explicit Vector(const TrilinosWrappers::MPI::Vector &v);
#endif

  explicit Vector(const size_type n);

  template <typename InputIterator>
  Vector(const InputIterator first, const InputIterator last);

  virtual ~Vector() override = default;

  void
  compress(::VectorOperation::values operation =
             ::VectorOperation::unknown) const;

  virtual void
  reinit(const size_type N, const bool omit_zeroing_entries = false);

  void
  grow_or_shrink(const size_type N);

  template <typename Number2>
  void
  reinit(const Vector<Number2> &V, const bool omit_zeroing_entries = false);

  virtual void
  swap(Vector<Number> &v);

  Vector<Number> &
  operator=(const Number s);

  Vector<Number> &
  operator=(const Vector<Number> &v);

  Vector<Number> &
  operator=(Vector<Number> &&v) noexcept;

  template <typename Number2>
  Vector<Number> &
  operator=(const Vector<Number2> &v);

  Vector<Number> &
  operator=(const BlockVector<Number> &v);

#ifdef DEAL_II_WITH_PETSC

  Vector<Number> &
  operator=(const PETScWrappers::VectorBase &v);
#endif


#ifdef DEAL_II_WITH_TRILINOS

  Vector<Number> &
  operator=(const TrilinosWrappers::MPI::Vector &v);
#endif

  template <typename Number2>
  bool
  operator==(const Vector<Number2> &v) const;

  template <typename Number2>
  bool
  operator!=(const Vector<Number2> &v) const;




  template <typename Number2>
  Number operator*(const Vector<Number2> &V) const;

  real_type
  norm_sqr() const;

  Number
  mean_value() const;

  real_type
  l1_norm() const;

  real_type
  l2_norm() const;

  real_type
  lp_norm(const real_type p) const;

  real_type
  linfty_norm() const;

  Number
  add_and_dot(const Number a, const Vector<Number> &V, const Vector<Number> &W);




  iterator
  begin();

  const_iterator
  begin() const;

  iterator
  end();

  const_iterator
  end() const;

  Number
  operator()(const size_type i) const;

  Number &
  operator()(const size_type i);

  Number operator[](const size_type i) const;

  Number &operator[](const size_type i);

  template <typename OtherNumber>
  void
  extract_subvector_to(const std::vector<size_type> &indices,
                       std::vector<OtherNumber> &    values) const;

  template <typename ForwardIterator, typename OutputIterator>
  void
  extract_subvector_to(ForwardIterator       indices_begin,
                       const ForwardIterator indices_end,
                       OutputIterator        values_begin) const;



  Vector<Number> &
  operator+=(const Vector<Number> &V);

  Vector<Number> &
  operator-=(const Vector<Number> &V);

  template <typename OtherNumber>
  void
  add(const std::vector<size_type> &  indices,
      const std::vector<OtherNumber> &values);

  template <typename OtherNumber>
  void
  add(const std::vector<size_type> &indices, const Vector<OtherNumber> &values);

  template <typename OtherNumber>
  void
  add(const size_type    n_elements,
      const size_type *  indices,
      const OtherNumber *values);

  void
  add(const Number s);

  void
  add(const Number          a,
      const Vector<Number> &V,
      const Number          b,
      const Vector<Number> &W);

  void
  add(const Number a, const Vector<Number> &V);

  void
  sadd(const Number s, const Vector<Number> &V);

  void
  sadd(const Number s, const Number a, const Vector<Number> &V);

  Vector<Number> &
  operator*=(const Number factor);

  Vector<Number> &
  operator/=(const Number factor);

  void
  scale(const Vector<Number> &scaling_factors);

  template <typename Number2>
  void
  scale(const Vector<Number2> &scaling_factors);

  void
  equ(const Number a, const Vector<Number> &u);

  template <typename Number2>
  void
  equ(const Number a, const Vector<Number2> &u);

  DEAL_II_DEPRECATED
  void
  ratio(const Vector<Number> &a, const Vector<Number> &b);

  void
  update_ghost_values() const;


  DEAL_II_DEPRECATED
  void
  print(const char *format = nullptr) const;

  void
  print(std::ostream &     out,
        const unsigned int precision  = 3,
        const bool         scientific = true,
        const bool         across     = true) const;

  DEAL_II_DEPRECATED
  void
  print(LogStream &        out,
        const unsigned int width  = 6,
        const bool         across = true) const;

  void
  block_write(std::ostream &out) const;

  void
  block_read(std::istream &in);

  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()

  

  bool
  in_local_range(const size_type global_index) const;

  IndexSet
  locally_owned_elements() const;

  size_type
  size() const;

  bool
  all_zero() const;

  bool
  is_non_negative() const;

  std::size_t
  memory_consumption() const;

protected:
  size_type vec_size;

  size_type max_vec_size;

  std::unique_ptr<Number[], decltype(&free)> values;

  mutable std::shared_ptr<parallel::internal::TBBPartitioner>
    thread_loop_partitioner;

  template <typename Number2>
  friend class Vector;

  template <typename Number2>
  friend class LAPACKFullMatrix;

  friend class VectorView<Number>;

private:
  void
  allocate(const size_type copy_n_el = 0);
};

/*----------------------- Inline functions ----------------------------------*/


#ifndef DOXYGEN


//------------------------ declarations for explicit specializations
template <>
Vector<int>::real_type
Vector<int>::lp_norm(const real_type) const;


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

template <typename Number>
inline Vector<Number>::Vector()
  : vec_size(0)
  , max_vec_size(0)
  , values(nullptr, &free)
{
  // virtual functions called in constructors and destructors never use the
  // override in a derived class
  // for clarity be explicit on which function is called
  Vector<Number>::reinit(0);
}



template <typename Number>
template <typename InputIterator>
Vector<Number>::Vector(const InputIterator first, const InputIterator last)
  : vec_size(0)
  , max_vec_size(0)
  , values(nullptr, &free)
{
  // allocate memory. do not initialize it, as we will copy over to it in a
  // second
  reinit(std::distance(first, last), true);
  std::copy(first, last, begin());
}



template <typename Number>
inline Vector<Number>::Vector(const size_type n)
  : vec_size(0)
  , max_vec_size(0)
  , values(nullptr, &free)
{
  // virtual functions called in constructors and destructors never use the
  // override in a derived class
  // for clarity be explicit on which function is called
  Vector<Number>::reinit(n, false);
}



template <typename Number>
inline typename Vector<Number>::size_type
Vector<Number>::size() const
{
  return vec_size;
}


template <typename Number>
inline bool
Vector<Number>::in_local_range(const size_type) const
{
  return true;
}



template <typename Number>
inline typename Vector<Number>::iterator
Vector<Number>::begin()
{
  return values.get();
}



template <typename Number>
inline typename Vector<Number>::const_iterator
Vector<Number>::begin() const
{
  return values.get();
}



template <typename Number>
inline typename Vector<Number>::iterator
Vector<Number>::end()
{
  return values.get() + vec_size;
}



template <typename Number>
inline typename Vector<Number>::const_iterator
Vector<Number>::end() const
{
  return values.get() + vec_size;
}



template <typename Number>
inline Number
Vector<Number>::operator()(const size_type i) const
{
  Assert(i < vec_size, ExcIndexRange(i, 0, vec_size));
  return values[i];
}



template <typename Number>
inline Number &
Vector<Number>::operator()(const size_type i)
{
  Assert(i < vec_size, ExcIndexRangeType<size_type>(i, 0, vec_size));
  return values[i];
}



template <typename Number>
inline Number Vector<Number>::operator[](const size_type i) const
{
  return operator()(i);
}



template <typename Number>
inline Number &Vector<Number>::operator[](const size_type i)
{
  return operator()(i);
}



template <typename Number>
template <typename OtherNumber>
inline void
Vector<Number>::extract_subvector_to(const std::vector<size_type> &indices,
                                     std::vector<OtherNumber> &    values) const
{
  for (size_type i = 0; i < indices.size(); ++i)
    values[i] = operator()(indices[i]);
}



template <typename Number>
template <typename ForwardIterator, typename OutputIterator>
inline void
Vector<Number>::extract_subvector_to(ForwardIterator       indices_begin,
                                     const ForwardIterator indices_end,
                                     OutputIterator        values_begin) const
{
  while (indices_begin != indices_end)
    {
      *values_begin = operator()(*indices_begin);
      indices_begin++;
      values_begin++;
    }
}



template <typename Number>
inline Vector<Number> &
Vector<Number>::operator/=(const Number factor)
{
  AssertIsFinite(factor);
  Assert(factor != Number(0.), ExcZero());

  this->operator*=(Number(1.) / factor);
  return *this;
}



template <typename Number>
template <typename OtherNumber>
inline void
Vector<Number>::add(const std::vector<size_type> &  indices,
                    const std::vector<OtherNumber> &values)
{
  Assert(indices.size() == values.size(),
         ExcDimensionMismatch(indices.size(), values.size()));
  add(indices.size(), indices.data(), values.data());
}



template <typename Number>
template <typename OtherNumber>
inline void
Vector<Number>::add(const std::vector<size_type> &indices,
                    const Vector<OtherNumber> &   values)
{
  Assert(indices.size() == values.size(),
         ExcDimensionMismatch(indices.size(), values.size()));
  add(indices.size(), indices.data(), values.values.get());
}



template <typename Number>
template <typename OtherNumber>
inline void
Vector<Number>::add(const size_type    n_indices,
                    const size_type *  indices,
                    const OtherNumber *values)
{
  for (size_type i = 0; i < n_indices; ++i)
    {
      Assert(indices[i] < vec_size, ExcIndexRange(indices[i], 0, vec_size));
      Assert(
        numbers::is_finite(values[i]),
        ExcMessage(
          "The given value is not finite but either infinite or Not A Number (NaN)"));

      this->values[indices[i]] += values[i];
    }
}



template <typename Number>
template <typename Number2>
inline bool
Vector<Number>::operator!=(const Vector<Number2> &v) const
{
  return !(*this == v);
}



template <typename Number>
inline void Vector<Number>::compress(::VectorOperation::values) const
{}



template <typename Number>
inline void
Vector<Number>::update_ghost_values() const
{}



// Moved from vector.templates.h as an inline function by Luca Heltai
// on 2009/04/12 to prevent strange compiling errors, after making
// swap virtual.
template <typename Number>
inline void
Vector<Number>::swap(Vector<Number> &v)
{
  std::swap(vec_size, v.vec_size);
  std::swap(max_vec_size, v.max_vec_size);
  std::swap(values, v.values);
}



template <typename Number>
template <class Archive>
inline void
Vector<Number>::save(Archive &ar, const unsigned int) const
{
  // forward to serialization function in the base class.
  ar &static_cast<const Subscriptor &>(*this);

  ar &vec_size &max_vec_size;
  ar &          boost::serialization::make_array(values.get(), max_vec_size);
}



template <typename Number>
template <class Archive>
inline void
Vector<Number>::load(Archive &ar, const unsigned int)
{
  // get rid of previous content
  values.reset();

  // the load stuff again from the archive
  ar &static_cast<Subscriptor &>(*this);
  ar &vec_size &max_vec_size;

  allocate();
  ar &boost::serialization::make_array(values.get(), max_vec_size);
}

#endif


template <typename Number>
inline void
swap(Vector<Number> &u, Vector<Number> &v)
{
  u.swap(v);
}


template <typename number>
inline std::ostream &
operator<<(std::ostream &os, const Vector<number> &v)
{
  v.print(os);
  return os;
}

template <typename number>
inline LogStream &
operator<<(LogStream &os, const Vector<number> &v)
{
  v.print(os);
  return os;
}

template <typename Number>
struct is_serial_vector<Vector<Number>> : std::true_type
{};


DEAL_II_NAMESPACE_CLOSE

#endif