petsc_parallel_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2016 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_parallel_vector_h
#define dealii__petsc_parallel_vector_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/petsc_vector_base.h>
#  include <deal.II/base/index_set.h>

DEAL_II_NAMESPACE_OPEN



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


namespace PETScWrappers
{
  namespace MPI
  {

    class Vector : public VectorBase
    {
    public:
      typedef types::global_dof_index size_type;

      static const bool supports_distributed_data = true;

      Vector ();

      explicit Vector (const MPI_Comm  &communicator,
                       const size_type  n,
                       const size_type  local_size);


      template <typename Number>
      explicit Vector (const MPI_Comm               &communicator,
                       const ::Vector<Number> &v,
                       const size_type               local_size);


      explicit Vector (const MPI_Comm     &communicator,
                       const VectorBase   &v,
                       const size_type     local_size);

      Vector (const IndexSet &local,
              const IndexSet &ghost,
              const MPI_Comm &communicator);

      explicit Vector (const IndexSet &local,
                       const MPI_Comm &communicator);

      void clear ();

      Vector &operator= (const Vector &v);

      Vector &operator= (const PETScWrappers::Vector &v);

      Vector &operator= (const PetscScalar s);

      template <typename number>
      Vector &operator= (const ::Vector<number> &v);

      void reinit (const MPI_Comm  &communicator,
                   const size_type  N,
                   const size_type  local_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);

      const MPI_Comm &get_mpi_communicator () const;

      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 size_type n,
                                  const size_type local_size);



      virtual void create_vector (const size_type n,
                                  const size_type local_size,
                                  const IndexSet &ghostnodes);


    private:
      MPI_Comm communicator;
    };


// ------------------ 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         local_size)
      :
      communicator (communicator)
    {
      Vector::create_vector (v.size(), local_size);

      *this = v;
    }



    inline
    Vector &
    Vector::operator= (const PetscScalar s)
    {
      VectorBase::operator= (s);

      return *this;
    }



    inline
    Vector &
    Vector::operator= (const Vector &v)
    {
      // make sure left- and right-hand side of the assignment are compress()'ed:
      Assert(v.last_action == VectorOperation::unknown,
             internal::VectorReference::ExcWrongMode (VectorOperation::unknown,
                                                      v.last_action));
      Assert(last_action == VectorOperation::unknown,
             internal::VectorReference::ExcWrongMode (VectorOperation::unknown,
                                                      last_action));


      if (v.size()==0)
        {
          // this happens if v has not been initialized to something useful:
          // Vector x,v;x=v;
          // we skip the code below and create a simple serial vector of
          // length 0

          int ierr;
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
          ierr = VecDestroy (vector);
#else
          ierr = VecDestroy (&vector);
#endif
          AssertThrow (ierr == 0, ExcPETScError(ierr));

          const int n = 0;
          ierr = VecCreateSeq (PETSC_COMM_SELF, n, &vector);
          AssertThrow (ierr == 0, ExcPETScError(ierr));
          ghosted = false;
          ghost_indices.clear();
          return *this;
        }

      // if the vectors have different sizes,
      // then first resize the present one
      if (size() != v.size())
        {
          if (v.has_ghost_elements())
            reinit( v.locally_owned_elements(), v.ghost_indices, v.communicator);
          else
            reinit (v.communicator, v.size(), v.local_size(), true);
        }

      const int ierr = VecCopy (v.vector, vector);
      AssertThrow (ierr == 0, ExcPETScError(ierr));

      if (has_ghost_elements())
        {
          int 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));
        }
      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;
    }



    inline
    const MPI_Comm &
    Vector::get_mpi_communicator () const
    {
      return communicator;
    }

#endif // DOXYGEN
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC

/*----------------------------   petsc_parallel_vector.h     ---------------------------*/

#endif
/*----------------------------   petsc_parallel_vector.h     ---------------------------*/