petsc_vector_base.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2004 - 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_petsc_vector_base_h
#define dealii_petsc_vector_base_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/subscriptor.h>
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_operation.h>
 
#  include <boost/serialization/split_member.hpp>
#  include <boost/serialization/utility.hpp>
 
#  include <petscvec.h>
 
#  include <utility>
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declaration
#  ifndef DOXYGEN
template <typename number>
class Vector;
 
namespace PETScWrappers
{
  class VectorBase;
}
#  endif
 
namespace PETScWrappers
{
  namespace internal
  {
    class VectorReference
    {
    public:
      using size_type = types::global_dof_index;
 
    private:
      VectorReference(const VectorBase &vector, const size_type index);
 
    public:
      /*
       * Copy constructor.
       */
      VectorReference(const VectorReference &vector) = default;
 
      const VectorReference &
      operator=(const VectorReference &r) const;
 
      VectorReference &
      operator=(const VectorReference &r);
 
      const VectorReference &
      operator=(const PetscScalar &s) const;
 
      const VectorReference &
      operator+=(const PetscScalar &s) const;
 
      const VectorReference &
      operator-=(const PetscScalar &s) const;
 
      const VectorReference &
      operator*=(const PetscScalar &s) const;
 
      const VectorReference &
      operator/=(const PetscScalar &s) const;
 
      PetscReal
      real() const;
 
      PetscReal
      imag() const;
 
      operator PetscScalar() const;
      DeclException3(
        ExcAccessToNonlocalElement,
        int,
        int,
        int,
        << "You tried to access element " << arg1
        << " of a distributed vector, but only elements in range [" << arg2
        << ',' << arg3 << "] are stored locally and can be accessed."
        << "\n\n"
        << "A common source for this kind of problem is that you "
        << "are passing a 'fully distributed' vector into a function "
        << "that needs read access to vector elements that correspond "
        << "to degrees of freedom on ghost cells (or at least to "
        << "'locally active' degrees of freedom that are not also "
        << "'locally owned'). You need to pass a vector that has these "
        << "elements as ghost entries.");
      DeclException2(ExcWrongMode,
                     int,
                     int,
                     << "You tried to do a "
                     << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
                     << " operation but the vector is currently in "
                     << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
                     << " mode. You first have to call 'compress()'.");
 
    private:
      const VectorBase &vector;
 
      const size_type index;
 
      // Make the vector class a friend, so that it can create objects of the
      // present type.
      friend class ::PETScWrappers::VectorBase;
    };
  } // namespace internal
  class VectorBase : public ReadVector<PetscScalar>, public Subscriptor
  {
  public:
    using value_type      = PetscScalar;
    using real_type       = PetscReal;
    using size_type       = types::global_dof_index;
    using reference       = internal::VectorReference;
    using const_reference = const internal::VectorReference;
 
    VectorBase();
 
    VectorBase(const VectorBase &v);
 
    explicit VectorBase(const Vec &v);
 
    virtual ~VectorBase() override;
 
    virtual void
    clear();
 
    void
    compress(const VectorOperation::values operation);
 
    VectorBase &
    operator=(const VectorBase &);
 
    VectorBase &
    operator=(const PetscScalar s);
 
    void
    reinit(Vec v);
 
    bool
    operator==(const VectorBase &v) const;
 
    bool
    operator!=(const VectorBase &v) const;
 
    size_type
    size() const override;
 
    size_type
    locally_owned_size() const;
 
    std::pair<size_type, size_type>
    local_range() const;
 
    bool
    in_local_range(const size_type index) const;
 
    IndexSet
    locally_owned_elements() const;
 
    bool
    has_ghost_elements() const;
 
    const IndexSet &
    ghost_elements() const;
 
    void
    update_ghost_values() const;
 
    reference
    operator()(const size_type index);
 
    PetscScalar
    operator()(const size_type index) const;
 
    reference
    operator[](const size_type index);
 
    PetscScalar
    operator[](const size_type index) const;
 
    void
    set(const std::vector<size_type>   &indices,
        const std::vector<PetscScalar> &values);
 
    void
    extract_subvector_to(const std::vector<size_type> &indices,
                         std::vector<PetscScalar>     &values) const;
 
    virtual void
    extract_subvector_to(
      const ArrayView<const types::global_dof_index> &indices,
      ArrayView<PetscScalar>                         &elements) const override;
 
    template <typename ForwardIterator, typename OutputIterator>
    void
    extract_subvector_to(const ForwardIterator indices_begin,
                         const ForwardIterator indices_end,
                         OutputIterator        values_begin) const;
 
    void
    add(const std::vector<size_type>   &indices,
        const std::vector<PetscScalar> &values);
 
    void
    add(const std::vector<size_type>        &indices,
        const ::Vector<PetscScalar> &values);
 
    void
    add(const size_type    n_elements,
        const size_type   *indices,
        const PetscScalar *values);
 
    PetscScalar
    operator*(const VectorBase &vec) const;
 
    real_type
    norm_sqr() const;
 
    PetscScalar
    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;
 
    PetscScalar
    add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
 
    bool
    all_zero() const;
 
    VectorBase &
    operator*=(const PetscScalar factor);
 
    VectorBase &
    operator/=(const PetscScalar factor);
 
    VectorBase &
    operator+=(const VectorBase &V);
 
    VectorBase &
    operator-=(const VectorBase &V);
 
    void
    add(const PetscScalar s);
 
    void
    add(const PetscScalar a, const VectorBase &V);
 
    void
    add(const PetscScalar a,
        const VectorBase &V,
        const PetscScalar b,
        const VectorBase &W);
 
    void
    sadd(const PetscScalar s, const VectorBase &V);
 
    void
    sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
 
    void
    scale(const VectorBase &scaling_factors);
 
    void
    equ(const PetscScalar a, const VectorBase &V);
 
    void
    write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
 
    void
    print(std::ostream      &out,
          const unsigned int precision  = 3,
          const bool         scientific = true,
          const bool         across     = true) const;
 
    template <class Archive>
    void
    save(Archive &ar, const unsigned int version) const;
 
    template <class Archive>
    void
    load(Archive &ar, const unsigned int version);
 
#  ifdef DOXYGEN
    template <class Archive>
    void
    serialize(Archive &archive, const unsigned int version);
#  else
    // This macro defines the serialize() method that is compatible with
    // the templated save() and load() method that have been implemented.
    BOOST_SERIALIZATION_SPLIT_MEMBER()
#  endif
 
    void
    swap(VectorBase &v) noexcept;
 
    operator const Vec &() const;
 
    Vec &
    petsc_vector();
 
    std::size_t
    memory_consumption() const;
 
    MPI_Comm
    get_mpi_communicator() const;
 
  protected:
    Vec vector;
 
    bool ghosted;
 
    IndexSet ghost_indices;
 
    mutable VectorOperation::values last_action;
 
    // Make the reference class a friend.
    friend class internal::VectorReference;
 
    void
    do_set_add_operation(const size_type    n_elements,
                         const size_type   *indices,
                         const PetscScalar *values,
                         const bool         add_values);
 
    void
    determine_ghost_indices();
  };
 
 
 
  // ------------------- inline and template functions --------------
 
  inline void
  swap(VectorBase &u, VectorBase &v) noexcept
  {
    u.swap(v);
  }
 
#  ifndef DOXYGEN
  namespace internal
  {
    inline VectorReference::VectorReference(const VectorBase &vector,
                                            const size_type   index)
      : vector(vector)
      , index(index)
    {}
 
 
    inline const VectorReference &
    VectorReference::operator=(const VectorReference &r) const
    {
      // as explained in the class
      // documentation, this is not the copy
      // operator. so simply pass on to the
      // "correct" assignment operator
      *this = static_cast<PetscScalar>(r);
 
      return *this;
    }
 
 
 
    inline VectorReference &
    VectorReference::operator=(const VectorReference &r)
    {
      // as explained in the class
      // documentation, this is not the copy
      // operator. so simply pass on to the
      // "correct" assignment operator
      *this = static_cast<PetscScalar>(r);
 
      return *this;
    }
 
 
 
    inline const VectorReference &
    VectorReference::operator=(const PetscScalar &value) const
    {
      Assert((vector.last_action == VectorOperation::insert) ||
               (vector.last_action == VectorOperation::unknown),
             ExcWrongMode(VectorOperation::insert, vector.last_action));
 
      Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
      const PetscInt petsc_i = index;
 
      const PetscErrorCode ierr =
        VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
      AssertThrow(ierr == 0, ExcPETScError(ierr));
 
      vector.last_action = VectorOperation::insert;
 
      return *this;
    }
 
 
 
    inline const VectorReference &
    VectorReference::operator+=(const PetscScalar &value) const
    {
      Assert((vector.last_action == VectorOperation::add) ||
               (vector.last_action == VectorOperation::unknown),
             ExcWrongMode(VectorOperation::add, vector.last_action));
 
      Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
      vector.last_action = VectorOperation::add;
 
      // we have to do above actions in any
      // case to be consistent with the MPI
      // communication model (see the
      // comments in the documentation of
      // PETScWrappers::MPI::Vector), but we
      // can save some work if the addend is
      // zero
      if (value == PetscScalar())
        return *this;
 
      // use the PETSc function to add something
      const PetscInt       petsc_i = index;
      const PetscErrorCode ierr =
        VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
      AssertThrow(ierr == 0, ExcPETScError(ierr));
 
 
      return *this;
    }
 
 
 
    inline const VectorReference &
    VectorReference::operator-=(const PetscScalar &value) const
    {
      Assert((vector.last_action == VectorOperation::add) ||
               (vector.last_action == VectorOperation::unknown),
             ExcWrongMode(VectorOperation::add, vector.last_action));
 
      Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
      vector.last_action = VectorOperation::add;
 
      // we have to do above actions in any
      // case to be consistent with the MPI
      // communication model (see the
      // comments in the documentation of
      // PETScWrappers::MPI::Vector), but we
      // can save some work if the addend is
      // zero
      if (value == PetscScalar())
        return *this;
 
      // use the PETSc function to
      // add something
      const PetscInt       petsc_i     = index;
      const PetscScalar    subtractand = -value;
      const PetscErrorCode ierr =
        VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
      AssertThrow(ierr == 0, ExcPETScError(ierr));
 
      return *this;
    }
 
 
 
    inline const VectorReference &
    VectorReference::operator*=(const PetscScalar &value) const
    {
      Assert((vector.last_action == VectorOperation::insert) ||
               (vector.last_action == VectorOperation::unknown),
             ExcWrongMode(VectorOperation::insert, vector.last_action));
 
      Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
      vector.last_action = VectorOperation::insert;
 
      // we have to do above actions in any
      // case to be consistent with the MPI
      // communication model (see the
      // comments in the documentation of
      // PETScWrappers::MPI::Vector), but we
      // can save some work if the factor is
      // one
      if (value == 1.)
        return *this;
 
      const PetscInt    petsc_i   = index;
      const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
 
      const PetscErrorCode ierr =
        VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
      AssertThrow(ierr == 0, ExcPETScError(ierr));
 
      return *this;
    }
 
 
 
    inline const VectorReference &
    VectorReference::operator/=(const PetscScalar &value) const
    {
      Assert((vector.last_action == VectorOperation::insert) ||
               (vector.last_action == VectorOperation::unknown),
             ExcWrongMode(VectorOperation::insert, vector.last_action));
 
      Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
      vector.last_action = VectorOperation::insert;
 
      // we have to do above actions in any
      // case to be consistent with the MPI
      // communication model (see the
      // comments in the documentation of
      // PETScWrappers::MPI::Vector), but we
      // can save some work if the factor is
      // one
      if (value == 1.)
        return *this;
 
      const PetscInt    petsc_i   = index;
      const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
 
      const PetscErrorCode ierr =
        VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
      AssertThrow(ierr == 0, ExcPETScError(ierr));
 
      return *this;
    }
 
 
 
    inline PetscReal
    VectorReference::real() const
    {
#    ifndef PETSC_USE_COMPLEX
      return static_cast<PetscScalar>(*this);
#    else
      return PetscRealPart(static_cast<PetscScalar>(*this));
#    endif
    }
 
 
 
    inline PetscReal
    VectorReference::imag() const
    {
#    ifndef PETSC_USE_COMPLEX
      return PetscReal(0);
#    else
      return PetscImaginaryPart(static_cast<PetscScalar>(*this));
#    endif
    }
 
  } // namespace internal
 
  inline bool
  VectorBase::in_local_range(const size_type index) const
  {
    PetscInt             begin, end;
    const PetscErrorCode ierr =
      VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
 
    return ((index >= static_cast<size_type>(begin)) &&
            (index < static_cast<size_type>(end)));
  }
 
 
  inline IndexSet
  VectorBase::locally_owned_elements() const
  {
    IndexSet is(size());
 
    // PETSc only allows for contiguous local ranges, so this is simple
    const std::pair<size_type, size_type> x = local_range();
    is.add_range(x.first, x.second);
    return is;
  }
 
 
 
  inline bool
  VectorBase::has_ghost_elements() const
  {
    return ghosted;
  }
 
 
  inline const IndexSet &
  VectorBase::ghost_elements() const
  {
    return ghost_indices;
  }
 
 
  inline void
  VectorBase::update_ghost_values() const
  {
    if (ghosted)
      {
        PetscErrorCode ierr;
 
        ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
        ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }
  }
 
 
 
  inline internal::VectorReference
  VectorBase::operator()(const size_type index)
  {
    return internal::VectorReference(*this, index);
  }
 
 
 
  inline PetscScalar
  VectorBase::operator()(const size_type index) const
  {
    return static_cast<PetscScalar>(internal::VectorReference(*this, index));
  }
 
 
 
  inline internal::VectorReference
  VectorBase::operator[](const size_type index)
  {
    return operator()(index);
  }
 
 
 
  inline PetscScalar
  VectorBase::operator[](const size_type index) const
  {
    return operator()(index);
  }
 
  inline MPI_Comm
  VectorBase::get_mpi_communicator() const
  {
    return PetscObjectComm(reinterpret_cast<PetscObject>(vector));
  }
 
  inline void
  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
                                   std::vector<PetscScalar>     &values) const
  {
    Assert(indices.size() <= values.size(),
           ExcDimensionMismatch(indices.size(), values.size()));
    extract_subvector_to(indices.begin(), indices.end(), values.begin());
  }
 
  inline void
  VectorBase::extract_subvector_to(
    const ArrayView<const types::global_dof_index> &indices,
    ArrayView<PetscScalar>                         &elements) const
  {
    AssertDimension(indices.size(), elements.size());
    extract_subvector_to(indices.begin(), indices.end(), elements.begin());
  }
 
 
  template <typename ForwardIterator, typename OutputIterator>
  inline void
  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
                                   const ForwardIterator indices_end,
                                   OutputIterator        values_begin) const
  {
    if (indices_begin == indices_end)
      return;
 
    // if we are dealing
    // with a parallel vector
    if (ghosted)
      {
        // there is the possibility
        // that the vector has
        // ghost elements. in that
        // case, we first need to
        // figure out which
        // elements we own locally,
        // then get a pointer to
        // the elements that are
        // stored here (both the
        // ones we own as well as
        // the ghost elements). in
        // this array, the locally
        // owned elements come
        // first followed by the
        // ghost elements whose
        // position we can get from
        // an index set
        PetscInt       begin, end;
        PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        Vec locally_stored_elements = nullptr;
        ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        PetscInt lsize;
        ierr = VecGetSize(locally_stored_elements, &lsize);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        const PetscScalar *ptr;
        ierr = VecGetArrayRead(locally_stored_elements, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        auto input  = indices_begin;
        auto output = values_begin;
        while (input != indices_end)
          {
            const auto index = static_cast<PetscInt>(*input);
            AssertIntegerConversion(index, *input);
            if (index >= begin && index < end)
              {
                // local entry
                *output = *(ptr + index - begin);
              }
            else
              {
                // ghost entry
                const auto ghost_index = ghost_indices.index_within_set(*input);
 
                AssertIndexRange(ghost_index + end - begin, lsize);
                *output = *(ptr + ghost_index + end - begin);
              }
 
            ++input;
            ++output;
          }
 
        ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }
    // if the vector is local or the
    // caller, then simply access the
    // element we are interested in
    else
      {
        PetscInt       begin, end;
        PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        const PetscScalar *ptr;
        ierr = VecGetArrayRead(vector, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
        auto input  = indices_begin;
        auto output = values_begin;
        while (input != indices_end)
          {
            const auto index = static_cast<PetscInt>(*input);
            AssertIntegerConversion(index, *input);
 
            Assert(index >= begin && index < end,
                   ExcMessage("You are accessing elements of a vector without "
                              "ghost elements that are not actually owned by "
                              "this vector. A typical case where this may "
                              "happen is if you are passing a non-ghosted "
                              "(completely distributed) vector to a function "
                              "that expects a vector that stores ghost "
                              "elements for all locally relevant or locally "
                              "active vector entries."));
 
            *output = *(ptr + index - begin);
 
            ++input;
            ++output;
          }
 
        ierr = VecRestoreArrayRead(vector, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }
  }
 
  template <class Archive>
  inline void
  VectorBase::save(Archive &ar, const unsigned int) const
  {
    // forward to serialization function in the base class.
    ar &static_cast<const Subscriptor &>(*this);
    ar &size();
    ar &local_range();
 
    const PetscScalar *array = nullptr;
    int                ierr  = VecGetArrayRead(*this, &array);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
 
    boost::serialization::array_wrapper<const PetscScalar> wrapper(
      array, locally_owned_size());
    ar &wrapper;
 
    ierr = VecRestoreArrayRead(*this, &array);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }
 
 
 
  template <class Archive>
  inline void
  VectorBase::load(Archive &ar, const unsigned int)
  {
    ar &static_cast<Subscriptor &>(*this);
 
    size_type                       size = 0;
    std::pair<size_type, size_type> local_range;
 
    ar &size;
    Assert(size == this->size(),
           ExcMessage("The serialized value of size (" + std::to_string(size) +
                      ") does not match the current size (" +
                      std::to_string(this->size()) + ")"));
    (void)size;
    ar &local_range;
    Assert(local_range == this->local_range(),
           ExcMessage("The serialized value of local_range (" +
                      std::to_string(local_range.first) + ", " +
                      std::to_string(local_range.second) +
                      ") does not match the current local_range (" +
                      std::to_string(this->local_range().first) + ", " +
                      std::to_string(this->local_range().second) + ")"));
    (void)local_range;
 
    PetscScalar *array = nullptr;
    int          ierr  = VecGetArray(petsc_vector(), &array);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
 
    boost::serialization::array_wrapper<PetscScalar> wrapper(
      array, locally_owned_size());
    ar &wrapper;
 
    ierr = VecRestoreArray(petsc_vector(), &array);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }
#  endif // DOXYGEN
} // namespace PETScWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_PETSC
 
#endif