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

#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/index_set.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/base/subscriptor.h>
#  include <deal.II/base/utilities.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 <Epetra_ConfigDefs.h>

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

DEAL_II_NAMESPACE_OPEN


namespace TrilinosWrappers
{
  class SparseMatrix;

  namespace internal
  {
    using size_type = ::types::global_dof_index;

    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 internal
#  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:
      using value_type      = TrilinosScalar;
      using real_type       = TrilinosScalar;
      using size_type       = ::types::global_dof_index;
      using iterator        = value_type *;
      using const_iterator  = const value_type *;
      using reference       = internal::VectorReference;
      using const_reference = const internal::VectorReference;

      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() override = 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);
    }
  } // namespace MPI

#  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 internal

  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 LinearOperatorImplementation
} /* 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