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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/subscriptor.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_operation.h>

#  include <vector>
#  include <utility>

#  include <petscvec.h>
#  include <deal.II/base/index_set.h>

DEAL_II_NAMESPACE_OPEN

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


namespace PETScWrappers
{
  // forward declaration
  class VectorBase;

  namespace internal
  {
    class VectorReference
    {
    public:
      typedef types::global_dof_index size_type;

    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;
    };
  }
  class VectorBase : public Subscriptor
  {
  public:
    typedef PetscScalar               value_type;
    typedef PetscReal                 real_type;
    typedef types::global_dof_index   size_type;
    typedef internal::VectorReference reference;
    typedef const internal::VectorReference const_reference;

    VectorBase ();

    VectorBase (const VectorBase &v);

    explicit VectorBase (const Vec &v);

    virtual ~VectorBase ();

    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
}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC

/*----------------------------   petsc_vector_base.h     ---------------------------*/

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