petsc_vector_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 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_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 <petscvec.h>

#    include <utility>
#    include <vector>

DEAL_II_NAMESPACE_OPEN

// forward declaration
template <typename number>
class Vector;


namespace PETScWrappers
{
  // forward declaration
  class VectorBase;

  namespace internal
  {
    class VectorReference
    {
    public:
      using size_type = types::global_dof_index;

    private:
      VectorReference(const VectorBase &vector, const size_type index);


    public:
      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 " << arg2
                     << " through " << arg3
                     << " are stored locally and can be accessed.");
      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;

      friend class ::PETScWrappers::VectorBase;
    };
  } // namespace internal
  class VectorBase : 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 PetscScalar s);

    bool
    operator==(const VectorBase &v) const;

    bool
    operator!=(const VectorBase &v) const;

    size_type
    size() const;

    size_type
    local_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;

    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;

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

    real_type
    min() const;

    real_type
    max() const;

    bool
    all_zero() const;

    bool
    is_non_negative() 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);

    DEAL_II_DEPRECATED
    void
    ratio(const VectorBase &a, const VectorBase &b);

    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;

    void
    swap(VectorBase &v);

    operator const Vec &() const;

    std::size_t
    memory_consumption() const;

    virtual const MPI_Comm &
    get_mpi_communicator() const;

  protected:
    Vec vector;

    bool ghosted;

    IndexSet ghost_indices;

    mutable VectorOperation::values last_action;

    friend class internal::VectorReference;

    bool obtained_ownership;

    void
    do_set_add_operation(const size_type    n_elements,
                         const size_type *  indices,
                         const PetscScalar *values,
                         const bool         add_values);

  private:
    VectorBase &
    operator=(const VectorBase &) = delete;
  };



  // ------------------- inline and template functions --------------

  inline void
  swap(VectorBase &u, VectorBase &v)
  {
    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 void
  VectorBase::update_ghost_values() const
  {}



  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 const MPI_Comm &
  VectorBase::get_mpi_communicator() const
  {
    static MPI_Comm comm;
    PetscObjectGetComm((PetscObject)vector, &comm);
    return comm;
  }

  inline void
  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
                                   std::vector<PetscScalar> &    values) const
  {
    extract_subvector_to(&(indices[0]),
                         &(indices[0]) + indices.size(),
                         &(values[0]));
  }

  template <typename ForwardIterator, typename OutputIterator>
  inline void
  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
                                   const ForwardIterator indices_end,
                                   OutputIterator        values_begin) const
  {
    const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
    if (n_idx == 0)
      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));

        PetscScalar *ptr;
        ierr = VecGetArray(locally_stored_elements, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        for (PetscInt i = 0; i < n_idx; ++i)
          {
            const unsigned int index = *(indices_begin + i);
            if (index >= static_cast<unsigned int>(begin) &&
                index < static_cast<unsigned int>(end))
              {
                // local entry
                *(values_begin + i) = *(ptr + index - begin);
              }
            else
              {
                // ghost entry
                const unsigned int ghostidx =
                  ghost_indices.index_within_set(index);

                Assert(ghostidx + end - begin < (unsigned int)lsize,
                       ExcInternalError());
                *(values_begin + i) = *(ptr + ghostidx + end - begin);
              }
          }

        ierr = VecRestoreArray(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));

        PetscScalar *ptr;
        ierr = VecGetArray(vector, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        for (PetscInt i = 0; i < n_idx; ++i)
          {
            const unsigned int index = *(indices_begin + i);

            Assert(index >= static_cast<unsigned int>(begin) &&
                     index < static_cast<unsigned int>(end),
                   ExcInternalError());

            *(values_begin + i) = *(ptr + index - begin);
          }

        ierr = VecRestoreArray(vector, &ptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }
  }

#    endif // DOXYGEN
} // namespace PETScWrappers

DEAL_II_NAMESPACE_CLOSE

#  endif // DEAL_II_WITH_PETSC

#endif
/*---------------------------- petsc_vector_base.h --------------------------*/