petsc_vector.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2023 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_h
#define dealii_petsc_vector_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/partitioner.h>
#  include <deal.II/base/subscriptor.h>
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/petsc_vector_base.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_operation.h>
#  include <deal.II/lac/vector_type_traits.h>
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace PETScWrappers
{
  namespace MPI
  {
    class Vector : public VectorBase
    {
    public:
      using size_type = types::global_dof_index;
 
      Vector();
 
      using VectorBase::VectorBase;
 
      explicit Vector(const MPI_Comm  communicator,
                      const size_type n,
                      const size_type locally_owned_size);
 
      template <typename Number>
      explicit Vector(const MPI_Comm                communicator,
                      const ::Vector<Number> &v,
                      const size_type               locally_owned_size);
 
      Vector(const IndexSet &local,
             const IndexSet &ghost,
             const MPI_Comm  communicator);
 
      explicit Vector(const IndexSet &local, const MPI_Comm communicator);
 
      Vector(const Vector &v);
 
      virtual void
      clear() override;
 
      Vector &
      operator=(const Vector &v);
 
      Vector &
      operator=(const PetscScalar s);
 
      template <typename number>
      Vector &
      operator=(const ::Vector<number> &v);
 
      using VectorBase::reinit;
 
      void
      reinit(const MPI_Comm  communicator,
             const size_type N,
             const size_type locally_owned_size,
             const bool      omit_zeroing_entries = false);
 
      void
      reinit(const Vector &v, const bool omit_zeroing_entries = false);
 
      void
      reinit(const IndexSet &local,
             const IndexSet &ghost,
             const MPI_Comm  communicator);
 
      void
      reinit(const IndexSet &local, const MPI_Comm communicator);
 
      void
      reinit(
        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
        const bool make_ghosted = true);
 
      void
      print(std::ostream &     out,
            const unsigned int precision  = 3,
            const bool         scientific = true,
            const bool         across     = true) const;
 
      bool
      all_zero() const;
 
    protected:
      virtual void
      create_vector(const MPI_Comm  comm,
                    const size_type n,
                    const size_type locally_owned_size);
 
 
 
      virtual void
      create_vector(const MPI_Comm  comm,
                    const size_type n,
                    const size_type locally_owned_size,
                    const IndexSet &ghostnodes);
    };
 
 
    // ------------------ template and inline functions -------------
 
 
    inline void
    swap(Vector &u, Vector &v)
    {
      u.swap(v);
    }
 
 
#  ifndef DOXYGEN
 
    template <typename number>
    Vector::Vector(const MPI_Comm                communicator,
                   const ::Vector<number> &v,
                   const size_type               locally_owned_size)
    {
      Vector::create_vector(communicator, v.size(), locally_owned_size);
 
      *this = v;
    }
 
 
 
    inline Vector &
    Vector::operator=(const PetscScalar s)
    {
      VectorBase::operator=(s);
 
      return *this;
    }
 
 
 
    template <typename number>
    inline Vector &
    Vector::operator=(const ::Vector<number> &v)
    {
      Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
 
      // FIXME: the following isn't necessarily fast, but this is due to
      // the fact that PETSc doesn't offer an inlined access operator.
      //
      // if someone wants to contribute some code: to make this code
      // faster, one could either first convert all values to PetscScalar,
      // and then set them all at once using VecSetValues. This has the
      // drawback that it could take quite some memory, if the vector is
      // large, and it would in addition allocate memory on the heap, which
      // is expensive. an alternative would be to split the vector into
      // chunks of, say, 128 elements, convert a chunk at a time and set it
      // in the output vector using VecSetValues. since 128 elements is
      // small enough, this could easily be allocated on the stack (as a
      // local variable) which would make the whole thing much more
      // efficient.
      //
      // a second way to make things faster is for the special case that
      // number==PetscScalar. we could then declare a specialization of
      // this template, and omit the conversion. the problem with this is
      // that the best we can do is to use VecSetValues, but this isn't
      // very efficient either: it wants to see an array of indices, which
      // in this case a) again takes up a whole lot of memory on the heap,
      // and b) is totally dumb since its content would simply be the
      // sequence 0,1,2,3,...,n. the best of all worlds would probably be a
      // function in PETSc that would take a pointer to an array of
      // PetscScalar values and simply copy n elements verbatim into the
      // vector...
      for (size_type i = 0; i < v.size(); ++i)
        (*this)(i) = v(i);
 
      compress(VectorOperation::insert);
 
      return *this;
    }
 
 
 
#  endif // DOXYGEN
  }      // namespace MPI
} // namespace PETScWrappers
 
namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename>
    class ReinitHelper;
 
    template <>
    class ReinitHelper<PETScWrappers::MPI::Vector>
    {
    public:
      template <typename Matrix>
      static void
      reinit_range_vector(const Matrix &              matrix,
                          PETScWrappers::MPI::Vector &v,
                          bool /*omit_zeroing_entries*/)
      {
        v.reinit(matrix.locally_owned_range_indices(),
                 matrix.get_mpi_communicator());
      }
 
      template <typename Matrix>
      static void
      reinit_domain_vector(const Matrix &              matrix,
                           PETScWrappers::MPI::Vector &v,
                           bool /*omit_zeroing_entries*/)
      {
        v.reinit(matrix.locally_owned_domain_indices(),
                 matrix.get_mpi_communicator());
      }
    };
 
  } // namespace LinearOperatorImplementation
} /* namespace internal */
 
template <>
struct is_serial_vector<PETScWrappers::MPI::Vector> : std::false_type
{};
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_PETSC
 
#endif