petsc_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2017 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_h
#define dealii__petsc_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/petsc_vector_base.h>
#  include <deal.II/lac/petsc_parallel_vector.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_type_traits.h>

DEAL_II_NAMESPACE_OPEN


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

    static const bool supports_distributed_data DEAL_II_DEPRECATED = false;

    Vector ();

    explicit Vector (const size_type n);

    template <typename Number>
    explicit Vector (const ::Vector<Number> &v);

    explicit Vector (const Vec &v);

    Vector (const Vector &v);

    explicit Vector (const MPI::Vector &v);

    void clear ();

    Vector &operator = (const Vector &v);

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

    Vector &operator = (const PetscScalar s);

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

    void reinit (const size_type N,
                 const bool      omit_zeroing_entries = false);

    void reinit (const Vector &v,
                 const bool    omit_zeroing_entries = false);

  protected:
    void create_vector (const size_type n);
  } DEAL_II_DEPRECATED;

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


  inline
  void swap (Vector &u, Vector &v)
  {
    u.swap (v);
  }


#ifndef DOXYGEN

  template <typename number>
  Vector::Vector (const ::Vector<number> &v)
  {
    Vector::create_vector (v.size());

    *this = v;
  }



  inline
  Vector::Vector (const Vec &v)
    :
    VectorBase(v)
  {}


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

    return *this;
  }


  inline
  Vector &
  Vector::operator = (const Vector &v)
  {
    // if the vectors have different sizes,
    // then first resize the present one
    if (size() != v.size())
      reinit (v.size(), true);

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

    return *this;
  }



  inline
  Vector &
  Vector::operator = (const MPI::Vector &v)
  {
    PetscErrorCode ierr;
    if (attained_ownership)
      {
        // the petsc function we call wants to
        // generate the vector itself, so destroy
        // the old one first
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
        ierr = VecDestroy (vector);
#else
        ierr = VecDestroy (&vector);
#endif
        AssertThrow (ierr == 0, ExcPETScError(ierr));
      }

    attained_ownership = true;

    // then do the gather
    // operation. <rant>petsc has changed its
    // interface again, and replaced a single
    // function call by several calls that
    // are hard to understand. gets me all
    // annoyed at their development
    // model</rant>
#if DEAL_II_PETSC_VERSION_LT(2,2,0)
    ierr = VecConvertMPIToSeqAll (static_cast<const Vec &>(v),
                                  &vector);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#else

    VecScatter ctx;

    ierr = VecScatterCreateToAll (static_cast<const Vec &>(v), &ctx, &vector);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#if ((PETSC_VERSION_MAJOR == 2) && \
     ((PETSC_VERSION_MINOR < 3) || \
      ((PETSC_VERSION_MINOR == 3) &&            \
       (PETSC_VERSION_SUBMINOR < 3))))
    ierr = VecScatterBegin (static_cast<const Vec &>(v), vector,
                            INSERT_VALUES, SCATTER_FORWARD, ctx);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = VecScatterEnd (static_cast<const Vec &>(v), vector,
                          INSERT_VALUES, SCATTER_FORWARD, ctx);

#else

    ierr = VecScatterBegin (ctx,static_cast<const Vec &>(v), vector,
                            INSERT_VALUES, SCATTER_FORWARD);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = VecScatterEnd (ctx, static_cast<const Vec &>(v), vector,
                          INSERT_VALUES, SCATTER_FORWARD);

#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    ierr = VecScatterDestroy (ctx);
#else
    ierr = VecScatterDestroy (&ctx);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));
#endif

    return *this;
  }



  template <typename number>
  inline
  Vector &
  Vector::operator = (const ::Vector<number> &v)
  {
    reinit (v.size(), true);
    // 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 internal
{
  namespace LinearOperator
  {
    template <typename> class ReinitHelper;

    template<>
    class ReinitHelper<PETScWrappers::Vector>
    {
    public:
      template <typename Matrix>
      static
      void reinit_range_vector (const Matrix &matrix,
                                PETScWrappers::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit( matrix.m(),
                  omit_zeroing_entries);
      }

      template <typename Matrix>
      static
      void reinit_domain_vector(const Matrix &matrix,
                                PETScWrappers::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit( matrix.n(),
                  omit_zeroing_entries);
      }
    };

  } /* namespace LinearOperator */
} /* namespace internal */


template <>
struct is_serial_vector< PETScWrappers::Vector > : std_cxx11::true_type
{
};


DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC

/*----------------------------   petsc_vector.h     ---------------------------*/

#endif
/*----------------------------   petsc_vector.h     ---------------------------*/