read_write_vector.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2015 - 2021 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_read_write_vector_h
#define dealii_read_write_vector_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/communication_pattern_base.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/parallel.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/lac/vector_operation.h>
 
#include <cstdlib>
#include <cstring>
#include <iomanip>
 
#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
#  include <deal.II/lac/trilinos_epetra_vector.h>
#  include <deal.II/lac/trilinos_tpetra_vector.h>
 
#  include <Epetra_MultiVector.h>
 
#endif
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace LinearAlgebra
{
  template <typename>
  class Vector;
  namespace distributed
  {
    template <typename, typename>
    class Vector;
  } // namespace distributed
} // namespace LinearAlgebra
 
#  ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace PETScWrappers
#  endif
 
#  ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace TrilinosWrappers
#  endif
 
#  ifdef DEAL_II_WITH_CUDA
namespace LinearAlgebra
{
  namespace CUDAWrappers
  {
    template <typename>
    class Vector;
  }
} // namespace LinearAlgebra
#  endif
#endif
 
namespace LinearAlgebra
{
  template <typename Number>
  class ReadWriteVector : public Subscriptor
  {
  public:
    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;
 
    ReadWriteVector();
 
    ReadWriteVector(const ReadWriteVector<Number> &in_vector);
 
    explicit ReadWriteVector(const size_type size);
 
    explicit ReadWriteVector(const IndexSet &locally_stored_indices);
 
    ~ReadWriteVector() override = default;
 
    virtual void
    reinit(const size_type size, const bool omit_zeroing_entries = false);
 
    template <typename Number2>
    void
    reinit(const ReadWriteVector<Number2> &in_vector,
           const bool                      omit_zeroing_entries = false);
 
    virtual void
    reinit(const IndexSet &locally_stored_indices,
           const bool      omit_zeroing_entries = false);
 
 
#ifdef DEAL_II_WITH_TRILINOS
#  ifdef DEAL_II_WITH_MPI
    void
    reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
#  endif
#endif
 
    template <typename Functor>
    void
    apply(const Functor &func);
 
    void
    swap(ReadWriteVector<Number> &v);
 
    ReadWriteVector<Number> &
    operator=(const ReadWriteVector<Number> &in_vector);
 
    template <typename Number2>
    ReadWriteVector<Number> &
    operator=(const ReadWriteVector<Number2> &in_vector);
 
    ReadWriteVector<Number> &
    operator=(const Number s);
 
    void
    import(const ::Vector<Number> &vec,
           VectorOperation::values       operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
 
    void
    import(const LinearAlgebra::Vector<Number> &vec,
           VectorOperation::values              operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
 
    template <typename MemorySpace>
    void
    import(const distributed::Vector<Number, MemorySpace> &vec,
           VectorOperation::values                         operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
 
#ifdef DEAL_II_WITH_PETSC
    void
    import(const PETScWrappers::MPI::Vector &petsc_vec,
           VectorOperation::values           operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
#endif
 
#ifdef DEAL_II_WITH_TRILINOS
    void
    import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
           VectorOperation::values              operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
 
#  ifdef DEAL_II_WITH_MPI
#    ifdef DEAL_II_TRILINOS_WITH_TPETRA
    void
    import(const TpetraWrappers::Vector<Number> &tpetra_vec,
           VectorOperation::values               operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
#    endif
 
    void
    import(const EpetraWrappers::Vector &epetra_vec,
           VectorOperation::values       operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
#  endif
#endif
 
#ifdef DEAL_II_WITH_CUDA
    void
    import(const CUDAWrappers::Vector<Number> &cuda_vec,
           VectorOperation::values             operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {});
#endif
 
    size_type
    size() const;
 
    DEAL_II_DEPRECATED_EARLY
    size_type
    n_elements() const;
 
    size_type
    locally_owned_size() const;
 
    const IndexSet &
    get_stored_elements() const;
 
    iterator
    begin();
 
    const_iterator
    begin() const;
 
    iterator
    end();
 
    const_iterator
    end() const;
 
 
 
    Number
    operator()(const size_type global_index) const;
 
    Number &
    operator()(const size_type global_index);
 
    Number operator[](const size_type global_index) const;
 
    Number &operator[](const size_type global_index);
 
    template <typename Number2>
    void
    extract_subvector_to(const std::vector<size_type> &indices,
                         std::vector<Number2> &        values) const;
 
    template <typename ForwardIterator, typename OutputIterator>
    void
    extract_subvector_to(ForwardIterator       indices_begin,
                         const ForwardIterator indices_end,
                         OutputIterator        values_begin) const;
 
    Number
    local_element(const size_type local_index) const;
 
    Number &
    local_element(const size_type local_index);
 
 
 
    template <typename Number2>
    void
    add(const std::vector<size_type> &indices,
        const std::vector<Number2> &  values);
 
    template <typename Number2>
    void
    add(const std::vector<size_type> &  indices,
        const ReadWriteVector<Number2> &values);
 
    template <typename Number2>
    void
    add(const size_type  n_elements,
        const size_type *indices,
        const Number2 *  values);
 
    void
    print(std::ostream &     out,
          const unsigned int precision  = 3,
          const bool         scientific = true) const;
 
    std::size_t
    memory_consumption() const;
 
  protected:
#ifdef DEAL_II_WITH_TRILINOS
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
    void
    import(
      const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
      const IndexSet &        locally_owned_elements,
      VectorOperation::values operation,
      const MPI_Comm &        mpi_comm,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern);
#  endif
 
    void
    import(const Epetra_MultiVector &multivector,
           const IndexSet &          locally_owned_elements,
           VectorOperation::values   operation,
           const MPI_Comm &          mpi_comm,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern);
#endif
 
    unsigned int
    global_to_local(const types::global_dof_index global_index) const
    {
      // the following will throw an exception if the global_index is not
      // in the remaining_elements
      return static_cast<unsigned int>(
        stored_elements.index_within_set(global_index));
    }
 
    void
    resize_val(const size_type new_allocated_size);
 
#if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
    TpetraWrappers::CommunicationPattern
    create_tpetra_comm_pattern(const IndexSet &source_index_set,
                               const MPI_Comm &mpi_comm);
#  endif
 
    EpetraWrappers::CommunicationPattern
    create_epetra_comm_pattern(const IndexSet &source_index_set,
                               const MPI_Comm &mpi_comm);
#endif
 
    IndexSet stored_elements;
 
    IndexSet source_stored_elements;
 
    std::shared_ptr<Utilities::MPI::CommunicationPatternBase> comm_pattern;
 
    std::unique_ptr<Number[], decltype(std::free) *> values;
 
    mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
      thread_loop_partitioner;
 
    // Make all other ReadWriteVector types friends.
    template <typename Number2>
    friend class ReadWriteVector;
 
  private:
    template <typename Functor>
    class FunctorTemplate
    {
    public:
      FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
 
      virtual void
      operator()(const size_type begin, const size_type end);
 
    private:
      ReadWriteVector &parent;
 
      const Functor &functor;
    };
  };
 
  /*---------------------------- Inline functions ---------------------------*/
 
#ifndef DOXYGEN
 
  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector()
    : Subscriptor()
    , 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
    ReadWriteVector<Number>::reinit(0, true);
  }
 
 
 
  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(
    const ReadWriteVector<Number> &v)
    : Subscriptor()
    , values(nullptr, free)
  {
    this->operator=(v);
  }
 
 
 
  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(const size_type size)
    : Subscriptor()
    , 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
    ReadWriteVector<Number>::reinit(size, false);
  }
 
 
 
  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(
    const IndexSet &locally_stored_indices)
    : Subscriptor()
    , 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
    ReadWriteVector<Number>::reinit(locally_stored_indices);
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::size_type
  ReadWriteVector<Number>::size() const
  {
    return stored_elements.size();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::size_type
  ReadWriteVector<Number>::n_elements() const
  {
    return stored_elements.n_elements();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::size_type
  ReadWriteVector<Number>::locally_owned_size() const
  {
    return stored_elements.n_elements();
  }
 
 
 
  template <typename Number>
  inline const IndexSet &
  ReadWriteVector<Number>::get_stored_elements() const
  {
    return stored_elements;
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::iterator
  ReadWriteVector<Number>::begin()
  {
    return values.get();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::begin() const
  {
    return values.get();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::iterator
  ReadWriteVector<Number>::end()
  {
    return values.get() + this->n_elements();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::end() const
  {
    return values.get() + this->n_elements();
  }
 
 
 
  template <typename Number>
  inline Number
  ReadWriteVector<Number>::operator()(const size_type global_index) const
  {
    return values[global_to_local(global_index)];
  }
 
 
 
  template <typename Number>
  inline Number &
  ReadWriteVector<Number>::operator()(const size_type global_index)
  {
    return values[global_to_local(global_index)];
  }
 
 
 
  template <typename Number>
  inline Number ReadWriteVector<Number>::
                operator[](const size_type global_index) const
  {
    return operator()(global_index);
  }
 
 
 
  template <typename Number>
  inline Number &ReadWriteVector<Number>::
                 operator[](const size_type global_index)
  {
    return operator()(global_index);
  }
 
 
 
  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::extract_subvector_to(
    const std::vector<size_type> &indices,
    std::vector<Number2> &        extracted_values) const
  {
    for (size_type i = 0; i < indices.size(); ++i)
      extracted_values[i] = operator()(indices[i]);
  }
 
 
 
  template <typename Number>
  template <typename ForwardIterator, typename OutputIterator>
  inline void
  ReadWriteVector<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 Number
  ReadWriteVector<Number>::local_element(const size_type local_index) const
  {
    AssertIndexRange(local_index, this->n_elements());
 
    return values[local_index];
  }
 
 
 
  template <typename Number>
  inline Number &
  ReadWriteVector<Number>::local_element(const size_type local_index)
  {
    AssertIndexRange(local_index, this->n_elements());
 
    return values[local_index];
  }
 
 
 
  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
                               const std::vector<Number2> &  values)
  {
    AssertDimension(indices.size(), values.size());
    add(indices.size(), indices.data(), values.data());
  }
 
 
 
  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const std::vector<size_type> &  indices,
                               const ReadWriteVector<Number2> &values)
  {
    const size_type size = indices.size();
    for (size_type i = 0; i < size; ++i)
      {
        Assert(
          numbers::is_finite(values[i]),
          ExcMessage(
            "The given value is not finite but either infinite or Not A Number (NaN)"));
        this->operator()(indices[i]) += values[indices[i]];
      }
  }
 
 
 
  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const size_type  n_indices,
                               const size_type *indices,
                               const Number2 *  values_to_add)
  {
    for (size_type i = 0; i < n_indices; ++i)
      {
        Assert(
          numbers::is_finite(values[i]),
          ExcMessage(
            "The given value is not finite but either infinite or Not A Number (NaN)"));
        this->operator()(indices[i]) += values_to_add[i];
      }
  }
 
 
 
  template <typename Number>
  template <typename Functor>
  inline ReadWriteVector<Number>::FunctorTemplate<Functor>::FunctorTemplate(
    ReadWriteVector<Number> &parent,
    const Functor &          functor)
    : parent(parent)
    , functor(functor)
  {}
 
 
 
  template <typename Number>
  template <typename Functor>
  void
  ReadWriteVector<Number>::FunctorTemplate<Functor>::
  operator()(const size_type begin, const size_type end)
  {
    for (size_type i = begin; i < end; ++i)
      functor(parent.values[i]);
  }
 
#endif // ifndef DOXYGEN
 
} // end of namespace LinearAlgebra
 
 
 
template <typename Number>
inline void
swap(LinearAlgebra::ReadWriteVector<Number> &u,
     LinearAlgebra::ReadWriteVector<Number> &v)
{
  u.swap(v);
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif