read_write_vector.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2015 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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/aligned_vector.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_stub.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/lac/read_vector.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
template <typename>
class Vector;
 
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 ReadVector<Number>
  {
  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
    void
    reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
#endif
 
    template <typename Functor>
    void
    apply(const Functor &func);
 
    void
    swap(ReadWriteVector<Number> &v) noexcept;
 
    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_elements(
      const ::Vector<Number> &vec,
      VectorOperation::values       operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    DEAL_II_DEPRECATED
    void
    import(const ::Vector<Number> &V,
           VectorOperation::values       operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
 
    template <typename MemorySpace>
    void
    import_elements(
      const distributed::Vector<Number, MemorySpace> &vec,
      VectorOperation::values                         operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    template <typename MemorySpace>
    DEAL_II_DEPRECATED void
    import(const distributed::Vector<Number, MemorySpace> &V,
           VectorOperation::values                         operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
 
 
#ifdef DEAL_II_WITH_PETSC
    void
    import_elements(
      const PETScWrappers::MPI::Vector &petsc_vec,
      VectorOperation::values           operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    DEAL_II_DEPRECATED
    void
    import(const PETScWrappers::MPI::Vector &V,
           VectorOperation::values           operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
#endif
 
#ifdef DEAL_II_WITH_TRILINOS
    void
    import_elements(
      const TrilinosWrappers::MPI::Vector &trilinos_vec,
      VectorOperation::values              operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    DEAL_II_DEPRECATED
    void
    import(const TrilinosWrappers::MPI::Vector &V,
           VectorOperation::values              operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
 
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
    template <typename MemorySpace, typename Dummy = Number>
    std::enable_if_t<std::is_same_v<Dummy, Number> &&
                     ::is_tpetra_type<Number>::value>
    import_elements(
      const TpetraWrappers::Vector<Number, MemorySpace> &tpetra_vec,
      VectorOperation::values                            operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    template <typename MemorySpace, typename Dummy = Number>
    DEAL_II_DEPRECATED std::enable_if_t<std::is_same_v<Dummy, Number> &&
                                        ::is_tpetra_type<Number>::value>
    import(const TpetraWrappers::Vector<Number, MemorySpace> &V,
           VectorOperation::values                            operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
#  endif
 
    void
    import_elements(
      const EpetraWrappers::Vector &epetra_vec,
      VectorOperation::values       operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
    DEAL_II_DEPRECATED
    void
    import(const EpetraWrappers::Vector &V,
           VectorOperation::values       operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
#endif
 
#ifdef DEAL_II_WITH_CUDA
    void
    import_elements(
      const CUDAWrappers::Vector<Number> &cuda_vec,
      VectorOperation::values             operation,
      const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
        &communication_pattern = {});
 
 
    DEAL_II_DEPRECATED
    void
    import(const CUDAWrappers::Vector<Number> &V,
           VectorOperation::values             operation,
           const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
             &communication_pattern = {})
    {
      import_elements(V, operation, communication_pattern);
    }
#endif
 
    size_type
    size() const override;
 
    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;
 
    virtual void
    extract_subvector_to(
      const ArrayView<const types::global_dof_index> &indices,
      ArrayView<Number>                              &entries) const override;
 
    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
    template <typename MemorySpace, typename Dummy = Number>
    std::enable_if_t<std::is_same_v<Dummy, Number> &&
                     ::is_tpetra_type<Number>::value>
    import_elements(
      const Tpetra::
        Vector<Number, int, types::signed_global_dof_index, MemorySpace>
                             &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_elements(
      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;
 
#ifdef DEAL_II_WITH_TRILINOS
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
    template <typename MemorySpace = ::MemorySpace::Host>
    TpetraWrappers::CommunicationPattern<MemorySpace>
    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;
 
    AlignedVector<Number> 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()
  {
    // 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()
  {
    this->operator=(v);
  }
 
 
 
  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(const size_type size)
  {
    // 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)
  {
    // 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>::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.begin();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::begin() const
  {
    return values.begin();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::iterator
  ReadWriteVector<Number>::end()
  {
    return values.end();
  }
 
 
 
  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::end() const
  {
    return values.end();
  }
 
 
 
  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>
  void
  ReadWriteVector<Number>::extract_subvector_to(
    const ArrayView<const types::global_dof_index> &indices,
    ArrayView<Number>                              &entries) const
  {
    AssertDimension(indices.size(), entries.size());
    for (unsigned int i = 0; i < indices.size(); ++i)
      entries[i] = (*this)[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->locally_owned_size());
 
    return values[local_index];
  }
 
 
 
  template <typename Number>
  inline Number &
  ReadWriteVector<Number>::local_element(const size_type local_index)
  {
    AssertIndexRange(local_index, this->locally_owned_size());
 
    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>
  inline unsigned int
  ReadWriteVector<Number>::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));
  }
 
 
 
  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) noexcept
{
  u.swap(v);
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif