trilinos_vector_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2015 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_base_h
#define dealii__trilinos_vector_base_h


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

#ifdef DEAL_II_WITH_TRILINOS

#include <deal.II/base/utilities.h>
#  include <deal.II/base/std_cxx11/shared_ptr.h>
#  include <deal.II/base/subscriptor.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector.h>

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

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#  define TrilinosScalar double
#  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"
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

DEAL_II_NAMESPACE_OPEN

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


namespace TrilinosWrappers
{
  // forward declaration
  class VectorBase;

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

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

    public:

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

      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:
      VectorBase   &vector;

      const size_type  index;

      friend class ::TrilinosWrappers::VectorBase;
    };
  }
  class VectorBase : 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;


    VectorBase ();

    VectorBase (const VectorBase &v);

    virtual ~VectorBase ();

    void clear ();

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

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

    bool is_compressed () const DEAL_II_DEPRECATED;

    VectorBase &
    operator = (const TrilinosScalar s);

    VectorBase &
    operator = (const VectorBase &v);

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

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

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

    TrilinosScalar operator * (const VectorBase &vec) const;

    real_type norm_sqr () const;

    TrilinosScalar mean_value () const;

    TrilinosScalar minimal_value () const DEAL_II_DEPRECATED;

    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 VectorBase    &V,
                                const VectorBase    &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;

    TrilinosScalar el (const size_type index) const DEAL_II_DEPRECATED;

    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);

    VectorBase &operator *= (const TrilinosScalar factor);

    VectorBase &operator /= (const TrilinosScalar factor);

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

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

    void add (const TrilinosScalar s);

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

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

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

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

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

    void sadd (const TrilinosScalar  s,
               const TrilinosScalar  a,
               const VectorBase     &V,
               const TrilinosScalar  b,
               const VectorBase     &W) DEAL_II_DEPRECATED;

    void sadd (const TrilinosScalar  s,
               const TrilinosScalar  a,
               const VectorBase     &V,
               const TrilinosScalar  b,
               const VectorBase     &W,
               const TrilinosScalar  c,
               const VectorBase     &X) DEAL_II_DEPRECATED;

    void scale (const VectorBase &scaling_factors);

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

    void equ (const TrilinosScalar  a,
              const VectorBase     &V,
              const TrilinosScalar  b,
              const VectorBase     &W) DEAL_II_DEPRECATED;

    void ratio (const VectorBase &a,
                const VectorBase &b) DEAL_II_DEPRECATED;



    const Epetra_MultiVector &trilinos_vector () const;

    Epetra_FEVector &trilinos_vector ();

    const Epetra_Map &vector_partitioner () const;

    void print (const char *format = 0) const DEAL_II_DEPRECATED;

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

    void swap (VectorBase &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_cxx11::shared_ptr<Epetra_FEVector> vector;

    std_cxx11::shared_ptr<Epetra_MultiVector> nonlocal_vector;

    friend class internal::VectorReference;
    friend class Vector;
    friend class MPI::Vector;
  };




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

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


#ifndef DOXYGEN

  namespace internal
  {
    inline
    VectorReference::VectorReference (VectorBase      &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
    const 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;
    }
  }



  inline
  bool
  VectorBase::is_compressed () const
  {
    return compressed;
  }



  inline
  bool
  VectorBase::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
  VectorBase::locally_owned_elements() const
  {
    IndexSet is (size());

    // easy case: local range is contiguous
    if (vector->Map().LinearMap())
      {
        const std::pair<size_type, size_type> x = local_range();
        is.add_range (x.first, x.second);
      }
    else if (vector->Map().NumMyElements() > 0)
      {
        const size_type n_indices = vector->Map().NumMyElements();
#ifndef DEAL_II_WITH_64BIT_INDICES
        unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements();
#else
        size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
#endif
        is.add_indices(vector_indices, vector_indices+n_indices);
        is.compress();
      }

    return is;
  }



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



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



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


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



  inline
  void VectorBase::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 VectorBase::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
  VectorBase::iterator
  VectorBase::begin()
  {
    return (*vector)[0];
  }



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



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



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



  inline
  void
  VectorBase::reinit (const VectorBase &v,
                      const bool        omit_zeroing_entries)
  {
    Assert (vector.get() != 0,
            ExcMessage("Vector has not been constructed properly."));

    if (omit_zeroing_entries == false ||
        vector_partitioner().SameAs(v.vector_partitioner())==false)
      vector.reset (new Epetra_FEVector(*v.vector));

    if (v.nonlocal_vector.get() != 0)
      nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
  }



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

    const int ierr = vector->PutScalar(s);

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

    if (nonlocal_vector.get() != 0)
      nonlocal_vector->PutScalar(0.);

    return *this;
  }



  inline
  void
  VectorBase::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[0], &values[0]);
  }



  inline
  void
  VectorBase::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[0], values.begin());
  }



  inline
  void
  VectorBase::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)
      vector->GlobalAssemble(Add);

    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
  VectorBase::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[0], &values[0]);
  }



  inline
  void
  VectorBase::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[0], values.begin());
  }



  inline
  void
  VectorBase::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)
          vector->GlobalAssemble(Insert);
        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() == 0)
          {
            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
  VectorBase::size_type
  VectorBase::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
  VectorBase::size_type
  VectorBase::local_size () const
  {
    return (size_type) vector->Map().NumMyElements();
  }



  inline
  std::pair<VectorBase::size_type, VectorBase::size_type>
  VectorBase::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
  VectorBase::operator * (const VectorBase &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
  VectorBase::real_type
  VectorBase::norm_sqr () const
  {
    const TrilinosScalar d = l2_norm();
    return d*d;
  }



  inline
  TrilinosScalar
  VectorBase::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
  VectorBase::minimal_value () const
  {
    return min();
  }



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

    return min_value;
  }



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

    return max_value;
  }



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

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

    return d;
  }



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

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

    return d;
  }



  inline
  VectorBase::real_type
  VectorBase::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
  VectorBase::real_type
  VectorBase::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
  VectorBase::add_and_dot (const TrilinosScalar a,
                           const VectorBase &V,
                           const VectorBase &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
  VectorBase &
  VectorBase::operator *= (const TrilinosScalar a)
  {
    AssertIsFinite(a);

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

    return *this;
  }



  inline
  VectorBase &
  VectorBase::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
  VectorBase &
  VectorBase::operator += (const VectorBase &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
  VectorBase &
  VectorBase::operator -= (const VectorBase &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
  VectorBase::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
  VectorBase::add (const TrilinosScalar  a,
                   const VectorBase     &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
  VectorBase::add (const TrilinosScalar  a,
                   const VectorBase     &v,
                   const TrilinosScalar  b,
                   const VectorBase     &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
  VectorBase::sadd (const TrilinosScalar  s,
                    const VectorBase     &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
  VectorBase::sadd (const TrilinosScalar  s,
                    const TrilinosScalar  a,
                    const VectorBase     &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;
        VectorBase tmp = v;
        tmp *= a;
        this->add(tmp, true);
      }
  }



  inline
  void
  VectorBase::sadd (const TrilinosScalar  s,
                    const TrilinosScalar  a,
                    const VectorBase     &v,
                    const TrilinosScalar  b,
                    const VectorBase     &w)
  {
    // 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()));
    Assert (size() == w.size(),
            ExcDimensionMismatch (size(), w.size()));
    AssertIsFinite(s);
    AssertIsFinite(a);
    AssertIsFinite(b);

    // 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() &&
        local_size() == w.local_size() && !w.has_ghost_elements())
      {
        Assert (this->vector->Map().SameAs(v.vector->Map())==true,
                ExcDifferentParallelPartitioning());
        Assert (this->vector->Map().SameAs(w.vector->Map())==true,
                ExcDifferentParallelPartitioning());
        const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
      }
    else
      {
        this->sadd( s, a, v);
        this->sadd(1., b, w);
      }
  }



  inline
  void
  VectorBase::sadd (const TrilinosScalar  s,
                    const TrilinosScalar  a,
                    const VectorBase     &v,
                    const TrilinosScalar  b,
                    const VectorBase     &w,
                    const TrilinosScalar  c,
                    const VectorBase     &x)
  {
    // 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()));
    Assert (size() == w.size(),
            ExcDimensionMismatch (size(), w.size()));
    Assert (size() == x.size(),
            ExcDimensionMismatch (size(), x.size()));
    AssertIsFinite(s);
    AssertIsFinite(a);
    AssertIsFinite(b);
    AssertIsFinite(c);

    // 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() &&
        local_size() == w.local_size() && !w.has_ghost_elements() &&
        local_size() == x.local_size() && !x.has_ghost_elements())
      {
        Assert (this->vector->Map().SameAs(v.vector->Map())==true,
                ExcDifferentParallelPartitioning());
        Assert (this->vector->Map().SameAs(w.vector->Map())==true,
                ExcDifferentParallelPartitioning());
        Assert (this->vector->Map().SameAs(x.vector->Map())==true,
                ExcDifferentParallelPartitioning());

        // Update member can only
        // input two other vectors so
        // do it in two steps
        const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));

        const int jerr = vector->Update(c, *(x.vector), 1.);
        Assert (jerr == 0, ExcTrilinosError(jerr));
        (void)jerr; // removes -Wunused-parameter warning in optimized mode
      }
    else
      {
        this->sadd( s, a, v);
        this->sadd(1., b, w);
        this->sadd(1., c, x);
      }
  }



  inline
  void
  VectorBase::scale (const VectorBase &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
  VectorBase::equ (const TrilinosScalar  a,
                   const VectorBase     &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
  void
  VectorBase::ratio (const VectorBase &v,
                     const VectorBase &w)
  {
    Assert (v.local_size() == w.local_size(),
            ExcDimensionMismatch (v.local_size(), w.local_size()));
    Assert (local_size() == w.local_size(),
            ExcDimensionMismatch (local_size(), w.local_size()));

    const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
                                                *(v.vector), 0.0);

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



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



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



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



  inline
  const MPI_Comm &
  VectorBase::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;
  }



#endif // DOXYGEN

}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_TRILINOS

/*----------------------------   trilinos_vector_base.h     ---------------------------*/

#endif
/*----------------------------   trilinos_vector_base.h     ---------------------------*/