trilinos_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 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_trilinos_vector_h
#define dealii_trilinos_vector_h


#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/base/subscriptor.h>
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/utilities.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_operation.h>
#  include <deal.II/lac/vector_type_traits.h>

#  include <vector>
#  include <utility>
#  include <memory>

#  include <Epetra_ConfigDefs.h>
#  ifdef DEAL_II_WITH_MPI // only if MPI is installed
#    include <mpi.h>
#    include <Epetra_MpiComm.h>
#  else
#    include <Epetra_SerialComm.h>
#  endif
#  include <Epetra_FEVector.h>
#  include <Epetra_Map.h>
#  include <Epetra_LocalMap.h>

DEAL_II_NAMESPACE_OPEN


namespace TrilinosWrappers
{
  class SparseMatrix;

  namespace internal
  {
    typedef ::types::global_dof_index size_type;

    class VectorReference
    {
    private:
      VectorReference (MPI::Vector     &vector,
                       const size_type  index);

    public:

      const VectorReference &
      operator= (const VectorReference &r) const;

      VectorReference &
      operator= (const VectorReference &r);

      const VectorReference &
      operator= (const TrilinosScalar &s) const;

      const VectorReference &
      operator+= (const TrilinosScalar &s) const;

      const VectorReference &
      operator-= (const TrilinosScalar &s) const;

      const VectorReference &
      operator*= (const TrilinosScalar &s) const;

      const VectorReference &
      operator/= (const TrilinosScalar &s) const;

      operator TrilinosScalar () const;

      DeclException1 (ExcTrilinosError,
                      int,
                      << "An error with error number " << arg1
                      << " occurred while calling a Trilinos function");

    private:
      MPI::Vector   &vector;

      const size_type  index;

      friend class ::TrilinosWrappers::MPI::Vector;
    };
  }
  namespace
  {
#ifndef DEAL_II_WITH_64BIT_INDICES
    // define a helper function that queries the global ID of local ID of
    // an Epetra_BlockMap object  by calling either the 32- or 64-bit
    // function necessary.
    inline
    int gid(const Epetra_BlockMap &map, int i)
    {
      return map.GID(i);
    }
#else
    // define a helper function that queries the global ID of local ID of
    // an Epetra_BlockMap object  by calling either the 32- or 64-bit
    // function necessary.
    inline
    long long int gid(const Epetra_BlockMap &map, int i)
    {
      return map.GID64(i);
    }
#endif
  }

  namespace MPI
  {
    class BlockVector;

    class Vector : public Subscriptor
    {
    public:
      typedef TrilinosScalar                  value_type;
      typedef TrilinosScalar                  real_type;
      typedef ::types::global_dof_index size_type;
      typedef value_type                     *iterator;
      typedef const value_type               *const_iterator;
      typedef internal::VectorReference       reference;
      typedef const internal::VectorReference const_reference;

      Vector ();

      Vector (const Vector &v);

      explicit Vector (const IndexSet &parallel_partitioning,
                       const MPI_Comm &communicator = MPI_COMM_WORLD);

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

      Vector (const IndexSet   &parallel_partitioning,
              const Vector &v,
              const MPI_Comm   &communicator = MPI_COMM_WORLD);

      template <typename Number>
      Vector (const IndexSet               &parallel_partitioning,
              const ::Vector<Number> &v,
              const MPI_Comm               &communicator = MPI_COMM_WORLD);

      Vector (Vector &&v) noexcept;

      ~Vector () = default;

      void clear ();

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

      void reinit (const IndexSet &parallel_partitioning,
                   const MPI_Comm &communicator = MPI_COMM_WORLD,
                   const bool      omit_zeroing_entries = false);

      void reinit (const IndexSet &locally_owned_entries,
                   const IndexSet &ghost_entries,
                   const MPI_Comm &communicator = MPI_COMM_WORLD,
                   const bool      vector_writable = false);

      void reinit (const BlockVector &v,
                   const bool         import_data = false);

      void compress (::VectorOperation::values operation);

      Vector &operator= (const TrilinosScalar s);

      Vector &operator= (const Vector &v);

      Vector &operator= (Vector &&v) noexcept;

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

      void import_nonlocal_data_for_fe
      (const ::TrilinosWrappers::SparseMatrix &matrix,
       const Vector                                 &vector);

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

      bool operator!= (const Vector &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;

      TrilinosScalar operator* (const Vector &vec) const;

      real_type norm_sqr () const;

      TrilinosScalar mean_value () const;

      TrilinosScalar min () const;

      TrilinosScalar max () const;

      real_type l1_norm () const;

      real_type l2_norm () const;

      real_type lp_norm (const TrilinosScalar p) const;

      real_type linfty_norm () const;

      TrilinosScalar add_and_dot (const TrilinosScalar a,
                                  const Vector    &V,
                                  const Vector    &W);

      bool all_zero () const;

      bool is_non_negative () const;



      reference
      operator() (const size_type index);

      TrilinosScalar
      operator() (const size_type index) const;

      reference
      operator[] (const size_type index);

      TrilinosScalar
      operator[] (const size_type index) const;

      void extract_subvector_to (const std::vector<size_type> &indices,
                                 std::vector<TrilinosScalar> &values) const;

      template <typename ForwardIterator, typename OutputIterator>
      void extract_subvector_to (ForwardIterator          indices_begin,
                                 const ForwardIterator    indices_end,
                                 OutputIterator           values_begin) const;

      iterator begin ();

      const_iterator begin () const;

      iterator end ();

      const_iterator end () const;




      void set (const std::vector<size_type>    &indices,
                const std::vector<TrilinosScalar>  &values);

      void set (const std::vector<size_type>        &indices,
                const ::Vector<TrilinosScalar> &values);

      void set (const size_type       n_elements,
                const size_type      *indices,
                const TrilinosScalar *values);

      void add (const std::vector<size_type>      &indices,
                const std::vector<TrilinosScalar> &values);

      void add (const std::vector<size_type>           &indices,
                const ::Vector<TrilinosScalar> &values);

      void add (const size_type       n_elements,
                const size_type      *indices,
                const TrilinosScalar *values);

      Vector &operator*= (const TrilinosScalar factor);

      Vector &operator/= (const TrilinosScalar factor);

      Vector &operator+= (const Vector &V);

      Vector &operator-= (const Vector &V);

      void add (const TrilinosScalar s);

      void add (const Vector &V,
                const bool        allow_different_maps = false);

      void add (const TrilinosScalar  a,
                const Vector     &V);

      void add (const TrilinosScalar  a,
                const Vector     &V,
                const TrilinosScalar  b,
                const Vector     &W);

      void sadd (const TrilinosScalar  s,
                 const Vector     &V);

      void sadd (const TrilinosScalar  s,
                 const TrilinosScalar  a,
                 const Vector     &V);

      void scale (const Vector &scaling_factors);

      void equ (const TrilinosScalar  a,
                const Vector     &V);


      const Epetra_MultiVector &trilinos_vector () const;

      Epetra_FEVector &trilinos_vector ();

      const Epetra_Map &vector_partitioner () const;

      void print (std::ostream       &out,
                  const unsigned int  precision  = 3,
                  const bool          scientific = true,
                  const bool          across     = true) const;

      void swap (Vector &v);

      std::size_t memory_consumption () const;

      const MPI_Comm &get_mpi_communicator () const;

      DeclException0 (ExcDifferentParallelPartitioning);

      DeclException1 (ExcTrilinosError,
                      int,
                      << "An error with error number " << arg1
                      << " occurred while calling a Trilinos function");

      DeclException4 (ExcAccessToNonLocalElement,
                      size_type, size_type, size_type, size_type,
                      << "You tried to access element " << arg1
                      << " of a distributed vector, but this element is not stored "
                      << "on the current processor. Note: There are "
                      << arg2 << " elements stored "
                      << "on the current processor from within the range "
                      << arg3 << " through " << arg4
                      << " but Trilinos vectors need not store contiguous "
                      << "ranges on each processor, and not every element in "
                      << "this range may in fact be stored locally.");

    private:
      Epetra_CombineMode last_action;

      bool compressed;

      bool has_ghosts;

      std::unique_ptr<Epetra_FEVector> vector;

      std::unique_ptr<Epetra_MultiVector> nonlocal_vector;

      IndexSet owned_elements;

      friend class internal::VectorReference;
    };




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


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

#ifndef DOXYGEN

  namespace internal
  {
    inline
    VectorReference::VectorReference (MPI::Vector      &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<TrilinosScalar> (r);

      return *this;
    }



    inline
    VectorReference &
    VectorReference::operator= (const VectorReference &r)
    {
      // as above
      *this = static_cast<TrilinosScalar> (r);

      return *this;
    }


    inline
    const VectorReference &
    VectorReference::operator= (const TrilinosScalar &value) const
    {
      vector.set (1, &index, &value);
      return *this;
    }



    inline
    const VectorReference &
    VectorReference::operator+= (const TrilinosScalar &value) const
    {
      vector.add (1, &index, &value);
      return *this;
    }



    inline
    const VectorReference &
    VectorReference::operator-= (const TrilinosScalar &value) const
    {
      TrilinosScalar new_value = -value;
      vector.add (1, &index, &new_value);
      return *this;
    }



    inline
    const VectorReference &
    VectorReference::operator*= (const TrilinosScalar &value) const
    {
      TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
      vector.set (1, &index, &new_value);
      return *this;
    }



    inline
    const VectorReference &
    VectorReference::operator/= (const TrilinosScalar &value) const
    {
      TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
      vector.set (1, &index, &new_value);
      return *this;
    }
  }

  namespace MPI
  {
    inline
    bool
    Vector::in_local_range (const size_type index) const
    {
      std::pair<size_type, size_type> range = local_range();

      return ((index >= range.first) && (index <  range.second));
    }



    inline
    IndexSet
    Vector::locally_owned_elements() const
    {
      Assert(owned_elements.size()==size(),
             ExcMessage("The locally owned elements have not been properly initialized!"
                        " This happens for example if this object has been initialized"
                        " with exactly one overlapping IndexSet."));
      return owned_elements;
    }



    inline
    bool
    Vector::has_ghost_elements() const
    {
      return has_ghosts;
    }



    inline
    void
    Vector::update_ghost_values () const
    {}



    inline
    internal::VectorReference
    Vector::operator() (const size_type index)
    {
      return internal::VectorReference (*this, index);
    }



    inline
    internal::VectorReference
    Vector::operator[] (const size_type index)
    {
      return operator() (index);
    }



    inline
    TrilinosScalar
    Vector::operator[] (const size_type index) const
    {
      return operator() (index);
    }



    inline
    void Vector::extract_subvector_to (const std::vector<size_type> &indices,
                                       std::vector<TrilinosScalar>  &values) const
    {
      for (size_type i = 0; i < indices.size(); ++i)
        values[i] = operator()(indices[i]);
    }



    template <typename ForwardIterator, typename OutputIterator>
    inline
    void Vector::extract_subvector_to (ForwardIterator          indices_begin,
                                       const ForwardIterator    indices_end,
                                       OutputIterator           values_begin) const
    {
      while (indices_begin != indices_end)
        {
          *values_begin = operator()(*indices_begin);
          indices_begin++;
          values_begin++;
        }
    }



    inline
    Vector::iterator
    Vector::begin()
    {
      return (*vector)[0];
    }



    inline
    Vector::iterator
    Vector::end()
    {
      return (*vector)[0]+local_size();
    }



    inline
    Vector::const_iterator
    Vector::begin() const
    {
      return (*vector)[0];
    }



    inline
    Vector::const_iterator
    Vector::end() const
    {
      return (*vector)[0]+local_size();
    }



    inline
    void
    Vector::set (const std::vector<size_type>      &indices,
                 const std::vector<TrilinosScalar>  &values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      Assert (indices.size() == values.size(),
              ExcDimensionMismatch(indices.size(),values.size()));

      set (indices.size(), indices.data(), values.data());
    }



    inline
    void
    Vector::set (const std::vector<size_type>           &indices,
                 const ::Vector<TrilinosScalar> &values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      Assert (indices.size() == values.size(),
              ExcDimensionMismatch(indices.size(),values.size()));

      set (indices.size(), indices.data(), values.begin());
    }



    inline
    void
    Vector::set (const size_type       n_elements,
                 const size_type      *indices,
                 const TrilinosScalar *values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      if (last_action == Add)
        {
          const int ierr = vector->GlobalAssemble(Add);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        }

      if (last_action != Insert)
        last_action = Insert;

      for (size_type i=0; i<n_elements; ++i)
        {
          const size_type row = indices[i];
          const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
          if (local_row != -1)
            (*vector)[0][local_row] = values[i];
          else
            {
              const int ierr = vector->ReplaceGlobalValues (1,
                                                            (const TrilinosWrappers::types::int_type *)(&row),
                                                            &values[i]);
              AssertThrow (ierr == 0, ExcTrilinosError(ierr));
              compressed = false;
            }
          // in set operation, do not use the pre-allocated vector for nonlocal
          // entries even if it exists. This is to ensure that we really only
          // set the elements touched by the set() method and not all contained
          // in the nonlocal entries vector (there is no way to distinguish them
          // on the receiving processor)
        }
    }



    inline
    void
    Vector::add (const std::vector<size_type>      &indices,
                 const std::vector<TrilinosScalar>  &values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (indices.size() == values.size(),
              ExcDimensionMismatch(indices.size(),values.size()));

      add (indices.size(), indices.data(), values.data());
    }



    inline
    void
    Vector::add (const std::vector<size_type>           &indices,
                 const ::Vector<TrilinosScalar> &values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (indices.size() == values.size(),
              ExcDimensionMismatch(indices.size(),values.size()));

      add (indices.size(), indices.data(), values.begin());
    }



    inline
    void
    Vector::add (const size_type       n_elements,
                 const size_type      *indices,
                 const TrilinosScalar *values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      if (last_action != Add)
        {
          if (last_action == Insert)
            {
              const int ierr = vector->GlobalAssemble(Insert);
              AssertThrow (ierr == 0, ExcTrilinosError(ierr));
            }
          last_action = Add;
        }

      for (size_type i=0; i<n_elements; ++i)
        {
          const size_type row = indices[i];
          const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
          if (local_row != -1)
            (*vector)[0][local_row] += values[i];
          else if (nonlocal_vector.get() == nullptr)
            {
              const int ierr = vector->SumIntoGlobalValues (1,
                                                            (const TrilinosWrappers::types::int_type *)(&row),
                                                            &values[i]);
              AssertThrow (ierr == 0, ExcTrilinosError(ierr));
              compressed = false;
            }
          else
            {
              // use pre-allocated vector for non-local entries if it exists for
              // addition operation
              const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
              Assert(my_row != -1,
                     ExcMessage("Attempted to write into off-processor vector entry "
                                "that has not be specified as being writable upon "
                                "initialization"));
              (*nonlocal_vector)[0][my_row] += values[i];
              compressed = false;
            }
        }
    }



    inline
    Vector::size_type
    Vector::size () const
    {
#ifndef DEAL_II_WITH_64BIT_INDICES
      return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
#else
      return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
#endif
    }



    inline
    Vector::size_type
    Vector::local_size () const
    {
      return (size_type) vector->Map().NumMyElements();
    }



    inline
    std::pair<Vector::size_type, Vector::size_type>
    Vector::local_range () const
    {
#ifndef DEAL_II_WITH_64BIT_INDICES
      const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
      const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
#else
      const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
      const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
#endif

      Assert (end-begin == vector->Map().NumMyElements(),
              ExcMessage ("This function only makes sense if the elements that this "
                          "vector stores on the current processor form a contiguous range. "
                          "This does not appear to be the case for the current vector."));

      return std::make_pair (begin, end);
    }



    inline
    TrilinosScalar
    Vector::operator* (const Vector &vec) const
    {
      Assert (vector->Map().SameAs(vec.vector->Map()),
              ExcDifferentParallelPartitioning());
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      TrilinosScalar result;

      const int ierr = vector->Dot(*(vec.vector), &result);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return result;
    }



    inline
    Vector::real_type
    Vector::norm_sqr () const
    {
      const TrilinosScalar d = l2_norm();
      return d*d;
    }



    inline
    TrilinosScalar
    Vector::mean_value () const
    {
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      TrilinosScalar mean;
      const int ierr = vector->MeanValue (&mean);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return mean;
    }



    inline
    TrilinosScalar
    Vector::min () const
    {
      TrilinosScalar min_value;
      const int ierr = vector->MinValue (&min_value);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return min_value;
    }



    inline
    TrilinosScalar
    Vector::max () const
    {
      TrilinosScalar max_value;
      const int ierr = vector->MaxValue (&max_value);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return max_value;
    }



    inline
    Vector::real_type
    Vector::l1_norm () const
    {
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      TrilinosScalar d;
      const int ierr = vector->Norm1 (&d);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return d;
    }



    inline
    Vector::real_type
    Vector::l2_norm () const
    {
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      TrilinosScalar d;
      const int ierr = vector->Norm2 (&d);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return d;
    }



    inline
    Vector::real_type
    Vector::lp_norm (const TrilinosScalar p) const
    {
      Assert (!has_ghost_elements(), ExcGhostsPresent());

      TrilinosScalar norm = 0;
      TrilinosScalar sum=0;
      const size_type n_local = local_size();

      // loop over all the elements because
      // Trilinos does not support lp norms
      for (size_type i=0; i<n_local; i++)
        sum += std::pow(std::fabs((*vector)[0][i]), p);

      norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));

      return norm;
    }



    inline
    Vector::real_type
    Vector::linfty_norm () const
    {
      // while we disallow the other
      // norm operations on ghosted
      // vectors, this particular norm
      // is safe to run even in the
      // presence of ghost elements
      TrilinosScalar d;
      const int ierr = vector->NormInf (&d);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return d;
    }



    inline
    TrilinosScalar
    Vector::add_and_dot (const TrilinosScalar a,
                         const Vector &V,
                         const Vector &W)
    {
      this->add(a, V);
      return *this * W;
    }



    // inline also scalar products, vector
    // additions etc. since they are all
    // representable by a single Trilinos
    // call. This reduces the overhead of the
    // wrapper class.
    inline
    Vector &
    Vector::operator*= (const TrilinosScalar a)
    {
      AssertIsFinite(a);

      const int ierr = vector->Scale(a);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return *this;
    }



    inline
    Vector &
    Vector::operator/= (const TrilinosScalar a)
    {
      AssertIsFinite(a);

      const TrilinosScalar factor = 1./a;

      AssertIsFinite(factor);

      const int ierr = vector->Scale(factor);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return *this;
    }



    inline
    Vector &
    Vector::operator+= (const Vector &v)
    {
      Assert (size() == v.size(),
              ExcDimensionMismatch(size(), v.size()));
      Assert (vector->Map().SameAs(v.vector->Map()),
              ExcDifferentParallelPartitioning());

      const int ierr = vector->Update (1.0, *(v.vector), 1.0);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return *this;
    }



    inline
    Vector &
    Vector::operator-= (const Vector &v)
    {
      Assert (size() == v.size(),
              ExcDimensionMismatch(size(), v.size()));
      Assert (vector->Map().SameAs(v.vector->Map()),
              ExcDifferentParallelPartitioning());

      const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      return *this;
    }



    inline
    void
    Vector::add (const TrilinosScalar s)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      AssertIsFinite(s);

      size_type n_local = local_size();
      for (size_type i=0; i<n_local; i++)
        (*vector)[0][i] += s;
    }



    inline
    void
    Vector::add (const TrilinosScalar  a,
                 const Vector     &v)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (local_size() == v.local_size(),
              ExcDimensionMismatch(local_size(), v.local_size()));

      AssertIsFinite(a);

      const int ierr = vector->Update(a, *(v.vector), 1.);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
    }



    inline
    void
    Vector::add (const TrilinosScalar  a,
                 const Vector     &v,
                 const TrilinosScalar  b,
                 const Vector     &w)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (local_size() == v.local_size(),
              ExcDimensionMismatch(local_size(), v.local_size()));
      Assert (local_size() == w.local_size(),
              ExcDimensionMismatch(local_size(), w.local_size()));

      AssertIsFinite(a);
      AssertIsFinite(b);

      const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);

      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
    }



    inline
    void
    Vector::sadd (const TrilinosScalar  s,
                  const Vector     &v)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (size() == v.size(),
              ExcDimensionMismatch (size(), v.size()));

      AssertIsFinite(s);

      // We assume that the vectors have the same Map
      // if the local size is the same and if the vectors are not ghosted
      if (local_size() == v.local_size() && !v.has_ghost_elements())
        {
          Assert (this->vector->Map().SameAs(v.vector->Map())==true,
                  ExcDifferentParallelPartitioning());
          const int ierr = vector->Update(1., *(v.vector), s);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        }
      else
        {
          (*this) *= s;
          this->add(v, true);
        }
    }



    inline
    void
    Vector::sadd (const TrilinosScalar  s,
                  const TrilinosScalar  a,
                  const Vector     &v)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (size() == v.size(),
              ExcDimensionMismatch (size(), v.size()));
      AssertIsFinite(s);
      AssertIsFinite(a);

      // We assume that the vectors have the same Map
      // if the local size is the same and if the vectors are not ghosted
      if (local_size() == v.local_size() && !v.has_ghost_elements())
        {
          Assert (this->vector->Map().SameAs(v.vector->Map())==true,
                  ExcDifferentParallelPartitioning());
          const int ierr = vector->Update(a, *(v.vector), s);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        }
      else
        {
          (*this) *= s;
          Vector tmp = v;
          tmp *= a;
          this->add(tmp, true);
        }
    }



    inline
    void
    Vector::scale (const Vector &factors)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      Assert (local_size() == factors.local_size(),
              ExcDimensionMismatch(local_size(), factors.local_size()));

      const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
    }



    inline
    void
    Vector::equ (const TrilinosScalar  a,
                 const Vector     &v)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert (!has_ghost_elements(), ExcGhostsPresent());
      AssertIsFinite(a);

      // If we don't have the same map, copy.
      if (vector->Map().SameAs(v.vector->Map())==false)
        {
          this->sadd(0., a, v);
        }
      else
        {
          // Otherwise, just update
          int ierr = vector->Update(a, *v.vector, 0.0);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));

          last_action = Zero;
        }

    }



    inline
    const Epetra_MultiVector &
    Vector::trilinos_vector () const
    {
      return static_cast<const Epetra_MultiVector &>(*vector);
    }



    inline
    Epetra_FEVector &
    Vector::trilinos_vector ()
    {
      return *vector;
    }



    inline
    const Epetra_Map &
    Vector::vector_partitioner () const
    {
      return static_cast<const Epetra_Map &>(vector->Map());
    }



    inline
    const MPI_Comm &
    Vector::get_mpi_communicator () const
    {
      static MPI_Comm comm;

#ifdef DEAL_II_WITH_MPI

      const Epetra_MpiComm *mpi_comm
        = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
      comm = mpi_comm->Comm();

#else

      comm = MPI_COMM_SELF;

#endif

      return comm;
    }

    template <typename number>
    Vector::Vector (const IndexSet               &parallel_partitioner,
                    const ::Vector<number> &v,
                    const MPI_Comm               &communicator)
    {
      *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
                     v);
      owned_elements = parallel_partitioner;
    }



    inline
    Vector &
    Vector::operator= (const TrilinosScalar s)
    {
      AssertIsFinite(s);

      int ierr = vector->PutScalar(s);
      AssertThrow (ierr == 0, ExcTrilinosError(ierr));

      if (nonlocal_vector.get() != nullptr)
        {
          ierr = nonlocal_vector->PutScalar(0.);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        }

      return *this;
    }
  } /* end of namespace MPI */

#endif /* DOXYGEN */

} /* end of namespace TrilinosWrappers */

namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename> class ReinitHelper;

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

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

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



template <>
struct is_serial_vector< TrilinosWrappers::MPI::Vector > : std::false_type
{
};


DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_TRILINOS

/*----------------------------   trilinos_vector.h     ---------------------------*/

#endif // dealii_trilinos_vector_h