]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove TrilinosWrappers::VectorBase
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 30 May 2017 22:25:23 +0000 (00:25 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Sat, 3 Jun 2017 16:56:51 +0000 (18:56 +0200)
include/deal.II/lac/trilinos_vector.h
include/deal.II/lac/trilinos_vector_base.h
source/lac/CMakeLists.txt
source/lac/trilinos_vector.cc
source/lac/trilinos_vector.inst.in [moved from source/lac/trilinos_vector_base.inst.in with 89% similarity]
source/lac/trilinos_vector_base.cc [deleted file]

index 60ffa95c6da12ab4cdd1951873cc13f23ed65a48..8a8bb874e81c211387bcbb537f7b83a13f75ae7c 100644 (file)
 #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/trilinos_vector_base.h>
 #  include <deal.II/lac/vector_type_traits.h>
 
+#  include <vector>
+#  include <utility>
+#  include <memory>
+
 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+#  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_ENABLE_EXTRA_DIAGNOSTICS
 
-#include <memory>
-
 DEAL_II_NAMESPACE_OPEN
 
 
-// forward declaration
-template <typename> class Vector;
-
 /**
  * @addtogroup TrilinosWrappers
  * @{
  */
+
+/**
+ * A namespace in which wrapper classes for Trilinos objects reside.
+ *
+ * @ingroup TrilinosWrappers
+ */
 namespace TrilinosWrappers
 {
   class SparseMatrix;
 
+  /**
+   * @cond internal
+   */
+
+  /**
+   * A namespace for internal implementation details of the TrilinosWrapper
+   * members.
+   *
+   * @ingroup TrilinosWrappers
+   */
+  namespace internal
+  {
+    /**
+     * Declare type for container size.
+     */
+    typedef dealii::types::global_dof_index size_type;
+
+    /**
+     * This class implements a wrapper for accessing the Trilinos vector in
+     * the same way as we access deal.II objects: it is initialized with a
+     * vector and an element within it, and has a conversion operatorto
+     * extract the scalar value of this element. It also has a variety of
+     * assignment operatorfor writing to this one element.
+     *
+     * @ingroup TrilinosWrappers
+     */
+    class VectorReference
+    {
+    private:
+      /**
+       * Constructor. It is made private so as to only allow the actual vector
+       * class to create it.
+       */
+      VectorReference (MPI::Vector     &vector,
+                       const size_type  index);
+
+    public:
+
+      /**
+       * This looks like a copy operator, but does something different than
+       * usual. In particular, it does not copy the member variables of this
+       * reference. Rather, it handles the situation where we have two vectors
+       * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
+       * both left and right hand side of the assignment have data type
+       * VectorReference, but what we really mean is to assign the vector
+       * elements represented by the two references. This operatorimplements
+       * this operation. Note also that this allows us to make the assignment
+       * operatorconst.
+       */
+      const VectorReference &
+      operator= (const VectorReference &r) const;
+
+      /**
+       * Same as above but for non-const reference objects.
+       */
+      VectorReference &
+      operator= (const VectorReference &r);
+
+      /**
+       * Set the referenced element of the vector to <tt>s</tt>.
+       */
+      const VectorReference &
+      operator= (const TrilinosScalar &s) const;
+
+      /**
+       * Add <tt>s</tt> to the referenced element of the vector->
+       */
+      const VectorReference &
+      operator+= (const TrilinosScalar &s) const;
+
+      /**
+       * Subtract <tt>s</tt> from the referenced element of the vector->
+       */
+      const VectorReference &
+      operator-= (const TrilinosScalar &s) const;
+
+      /**
+       * Multiply the referenced element of the vector by <tt>s</tt>.
+       */
+      const VectorReference &
+      operator*= (const TrilinosScalar &s) const;
+
+      /**
+       * Divide the referenced element of the vector by <tt>s</tt>.
+       */
+      const VectorReference &
+      operator/= (const TrilinosScalar &s) const;
+
+      /**
+       * Convert the reference to an actual value, i.e. return the value of
+       * the referenced element of the vector.
+       */
+      operator TrilinosScalar () const;
+
+      /**
+       * Exception
+       */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occurred while calling a Trilinos function");
+
+    private:
+      /**
+       * Point to the vector we are referencing.
+       */
+      MPI::Vector   &vector;
+
+      /**
+       * Index of the referenced element of the vector.
+       */
+      const size_type  index;
+
+      /**
+       * Make the vector class a friend, so that it can create objects of the
+       * present type.
+       */
+      friend class ::dealii::TrilinosWrappers::MPI::Vector;
+    };
+  }
+  /**
+   * @endcond
+   */
+
   namespace
   {
 #ifndef DEAL_II_WITH_64BIT_INDICES
@@ -75,11 +211,9 @@ namespace TrilinosWrappers
 
   /**
    * Namespace for Trilinos vector classes that work in parallel over MPI.
-   * This namespace is restricted to vectors only, whereas matrices are always
-   * MPI based when run on more than one processor.
    *
    * @ingroup TrilinosWrappers
-   * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+   * @author Martin Kronbichler, Wolfgang Bangerth, Daniel Arndt, 2008, 2017
    */
   namespace MPI
   {
@@ -87,9 +221,19 @@ namespace TrilinosWrappers
 
     /**
      * This class implements a wrapper to use the Trilinos distributed vector
-     * class Epetra_FEVector. This class is derived from the
-     * TrilinosWrappers::VectorBase class and provides all functionality
-     * included there.
+     * class Epetra_FEVector, the (parallel) partitioning of which
+     * is governed by an Epetra_Map.
+     * The Epetra_FEVector is precisely the kind of vector
+     * we deal with all the time - we probably get it from some assembly
+     * process, where also entries not locally owned might need to written and
+     * hence need to be forwarded to the owner. The only requirement for this
+     * class to work is that Trilinos is installed with the same compiler as is
+     * used for compilation of deal.II.
+     *
+     * The interface of this class is modeled after the existing Vector class in
+     * deal.II. It has almost the same member functions, and is often
+     * exchangeable. However, since Trilinos only supports a single scalar type
+     * (double), it is not templated, and only works with that type.
      *
      * Note that Trilinos only guarantees that operations do what you expect
      * if the function @p GlobalAssemble has been called after vector assembly
@@ -159,7 +303,7 @@ namespace TrilinosWrappers
      * code that may for example look like this:
      *
      * @code
-     * TrilinosWrappers::VectorBase vector;
+     * TrilinosWrappers::MPI::Vector vector;
      * // do some write operations on the vector
      * for (size_type i=0; i<vector->size(); ++i)
      *   vector(i) = i;
@@ -250,18 +394,27 @@ namespace TrilinosWrappers
      *
      * @ingroup TrilinosWrappers
      * @ingroup Vectors
-     * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009
+     * @author Martin Kronbichler, Wolfgang Bangerth, Daniel Arndt,
+     *         2008, 2009, 2017
      */
-    class Vector : public VectorBase
+    class Vector : public Subscriptor
     {
     public:
       /**
-       * Declare type for container size.
+       * Declare some of the standard types used in all containers. These types
+       * parallel those in the <tt>C</tt> standard libraries
+       * <tt>vector<...></tt> class.
        */
+      typedef TrilinosScalar                  value_type;
+      typedef TrilinosScalar                  real_type;
       typedef dealii::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;
 
       /**
-       * @name Basic constructors and initialization.
+       * @name 1: Basic Object-handling
        */
       //@{
       /**
@@ -276,116 +429,6 @@ namespace TrilinosWrappers
        */
       Vector (const Vector &v);
 
-      /**
-       * Move constructor. Creates a new vector by stealing the internal data
-       * of the vector @p v.
-       */
-      Vector (Vector &&v);
-
-      /**
-       * Destructor.
-       */
-      ~Vector ();
-
-      /**
-       * Reinit functionality. This function sets the calling vector to the
-       * dimension and the parallel distribution of the input vector, but does
-       * not copy the elements in <tt>v</tt>. If <tt>omit_zeroing_entries</tt>
-       * is not <tt>true</tt>, the elements in the vector are initialized with
-       * zero. If it is set to <tt>true</tt>, the vector entries are in an
-       * unspecified state and the user has to set all elements. In the
-       * current implementation, this method does not touch the vector entries
-       * in case the vector layout is unchanged from before, otherwise entries
-       * are set to zero.  Note that this behavior might change between
-       * releases without notification.
-       *
-       * This function has a third argument, <tt>allow_different_maps</tt>,
-       * that allows for an exchange of data between two equal-sized vectors
-       * (but being distributed differently among the processors). A trivial
-       * application of this function is to generate a replication of a whole
-       * vector on each machine, when the calling vector is built according to
-       * the localized vector class TrilinosWrappers::VectorBase, and <tt>v</tt>
-       * is a distributed vector. In this case, the variable
-       * <tt>omit_zeroing_entries</tt> needs to be set to <tt>false</tt>,
-       * since it does not make sense to exchange data between differently
-       * parallelized vectors without touching the elements.
-       */
-      void reinit (const VectorBase &v,
-                   const bool        omit_zeroing_entries = false,
-                   const bool        allow_different_maps = false);
-
-      /**
-       * Create vector by merging components from a block vector.
-       */
-      void reinit (const BlockVector &v,
-                   const bool         import_data = false);
-
-      /**
-       * Set all components of the vector to the given number @p s. Simply
-       * pass this down to the base class, but we still need to declare this
-       * function to make the example given in the discussion about making the
-       * constructor explicit work.
-       */
-      Vector &operator= (const TrilinosScalar s);
-
-      /**
-       * Copy the given vector. Resize the present vector if necessary. In
-       * this case, also the Epetra_Map that designs the parallel partitioning
-       * is taken from the input vector.
-       */
-      Vector &operator= (const Vector &v);
-
-      /**
-       * Move the given vector. This operator replaces the present vector with
-       * @p v by efficiently swapping the internal data structures.
-       */
-      Vector &operator= (Vector &&v);
-
-      /**
-       * Copy operator from a given localized vector (present on all
-       * processes) in TrilinosWrappers format to the current distributed
-       * vector. This function assumes that the calling vector (left hand
-       * object) already is of the same size as the right hand side vector.
-       * Otherwise, an exception will be thrown.
-       */
-      Vector &operator= (const ::dealii::TrilinosWrappers::VectorBase &v);
-
-      /**
-       * Another copy function. This one takes a deal.II vector and copies it
-       * into a TrilinosWrapper vector. Note that since we do not provide any
-       * Epetra_map that tells about the partitioning of the vector among the
-       * MPI processes, the size of the TrilinosWrapper vector has to be the
-       * same as the size of the input vector. In order to change the map, use
-       * the reinit(const Epetra_Map &input_map) function.
-       */
-      template <typename Number>
-      Vector &operator= (const ::dealii::Vector<Number> &v);
-
-      /**
-       * This reinit function is meant to be used for parallel calculations
-       * where some non-local data has to be used. The typical situation where
-       * one needs this function is the call of the
-       * FEValues<dim>::get_function_values function (or of some derivatives)
-       * in parallel. Since it is usually faster to retrieve the data in
-       * advance, this function can be called before the assembly forks out to
-       * the different processors. What this function does is the following:
-       * It takes the information in the columns of the given matrix and looks
-       * which data couples between the different processors. That data is
-       * then queried from the input vector. Note that you should not write to
-       * the resulting vector any more, since the some data can be stored
-       * several times on different processors, leading to unpredictable
-       * results. In particular, such a vector cannot be used for matrix-
-       * vector products as for example done during the solution of linear
-       * systems.
-       */
-      void import_nonlocal_data_for_fe
-      (const dealii::TrilinosWrappers::SparseMatrix &matrix,
-       const Vector                                 &vector);
-//@}
-      /**
-       * @name Initialization with an IndexSet
-       */
-//@{
       /**
        * This constructor takes an IndexSet that defines how to distribute the
        * individual components among the MPI processors. Since it also
@@ -437,7 +480,7 @@ namespace TrilinosWrappers
        * @ref GlossGhostedVector "vectors with ghost elements"
        */
       Vector (const IndexSet   &parallel_partitioning,
-              const VectorBase &v,
+              const Vector &v,
               const MPI_Comm   &communicator = MPI_COMM_WORLD);
 
       /**
@@ -457,6 +500,50 @@ namespace TrilinosWrappers
               const dealii::Vector<Number> &v,
               const MPI_Comm               &communicator = MPI_COMM_WORLD);
 
+      /**
+       * Move constructor. Creates a new vector by stealing the internal data
+       * of the vector @p v.
+       */
+      Vector (Vector &&v);
+
+      /**
+       * Destructor.
+       */
+      ~Vector ();
+
+      /**
+       * Release all memory and return to a state just like after having called
+       * the default constructor.
+       */
+      void clear ();
+
+      /**
+       * Reinit functionality. This function sets the calling vector to the
+       * dimension and the parallel distribution of the input vector, but does
+       * not copy the elements in <tt>v</tt>. If <tt>omit_zeroing_entries</tt>
+       * is not <tt>true</tt>, the elements in the vector are initialized with
+       * zero. If it is set to <tt>true</tt>, the vector entries are in an
+       * unspecified state and the user has to set all elements. In the
+       * current implementation, this method does not touch the vector entries
+       * in case the vector layout is unchanged from before, otherwise entries
+       * are set to zero.  Note that this behavior might change between
+       * releases without notification.
+       *
+       * This function has a third argument, <tt>allow_different_maps</tt>,
+       * that allows for an exchange of data between two equal-sized vectors
+       * (but being distributed differently among the processors). A trivial
+       * application of this function is to generate a replication of a whole
+       * vector on each machine, when the calling vector is built with a map
+       * consisting of all indices on each process, and <tt>v</tt>
+       * is a distributed vector. In this case, the variable
+       * <tt>omit_zeroing_entries</tt> needs to be set to <tt>false</tt>,
+       * since it does not make sense to exchange data between differently
+       * parallelized vectors without touching the elements.
+       */
+      void reinit (const Vector &v,
+                   const bool        omit_zeroing_entries = false,
+                   const bool        allow_different_maps = false);
+
       /**
        * Reinit functionality. This function destroys the old vector content
        * and generates a new one based on the input partitioning.  The flag
@@ -467,7 +554,6 @@ namespace TrilinosWrappers
        * the entries to zero, but this might change between releases without
        * notification.
        *
-       *
        * Depending on whether the @p parallel_partitioning argument uniquely
        * subdivides elements among processors or not, the resulting vector may
        * or may not have ghost elements. See the general documentation of this
@@ -513,85 +599,1562 @@ namespace TrilinosWrappers
                    const IndexSet &ghost_entries,
                    const MPI_Comm &communicator = MPI_COMM_WORLD,
                    const bool      vector_writable = false);
-//@}
-    };
-
 
+      /**
+       * Create vector by merging components from a block vector.
+       */
+      void reinit (const BlockVector &v,
+                   const bool         import_data = false);
 
+      /**
+       * Compress the underlying representation of the Trilinos object, i.e.
+       * flush the buffers of the vector object if it has any. This function is
+       * necessary after writing into a vector element-by-element and before
+       * anything else can be done on it.
+       *
+       * The (defaulted) argument can be used to specify the compress mode
+       * (<code>Add</code> or <code>Insert</code>) in case the vector has not
+       * been written to since the last time this function was called. The
+       * argument is ignored if the vector has been added or written to since
+       * the last time compress() was called.
+       *
+       * See
+       * @ref GlossCompress "Compressing distributed objects"
+       * for more information.
+       */
+      void compress (::dealii::VectorOperation::values operation);
 
-// ------------------- inline and template functions --------------
+      /**
+       * Set all components of the vector to the given number @p s. Simply
+       * pass this down to the base class, but we still need to declare this
+       * function to make the example given in the discussion about making the
+       * constructor explicit work.
+       * the constructor explicit work.
+       *
+       * Since the semantics of assigning a scalar to a vector are not
+       * immediately clear, this operatorcan only be used if you want
+       * to set the entire vector to zero. This allows the intuitive notation
+       * <tt>v=0</tt>.
+       */
+      Vector &operator= (const TrilinosScalar s);
 
+      /**
+       * Copy the given vector. Resize the present vector if necessary. In
+       * this case, also the Epetra_Map that designs the parallel partitioning
+       * is taken from the input vector.
+       */
+      Vector &operator= (const Vector &v);
 
-    /**
-     * Global function @p swap which overloads the default implementation of
-     * the C++ standard library which uses a temporary object. The function
-     * simply exchanges the data of the two vectors.
-     *
-     * @relates TrilinosWrappers::MPI::Vector
-     * @author Martin Kronbichler, Wolfgang Bangerth, 2008
-     */
-    inline
-    void swap (Vector &u, Vector &v)
-    {
-      u.swap (v);
-    }
+      /**
+       * Move the given vector. This operatorreplaces the present vector with
+       * @p v by efficiently swapping the internal data structures.
+       */
+      Vector &operator= (Vector &&v);
 
+      /**
+       * Another copy function. This one takes a deal.II vector and copies it
+       * into a TrilinosWrapper vector. Note that since we do not provide any
+       * Epetra_map that tells about the partitioning of the vector among the
+       * MPI processes, the size of the TrilinosWrapper vector has to be the
+       * same as the size of the input vector.
+       */
+      template <typename Number>
+      Vector &operator= (const ::dealii::Vector<Number> &v);
 
-#ifndef DOXYGEN
+      /**
+       * This reinit function is meant to be used for parallel calculations
+       * where some non-local data has to be used. The typical situation where
+       * one needs this function is the call of the
+       * FEValues<dim>::get_function_values function (or of some derivatives)
+       * in parallel. Since it is usually faster to retrieve the data in
+       * advance, this function can be called before the assembly forks out to
+       * the different processors. What this function does is the following:
+       * It takes the information in the columns of the given matrix and looks
+       * which data couples between the different processors. That data is
+       * then queried from the input vector. Note that you should not write to
+       * the resulting vector any more, since the some data can be stored
+       * several times on different processors, leading to unpredictable
+       * results. In particular, such a vector cannot be used for matrix-
+       * vector products as for example done during the solution of linear
+       * systems.
+       */
+      void import_nonlocal_data_for_fe
+      (const dealii::TrilinosWrappers::SparseMatrix &matrix,
+       const Vector                                 &vector);
+
+      /**
+       * Test for equality. This function assumes that the present vector and
+       * the one to compare with have the same size already, since comparing
+       * vectors of different sizes makes not much sense anyway.
+       */
+      bool operator== (const Vector &v) const;
+
+      /**
+       * Test for inequality. This function assumes that the present vector and
+       * the one to compare with have the same size already, since comparing
+       * vectors of different sizes makes not much sense anyway.
+       */
+      bool operator!= (const Vector &v) const;
+
+      /**
+       * Return the global dimension of the vector.
+       */
+      size_type size () const;
+
+      /**
+       * Return the local dimension of the vector, i.e. the number of elements
+       * stored on the present MPI process. For sequential vectors, this number
+       * is the same as size(), but for parallel vectors it may be smaller.
+       *
+       * To figure out which elements exactly are stored locally, use
+       * local_range().
+       *
+       * If the vector contains ghost elements, they are included in this
+       * number.
+       */
+      size_type local_size () const;
+
+      /**
+       * Return a pair of indices indicating which elements of this vector are
+       * stored locally. The first number is the index of the first element
+       * stored, the second the index of the one past the last one that is
+       * stored locally. If this is a sequential vector, then the result will be
+       * the pair <code>(0,N)</code>, otherwise it will be a pair
+       * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
+       * <code>i</code> is the first element of the vector stored on this
+       * processor, corresponding to the half open interval $[i,i+n)$
+       *
+       * @note The description above is true most of the time, but not always.
+       * In particular, Trilinos vectors need not store contiguous ranges of
+       * elements such as $[i,i+n)$. Rather, it can store vectors where the
+       * elements are distributed in an arbitrary way across all processors and
+       * each processor simply stores a particular subset, not necessarily
+       * contiguous. In this case, this function clearly makes no sense since it
+       * could, at best, return a range that includes all elements that are
+       * stored locally. Thus, the function only succeeds if the locally stored
+       * range is indeed contiguous. It will trigger an assertion if the local
+       * portion of the vector is not contiguous.
+       */
+      std::pair<size_type, size_type> local_range () const;
+
+      /**
+       * Return whether @p index is in the local range or not, see also
+       * local_range().
+       *
+       * @note The same limitation for the applicability of this function
+       * applies as listed in the documentation of local_range().
+       */
+      bool in_local_range (const size_type index) const;
+
+      /**
+       * Return an index set that describes which elements of this vector are
+       * owned by the current processor. Note that this index set does not
+       * include elements this vector may store locally as ghost elements but
+       * that are in fact owned by another processor. As a consequence, the
+       * index sets returned on different processors if this is a distributed
+       * vector will form disjoint sets that add up to the complete index set.
+       * Obviously, if a vector is created on only one processor, then the
+       * result would satisfy
+       * @code
+       *   vec.locally_owned_elements() == complete_index_set (vec.size())
+       * @endcode
+       */
+      IndexSet locally_owned_elements () const;
+
+      /**
+       * Return if the vector contains ghost elements. This answer is true if
+       * there are ghost elements on at least one process.
+       *
+       * @see
+       * @ref GlossGhostedVector "vectors with ghost elements"
+       */
+      bool has_ghost_elements() const;
+
+      /**
+       * This function only exists for compatibility with the @p
+       * LinearAlgebra::distributed::Vector class and does nothing: this class
+       * implements ghost value updates in a different way that is a better fit
+       * with the underlying Trilinos vector object.
+       */
+      void update_ghost_values () const;
+
+      /**
+       * Return the scalar (inner) product of two vectors. The vectors must have
+       * the same size.
+       */
+      TrilinosScalar operator* (const Vector &vec) const;
+
+      /**
+       * Return square of the $l_2$-norm.
+       */
+      real_type norm_sqr () const;
+
+      /**
+       * Mean value of the elements of this vector.
+       */
+      TrilinosScalar mean_value () const;
+
+      /**
+       * Compute the minimal value of the elements of this vector.
+       */
+      TrilinosScalar min () const;
+
+      /**
+       * Compute the maximal value of the elements of this vector.
+       */
+      TrilinosScalar max () const;
+
+      /**
+       * $l_1$-norm of the vector.  The sum of the absolute values.
+       */
+      real_type l1_norm () const;
+
+      /**
+       * $l_2$-norm of the vector.  The square root of the sum of the squares of
+       * the elements.
+       */
+      real_type l2_norm () const;
+
+      /**
+       * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
+       * <i>p</i>th powers of the absolute values of the elements.
+       */
+      real_type lp_norm (const TrilinosScalar p) const;
+
+      /**
+       * Maximum absolute value of the elements.
+       */
+      real_type linfty_norm () const;
+
+      /**
+       * Performs a combined operation of a vector addition and a subsequent
+       * inner product, returning the value of the inner product. In other
+       * words, the result of this function is the same as if the user called
+       * @code
+       * this->add(a, V);
+       * return_value = *this * W;
+       * @endcode
+       *
+       * The reason this function exists is for compatibility with deal.II's own
+       * vector classes which can implement this functionality with less memory
+       * transfer. However, for Trilinos vectors such a combined operation is
+       * not natively supported and thus the cost is completely equivalent as
+       * calling the two methods separately.
+       */
+      TrilinosScalar add_and_dot (const TrilinosScalar a,
+                                  const Vector    &V,
+                                  const Vector    &W);
+
+      /**
+       * Return whether the vector contains only elements with value zero. This
+       * is a collective operation. This function is expensive, because
+       * potentially all elements have to be checked.
+       */
+      bool all_zero () const;
+
+      /**
+       * Return @p true if the vector has no negative entries, i.e. all entries
+       * are zero or positive. This function is used, for example, to check
+       * whether refinement indicators are really all positive (or zero).
+       */
+      bool is_non_negative () const;
+      //@}
+
+
+      /**
+       * @name 2: Data-Access
+       */
+      //@{
+
+      /**
+       * Provide access to a given element, both read and write.
+       *
+       * When using a vector distributed with MPI, this operation only makes
+       * sense for elements that are actually present on the calling processor.
+       * Otherwise, an exception is thrown. This is different from the
+       * <code>el()</code> function below that always succeeds (but returns zero
+       * on non-local elements).
+       */
+      reference
+      operator() (const size_type index);
+
+      /**
+       * Provide read-only access to an element.
+       *
+       * When using a vector distributed with MPI, this operation only makes
+       * sense for elements that are actually present on the calling processor.
+       * Otherwise, an exception is thrown. This is different from the
+       * <code>el()</code> function below that always succeeds (but returns zero
+       * on non-local elements).
+       */
+      TrilinosScalar
+      operator() (const size_type index) const;
+
+      /**
+       * Provide access to a given element, both read and write.
+       *
+       * Exactly the same as operator().
+       */
+      reference
+      operator[] (const size_type index);
+
+      /**
+       * Provide read-only access to an element.
+       *
+       * Exactly the same as operator().
+       */
+      TrilinosScalar
+      operator[] (const size_type index) const;
+
+      /**
+       * A collective get operation: instead of getting individual elements of a
+       * vector, this function allows to get a whole set of elements at once.
+       * The indices of the elements to be read are stated in the first
+       * argument, the corresponding values are returned in the second.
+       */
+      void extract_subvector_to (const std::vector<size_type> &indices,
+                                 std::vector<TrilinosScalar> &values) const;
+
+      /**
+       * Just as the above, but with pointers.  Useful in minimizing copying of
+       * data around.
+       */
+      template <typename ForwardIterator, typename OutputIterator>
+      void extract_subvector_to (ForwardIterator          indices_begin,
+                                 const ForwardIterator    indices_end,
+                                 OutputIterator           values_begin) const;
+
+      /**
+       * Make the Vector class a bit like the <tt>vector<></tt> class of the C++
+       * standard library by returning iterators to the start and end of the
+       * locally owned elements of this vector. The ordering of local elements
+       * corresponds to the one given by the global indices in case the vector
+       * is constructed from an IndexSet or other methods in deal.II (note that
+       * an Epetra_Map can contain elements in arbitrary orders, though).
+       *
+       * It holds that end() - begin() == local_size().
+       */
+      iterator begin ();
+
+      /**
+       * Return constant iterator to the start of the locally owned elements of
+       * the vector.
+       */
+      const_iterator begin () const;
+
+      /**
+       * Return an iterator pointing to the element past the end of the array of
+       * locally owned entries.
+       */
+      iterator end ();
+
+      /**
+       * Return a constant iterator pointing to the element past the end of the
+       * array of the locally owned entries.
+       */
+      const_iterator end () const;
+
+      //@}
+
+
+      /**
+       * @name 3: Modification of vectors
+       */
+      //@{
+
+      /**
+       * A collective set operation: instead of setting individual elements of a
+       * vector, this function allows to set a whole set of elements at once.
+       * The indices of the elements to be set are stated in the first argument,
+       * the corresponding values in the second.
+       */
+      void set (const std::vector<size_type>    &indices,
+                const std::vector<TrilinosScalar>  &values);
+
+      /**
+       * This is a second collective set operation. As a difference, this
+       * function takes a deal.II vector of values.
+       */
+      void set (const std::vector<size_type>        &indices,
+                const ::dealii::Vector<TrilinosScalar> &values);
+
+      /**
+       * This collective set operation is of lower level and can handle anything
+       * else &mdash; the only thing you have to provide is an address where all
+       * the indices are stored and the number of elements to be set.
+       */
+      void set (const size_type       n_elements,
+                const size_type      *indices,
+                const TrilinosScalar *values);
+
+      /**
+       * A collective add operation: This function adds a whole set of values
+       * stored in @p values to the vector components specified by @p indices.
+       */
+      void add (const std::vector<size_type>      &indices,
+                const std::vector<TrilinosScalar> &values);
+
+      /**
+       * This is a second collective add operation. As a difference, this
+       * function takes a deal.II vector of values.
+       */
+      void add (const std::vector<size_type>           &indices,
+                const ::dealii::Vector<TrilinosScalar> &values);
+
+      /**
+       * Take an address where <tt>n_elements</tt> are stored contiguously and
+       * add them into the vector. Handles all cases which are not covered by
+       * the other two <tt>add()</tt> functions above.
+       */
+      void add (const size_type       n_elements,
+                const size_type      *indices,
+                const TrilinosScalar *values);
+
+      /**
+       * Multiply the entire vector by a fixed factor.
+       */
+      Vector &operator*= (const TrilinosScalar factor);
+
+      /**
+       * Divide the entire vector by a fixed factor.
+       */
+      Vector &operator/= (const TrilinosScalar factor);
+
+      /**
+       * Add the given vector to the present one.
+       */
+      Vector &operator+= (const Vector &V);
+
+      /**
+       * Subtract the given vector from the present one.
+       */
+      Vector &operator-= (const Vector &V);
+
+      /**
+       * Addition of @p s to all components. Note that @p s is a scalar and not
+       * a vector.
+       */
+      void add (const TrilinosScalar s);
+
+      /**
+       * Simple vector addition, equal to the <tt>operator+=</tt>.
+       *
+       * Though, if the second argument <tt>allow_different_maps</tt> is set,
+       * then it is possible to add data from a vector that uses a different
+       * map, i.e., a vector whose elements are split across processors
+       * differently. This may include vectors with ghost elements, for example.
+       * In general, however, adding vectors with a different element-to-
+       * processor map requires communicating data among processors and,
+       * consequently, is a slower operation than when using vectors using the
+       * same map.
+       */
+      void add (const Vector &V,
+                const bool        allow_different_maps = false);
+
+      /**
+       * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+       */
+      void add (const TrilinosScalar  a,
+                const Vector     &V);
+
+      /**
+       * Multiple addition of scaled vectors, i.e. <tt>*this += a*V + b*W</tt>.
+       */
+      void add (const TrilinosScalar  a,
+                const Vector     &V,
+                const TrilinosScalar  b,
+                const Vector     &W);
+
+      /**
+       * Scaling and simple vector addition, i.e.  <tt>*this = s*(*this) +
+       * V</tt>.
+       */
+      void sadd (const TrilinosScalar  s,
+                 const Vector     &V);
+
+      /**
+       * Scaling and simple addition, i.e.  <tt>*this = s*(*this) + a*V</tt>.
+       */
+      void sadd (const TrilinosScalar  s,
+                 const TrilinosScalar  a,
+                 const Vector     &V);
+
+      /**
+       * Scale each element of this vector by the corresponding element in the
+       * argument. This function is mostly meant to simulate multiplication (and
+       * immediate re-assignment) by a diagonal scaling matrix.
+       */
+      void scale (const Vector &scaling_factors);
+
+      /**
+       * Assignment <tt>*this = a*V</tt>.
+       */
+      void equ (const TrilinosScalar  a,
+                const Vector     &V);
+      //@}
+
+      /**
+       * @name 4: Mixed stuff
+       */
+      //@{
+
+      /**
+       * Return a const reference to the underlying Trilinos Epetra_MultiVector
+       * class.
+       */
+      const Epetra_MultiVector &trilinos_vector () const;
+
+      /**
+       * Return a (modifyable) reference to the underlying Trilinos
+       * Epetra_FEVector class.
+       */
+      Epetra_FEVector &trilinos_vector ();
+
+      /**
+       * Return a const reference to the underlying Trilinos Epetra_Map that
+       * sets the parallel partitioning of the vector.
+       */
+      const Epetra_Map &vector_partitioner () const;
+
+      /**
+       * Print to a stream. @p precision denotes the desired precision with
+       * which values shall be printed, @p scientific whether scientific
+       * notation shall be used. If @p across is @p true then the vector is
+       * printed in a line, while if @p false then the elements are printed on a
+       * separate line each.
+       */
+      void print (std::ostream       &out,
+                  const unsigned int  precision  = 3,
+                  const bool          scientific = true,
+                  const bool          across     = true) const;
+
+      /**
+       * Swap the contents of this vector and the other vector @p v. One could
+       * do this operation with a temporary variable and copying over the data
+       * elements, but this function is significantly more efficient since it
+       * only swaps the pointers to the data of the two vectors and therefore
+       * does not need to allocate temporary storage and move data around. Note
+       * that the vectors need to be of the same size and base on the same map.
+       *
+       * This function is analogous to the the @p swap function of all C++
+       * standard containers. Also, there is a global function
+       * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+       * analogy to standard functions.
+       */
+      void swap (Vector &v);
+
+      /**
+       * Estimate for the memory consumption in bytes.
+       */
+      std::size_t memory_consumption () const;
+
+      /**
+       * Return a reference to the MPI communicator object in use with this
+       * object.
+       */
+      const MPI_Comm &get_mpi_communicator () const;
+      //@}
+
+      /**
+       * Exception
+       */
+      DeclException0 (ExcDifferentParallelPartitioning);
+
+      /**
+       * Exception
+       */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occurred while calling a Trilinos function");
+
+      /**
+       * Exception
+       */
+      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:
+      /**
+       * Trilinos doesn't allow to mix additions to matrix entries and
+       * overwriting them (to make synchronisation of parallel computations
+       * simpler). The way we do it is to, for each access operation, store
+       * whether it is an insertion or an addition. If the previous one was of
+       * different type, then we first have to flush the Trilinos buffers;
+       * otherwise, we can simply go on.  Luckily, Trilinos has an object for
+       * this which does already all the parallel communications in such a case,
+       * so we simply use their model, which stores whether the last operation
+       * was an addition or an insertion.
+       */
+      Epetra_CombineMode last_action;
+
+      /**
+       * A boolean variable to hold information on whether the vector is
+       * compressed or not.
+       */
+      bool compressed;
+
+      /**
+       * Whether this vector has ghost elements. This is true on all processors
+       * even if only one of them has any ghost elements.
+       */
+      bool has_ghosts;
+
+      /**
+       * Pointer to the actual Epetra vector object. This may represent a vector
+       * that is in fact distributed among multiple processors. The object
+       * requires an existing Epetra_Map for storing data when setting it up.
+       */
+      std::shared_ptr<Epetra_FEVector> vector;
+
+      /**
+       * A vector object in Trilinos to be used for collecting the non-local
+       * elements if the vector was constructed with an additional IndexSet
+       * describing ghost elements.
+       */
+      std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
+
+      /**
+       * An IndexSet storing the indices this vector owns exclusively.
+       */
+      IndexSet owned_elements;
+
+      /**
+       * Make the reference class a friend.
+       */
+      friend class internal::VectorReference;
+      friend class Vector;
+      friend class MPI::Vector;
+    };
+
+
+
+
+// ------------------- inline and template functions --------------
+
+
+    /**
+     * Global function @p swap which overloads the default implementation of
+     * the C++ standard library which uses a temporary object. The function
+     * simply exchanges the data of the two vectors.
+     *
+     * @relates TrilinosWrappers::MPI::Vector
+     * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+     */
+    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[0], &values[0]);
+    }
+
+
+
+    inline
+    void
+    Vector::set (const std::vector<size_type>           &indices,
+                 const ::dealii::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
+    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[0], &values[0]);
+    }
+
+
+
+    inline
+    void
+    Vector::add (const std::vector<size_type>           &indices,
+                 const ::dealii::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
+    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);
+    }
 
-    template <typename number>
-    Vector::Vector (const IndexSet               &parallel_partitioner,
-                    const dealii::Vector<number> &v,
-                    const MPI_Comm               &communicator)
+
+
+    inline
+    TrilinosScalar
+    Vector::operator* (const Vector &vec) const
     {
-      *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
-                     v);
-      owned_elements = parallel_partitioner;
+      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 s)
+    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)
     {
-      VectorBase::operator= (s);
+      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;
     }
 
 
-    template <typename Number>
+
+    inline
     Vector &
-    Vector::operator= (const ::dealii::Vector<Number> &v)
+    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 (size() != v.size())
+      // 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
         {
-          vector.reset (new Epetra_FEVector(Epetra_Map
-                                            (static_cast<TrilinosWrappers::types::int_type>(v.size()), 0,
+          // 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
-                                             Epetra_MpiComm(MPI_COMM_SELF)
+
+      const Epetra_MpiComm *mpi_comm
+        = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
+      comm = mpi_comm->Comm();
+
 #else
-                                             Epetra_SerialComm()
+
+      comm = MPI_COMM_SELF;
+
 #endif
-                                            )));
-        }
 
-      const Epetra_Map &map = vector_partitioner();
-      const int size = map.NumMyElements();
+      return comm;
+    }
+
+    template <typename number>
+    Vector::Vector (const IndexSet               &parallel_partitioner,
+                    const dealii::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));
 
-      // Need to copy out values, since the deal.II might not use doubles, so
-      // that a direct access is not possible.
-      for (int i=0; i<size; ++i)
-        (*vector)[0][i] = v(gid(map,i));
+      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 MPI */
 } /* end of namespace TrilinosWrappers */
 
 /*@}*/
@@ -630,51 +2193,10 @@ namespace internal
       }
     };
 
-    /**
-     * A helper class used internally in linear_operator.h. Specialization for
-     * TrilinosWrappers::VectorBase.
-     */
-    template<>
-    class ReinitHelper<TrilinosWrappers::VectorBase>
-    {
-    public:
-      template <typename Matrix>
-      static
-      void reinit_range_vector (const Matrix &matrix,
-                                TrilinosWrappers::VectorBase &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::VectorBase &v,
-                                bool omit_zeroing_entries)
-      {
-        v.reinit(matrix.locally_owned_domain_indices(),
-                 matrix.get_mpi_communicator(),
-                 omit_zeroing_entries);
-      }
-    };
-
   } /* namespace LinearOperator */
 } /* namespace internal */
 
 
-/**
- * Declare dealii::TrilinosWrappers::VectorBase as serial vector.
- *
- * @author Uwe Koecher, 2017
- */
-template <>
-struct is_serial_vector< TrilinosWrappers::VectorBase > : std::true_type
-{
-};
-
 
 /**
  * Declare dealii::TrilinosWrappers::MPI::Vector as distributed vector.
@@ -693,5 +2215,4 @@ DEAL_II_NAMESPACE_CLOSE
 
 /*----------------------------   trilinos_vector.h     ---------------------------*/
 
-#endif
-/*----------------------------   trilinos_vector.h     ---------------------------*/
+#endif // dealii__trilinos_vector_h
index 523db9b0d9854707ad67fea832feb9df55ea3fa3..10b5f44f6bf0ccea6456d528f77deb0d9c33838d 100644 (file)
 #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/subscriptor.h>
-#  include <deal.II/lac/exceptions.h>
-#  include <deal.II/lac/vector.h>
-#  include <deal.II/base/mpi.h>
-
-#  include <vector>
-#  include <utility>
-#  include <memory>
-
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-#  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
+#include <deal.II/lac/trilinos_vector.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-// forward declaration
-template <typename number> class Vector;
-
-
-/**
- * @addtogroup TrilinosWrappers
- * @{
- */
-
-/**
- * A namespace in which wrapper classes for Trilinos objects reside.
- *
- * @ingroup TrilinosWrappers
- */
+#ifdef DEAL_II_WITH_TRILINOS
 namespace TrilinosWrappers
 {
-  // forward declaration
-  class VectorBase;
-  namespace MPI
-  {
-    class Vector;
-  }
-
-  /**
-   * @cond internal
-   */
-
-  /**
-   * A namespace for internal implementation details of the TrilinosWrapper
-   * members.
-   *
-   * @ingroup TrilinosWrappers
-   */
-  namespace internal
-  {
-    /**
-     * Declare type for container size.
-     */
-    typedef dealii::types::global_dof_index size_type;
-
-    /**
-     * This class implements a wrapper for accessing the Trilinos vector in
-     * the same way as we access deal.II objects: it is initialized with a
-     * vector and an element within it, and has a conversion operator to
-     * extract the scalar value of this element. It also has a variety of
-     * assignment operator for writing to this one element.
-     *
-     * @ingroup TrilinosWrappers
-     */
-    class VectorReference
-    {
-    private:
-      /**
-       * Constructor. It is made private so as to only allow the actual vector
-       * class to create it.
-       */
-      VectorReference (VectorBase     &vector,
-                       const size_type  index);
-
-    public:
-
-      /**
-       * This looks like a copy operator, but does something different than
-       * usual. In particular, it does not copy the member variables of this
-       * reference. Rather, it handles the situation where we have two vectors
-       * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
-       * both left and right hand side of the assignment have data type
-       * VectorReference, but what we really mean is to assign the vector
-       * elements represented by the two references. This operator implements
-       * this operation. Note also that this allows us to make the assignment
-       * operator const.
-       */
-      const VectorReference &
-      operator = (const VectorReference &r) const;
-
-      /**
-       * Same as above but for non-const reference objects.
-       */
-      VectorReference &
-      operator = (const VectorReference &r);
-
-      /**
-       * Set the referenced element of the vector to <tt>s</tt>.
-       */
-      const VectorReference &
-      operator = (const TrilinosScalar &s) const;
-
-      /**
-       * Add <tt>s</tt> to the referenced element of the vector->
-       */
-      const VectorReference &
-      operator += (const TrilinosScalar &s) const;
-
-      /**
-       * Subtract <tt>s</tt> from the referenced element of the vector->
-       */
-      const VectorReference &
-      operator -= (const TrilinosScalar &s) const;
-
-      /**
-       * Multiply the referenced element of the vector by <tt>s</tt>.
-       */
-      const VectorReference &
-      operator *= (const TrilinosScalar &s) const;
-
-      /**
-       * Divide the referenced element of the vector by <tt>s</tt>.
-       */
-      const VectorReference &
-      operator /= (const TrilinosScalar &s) const;
-
-      /**
-       * Convert the reference to an actual value, i.e. return the value of
-       * the referenced element of the vector.
-       */
-      operator TrilinosScalar () const;
-
-      /**
-       * Exception
-       */
-      DeclException1 (ExcTrilinosError,
-                      int,
-                      << "An error with error number " << arg1
-                      << " occurred while calling a Trilinos function");
-
-    private:
-      /**
-       * Point to the vector we are referencing.
-       */
-      VectorBase   &vector;
-
-      /**
-       * Index of the referenced element of the vector.
-       */
-      const size_type  index;
-
-      /**
-       * Make the vector class a friend, so that it can create objects of the
-       * present type.
-       */
-      friend class ::dealii::TrilinosWrappers::VectorBase;
-    };
-  }
-  /**
-   * @endcond
-   */
-
-
-  /**
-   * Base class for the the distributed memory vector MPI::Vector that is
-   * designed for use in either serial implementations or as a localized copy
-   * on each processor.  The implementation of this class is based on the
-   * Trilinos vector class Epetra_FEVector, the (parallel) partitioning of which
-   * is governed by an Epetra_Map. This means that the vector type is generic
-   * and can be done in this base class, while the definition of the partition
-   * map (and hence, the constructor and reinit function) will have to be done
-   * in the derived classes. The Epetra_FEVector is precisely the kind of vector
-   * we deal with all the time - we probably get it from some assembly
-   * process, where also entries not locally owned might need to written and
-   * hence need to be forwarded to the owner. The only requirement for this
-   * class to work is that Trilinos is installed with the same compiler as is
-   * used for compilation of deal.II.
-   *
-   * The interface of this class is modeled after the existing Vector class in
-   * deal.II. It has almost the same member functions, and is often
-   * exchangeable. However, since Trilinos only supports a single scalar type
-   * (double), it is not templated, and only works with that type.
-   *
-   * Note that Trilinos only guarantees that operations do what you expect if
-   * the function @p GlobalAssemble has been called after vector assembly in
-   * order to distribute the data. Therefore, you need to call
-   * Vector::compress() before you actually use the vectors.
-   *
-   * @ingroup TrilinosWrappers
-   * @ingroup Vectors
-   * @author Martin Kronbichler, 2008
-   */
-  class VectorBase : public Subscriptor
-  {
-  public:
-    /**
-     * Declare some of the standard types used in all containers. These types
-     * parallel those in the <tt>C</tt> standard libraries
-     * <tt>vector<...></tt> class.
-     */
-    typedef TrilinosScalar                  value_type;
-    typedef TrilinosScalar                  real_type;
-    typedef dealii::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;
-
-    /**
-     * @name 1: Basic Object-handling
-     */
-    //@{
-
-    /**
-     * Default constructor that generates an empty (zero size) vector. The
-     * function <tt>reinit()</tt> will have to give the vector the correct
-     * size and distribution among processes in case of an MPI run.
-     */
-    VectorBase ();
-
-    /**
-     * Copy constructor. Sets the dimension to that of the given vector, and
-     * copies all the elements.
-     */
-    VectorBase (const VectorBase &v);
-
-    /**
-     * Destructor
-     */
-    virtual ~VectorBase ();
-
-    /**
-     * Release all memory and return to a state just like after having called
-     * the default constructor.
-     */
-    void clear ();
-
-    /**
-     * Reinit functionality, sets the dimension and possibly the parallel
-     * partitioning (Epetra_Map) of the calling vector to the settings of the
-     * input vector.
-     */
-    void reinit (const VectorBase &v,
-                 const bool        omit_zeroing_entries = false);
-
-    /**
-     * Compress the underlying representation of the Trilinos object, i.e.
-     * flush the buffers of the vector object if it has any. This function is
-     * necessary after writing into a vector element-by-element and before
-     * anything else can be done on it.
-     *
-     * The (defaulted) argument can be used to specify the compress mode
-     * (<code>Add</code> or <code>Insert</code>) in case the vector has not
-     * been written to since the last time this function was called. The
-     * argument is ignored if the vector has been added or written to since
-     * the last time compress() was called.
-     *
-     * See
-     * @ref GlossCompress "Compressing distributed objects"
-     * for more information.
-     */
-    void compress (::dealii::VectorOperation::values operation);
-
-    /**
-     * Set all components of the vector to the given number @p s. Simply pass
-     * this down to the Trilinos Epetra object, but we still need to declare
-     * this function to make the example given in the discussion about making
-     * the constructor explicit work.
-     *
-     * Since the semantics of assigning a scalar to a vector are not
-     * immediately clear, this operator should really only be used if you want
-     * to set the entire vector to zero. This allows the intuitive notation
-     * <tt>v=0</tt>. Assigning other values is deprecated and may be
-     * disallowed in the future.
-     */
-    VectorBase &
-    operator = (const TrilinosScalar s);
-
-    /**
-     * Copy function. This function takes a VectorBase vector and copies all
-     * the elements. The target vector will have the same parallel
-     * distribution as the calling vector.
-     */
-    VectorBase &
-    operator = (const VectorBase &v);
-
-    /**
-     * Another copy function. This one takes a deal.II vector and copies it
-     * into a TrilinosWrapper vector. Note that since we do not provide any
-     * Epetra_map that tells about the partitioning of the vector among the
-     * MPI processes, the size of the TrilinosWrapper vector has to be the
-     * same as the size of the input vector. In order to change the map, use
-     * the reinit(const Epetra_Map &input_map) function.
-     */
-    template <typename Number>
-    VectorBase &
-    operator = (const ::dealii::Vector<Number> &v);
-
-    /**
-     * Test for equality. This function assumes that the present vector and
-     * the one to compare with have the same size already, since comparing
-     * vectors of different sizes makes not much sense anyway.
-     */
-    bool operator == (const VectorBase &v) const;
-
-    /**
-     * Test for inequality. This function assumes that the present vector and
-     * the one to compare with have the same size already, since comparing
-     * vectors of different sizes makes not much sense anyway.
-     */
-    bool operator != (const VectorBase &v) const;
-
-    /**
-     * Return the global dimension of the vector.
-     */
-    size_type size () const;
-
-    /**
-     * Return the local dimension of the vector, i.e. the number of elements
-     * stored on the present MPI process. For sequential vectors, this number
-     * is the same as size(), but for parallel vectors it may be smaller.
-     *
-     * To figure out which elements exactly are stored locally, use
-     * local_range().
-     *
-     * If the vector contains ghost elements, they are included in this
-     * number.
-     */
-    size_type local_size () const;
-
-    /**
-     * Return a pair of indices indicating which elements of this vector are
-     * stored locally. The first number is the index of the first element
-     * stored, the second the index of the one past the last one that is
-     * stored locally. If this is a sequential vector, then the result will be
-     * the pair <code>(0,N)</code>, otherwise it will be a pair
-     * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
-     * <code>i</code> is the first element of the vector stored on this
-     * processor, corresponding to the half open interval $[i,i+n)$
-     *
-     * @note The description above is true most of the time, but not always.
-     * In particular, Trilinos vectors need not store contiguous ranges of
-     * elements such as $[i,i+n)$. Rather, it can store vectors where the
-     * elements are distributed in an arbitrary way across all processors and
-     * each processor simply stores a particular subset, not necessarily
-     * contiguous. In this case, this function clearly makes no sense since it
-     * could, at best, return a range that includes all elements that are
-     * stored locally. Thus, the function only succeeds if the locally stored
-     * range is indeed contiguous. It will trigger an assertion if the local
-     * portion of the vector is not contiguous.
-     */
-    std::pair<size_type, size_type> local_range () const;
-
-    /**
-     * Return whether @p index is in the local range or not, see also
-     * local_range().
-     *
-     * @note The same limitation for the applicability of this function
-     * applies as listed in the documentation of local_range().
-     */
-    bool in_local_range (const size_type index) const;
-
-    /**
-     * Return an index set that describes which elements of this vector are
-     * owned by the current processor. Note that this index set does not
-     * include elements this vector may store locally as ghost elements but
-     * that are in fact owned by another processor. As a consequence, the
-     * index sets returned on different processors if this is a distributed
-     * vector will form disjoint sets that add up to the complete index set.
-     * Obviously, if a vector is created on only one processor, then the
-     * result would satisfy
-     * @code
-     *   vec.locally_owned_elements() == complete_index_set (vec.size())
-     * @endcode
-     */
-    IndexSet locally_owned_elements () const;
-
-    /**
-     * Return if the vector contains ghost elements. This answer is true if
-     * there are ghost elements on at least one process.
-     *
-     * @see
-     * @ref GlossGhostedVector "vectors with ghost elements"
-     */
-    bool has_ghost_elements() const;
-
-    /**
-     * This function only exists for compatibility with the @p
-     * LinearAlgebra::distributed::Vector class and does nothing: this class
-     * implements ghost value updates in a different way that is a better fit
-     * with the underlying Trilinos vector object.
-     */
-    void update_ghost_values () const;
-
-    /**
-     * Return the scalar (inner) product of two vectors. The vectors must have
-     * the same size.
-     */
-    TrilinosScalar operator * (const VectorBase &vec) const;
-
-    /**
-     * Return square of the $l_2$-norm.
-     */
-    real_type norm_sqr () const;
-
-    /**
-     * Mean value of the elements of this vector.
-     */
-    TrilinosScalar mean_value () const;
-
-    /**
-     * Compute the minimal value of the elements of this vector.
-     */
-    TrilinosScalar min () const;
-
-    /**
-     * Compute the maximal value of the elements of this vector.
-     */
-    TrilinosScalar max () const;
-
-    /**
-     * $l_1$-norm of the vector.  The sum of the absolute values.
-     */
-    real_type l1_norm () const;
-
-    /**
-     * $l_2$-norm of the vector.  The square root of the sum of the squares of
-     * the elements.
-     */
-    real_type l2_norm () const;
-
-    /**
-     * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
-     * <i>p</i>th powers of the absolute values of the elements.
-     */
-    real_type lp_norm (const TrilinosScalar p) const;
-
-    /**
-     * Maximum absolute value of the elements.
-     */
-    real_type linfty_norm () const;
-
-    /**
-     * Performs a combined operation of a vector addition and a subsequent
-     * inner product, returning the value of the inner product. In other
-     * words, the result of this function is the same as if the user called
-     * @code
-     * this->add(a, V);
-     * return_value = *this * W;
-     * @endcode
-     *
-     * The reason this function exists is for compatibility with deal.II's own
-     * vector classes which can implement this functionality with less memory
-     * transfer. However, for Trilinos vectors such a combined operation is
-     * not natively supported and thus the cost is completely equivalent as
-     * calling the two methods separately.
-     */
-    TrilinosScalar add_and_dot (const TrilinosScalar a,
-                                const VectorBase    &V,
-                                const VectorBase    &W);
-
-    /**
-     * Return whether the vector contains only elements with value zero. This
-     * is a collective operation. This function is expensive, because
-     * potentially all elements have to be checked.
-     */
-    bool all_zero () const;
-
-    /**
-     * Return @p true if the vector has no negative entries, i.e. all entries
-     * are zero or positive. This function is used, for example, to check
-     * whether refinement indicators are really all positive (or zero).
-     */
-    bool is_non_negative () const;
-    //@}
-
-
-    /**
-     * @name 2: Data-Access
-     */
-    //@{
-
-    /**
-     * Provide access to a given element, both read and write.
-     *
-     * When using a vector distributed with MPI, this operation only makes
-     * sense for elements that are actually present on the calling processor.
-     * Otherwise, an exception is thrown. This is different from the
-     * <code>el()</code> function below that always succeeds (but returns zero
-     * on non-local elements).
-     */
-    reference
-    operator () (const size_type index);
-
-    /**
-     * Provide read-only access to an element.
-     *
-     * When using a vector distributed with MPI, this operation only makes
-     * sense for elements that are actually present on the calling processor.
-     * Otherwise, an exception is thrown. This is different from the
-     * <code>el()</code> function below that always succeeds (but returns zero
-     * on non-local elements).
-     */
-    TrilinosScalar
-    operator () (const size_type index) const;
-
-    /**
-     * Provide access to a given element, both read and write.
-     *
-     * Exactly the same as operator().
-     */
-    reference
-    operator [] (const size_type index);
-
-    /**
-     * Provide read-only access to an element.
-     *
-     * Exactly the same as operator().
-     */
-    TrilinosScalar
-    operator [] (const size_type index) const;
-
-    /**
-     * A collective get operation: instead of getting individual elements of a
-     * vector, this function allows to get a whole set of elements at once.
-     * The indices of the elements to be read are stated in the first
-     * argument, the corresponding values are returned in the second.
-     */
-    void extract_subvector_to (const std::vector<size_type> &indices,
-                               std::vector<TrilinosScalar> &values) const;
-
-    /**
-     * Just as the above, but with pointers.  Useful in minimizing copying of
-     * data around.
-     */
-    template <typename ForwardIterator, typename OutputIterator>
-    void extract_subvector_to (ForwardIterator          indices_begin,
-                               const ForwardIterator    indices_end,
-                               OutputIterator           values_begin) const;
-
-    /**
-     * Make the Vector class a bit like the <tt>vector<></tt> class of the C++
-     * standard library by returning iterators to the start and end of the
-     * locally owned elements of this vector. The ordering of local elements
-     * corresponds to the one given by the global indices in case the vector
-     * is constructed from an IndexSet or other methods in deal.II (note that
-     * an Epetra_Map can contain elements in arbitrary orders, though).
-     *
-     * It holds that end() - begin() == local_size().
-     */
-    iterator begin ();
-
-    /**
-     * Return constant iterator to the start of the locally owned elements of
-     * the vector.
-     */
-    const_iterator begin () const;
-
-    /**
-     * Return an iterator pointing to the element past the end of the array of
-     * locally owned entries.
-     */
-    iterator end ();
-
-    /**
-     * Return a constant iterator pointing to the element past the end of the
-     * array of the locally owned entries.
-     */
-    const_iterator end () const;
-
-    //@}
-
-
-    /**
-     * @name 3: Modification of vectors
-     */
-    //@{
-
-    /**
-     * A collective set operation: instead of setting individual elements of a
-     * vector, this function allows to set a whole set of elements at once.
-     * The indices of the elements to be set are stated in the first argument,
-     * the corresponding values in the second.
-     */
-    void set (const std::vector<size_type>    &indices,
-              const std::vector<TrilinosScalar>  &values);
-
-    /**
-     * This is a second collective set operation. As a difference, this
-     * function takes a deal.II vector of values.
-     */
-    void set (const std::vector<size_type>        &indices,
-              const ::dealii::Vector<TrilinosScalar> &values);
-
-    /**
-     * This collective set operation is of lower level and can handle anything
-     * else &mdash; the only thing you have to provide is an address where all
-     * the indices are stored and the number of elements to be set.
-     */
-    void set (const size_type       n_elements,
-              const size_type      *indices,
-              const TrilinosScalar *values);
-
-    /**
-     * A collective add operation: This function adds a whole set of values
-     * stored in @p values to the vector components specified by @p indices.
-     */
-    void add (const std::vector<size_type>      &indices,
-              const std::vector<TrilinosScalar> &values);
-
-    /**
-     * This is a second collective add operation. As a difference, this
-     * function takes a deal.II vector of values.
-     */
-    void add (const std::vector<size_type>           &indices,
-              const ::dealii::Vector<TrilinosScalar> &values);
-
-    /**
-     * Take an address where <tt>n_elements</tt> are stored contiguously and
-     * add them into the vector. Handles all cases which are not covered by
-     * the other two <tt>add()</tt> functions above.
-     */
-    void add (const size_type       n_elements,
-              const size_type      *indices,
-              const TrilinosScalar *values);
-
-    /**
-     * Multiply the entire vector by a fixed factor.
-     */
-    VectorBase &operator *= (const TrilinosScalar factor);
-
-    /**
-     * Divide the entire vector by a fixed factor.
-     */
-    VectorBase &operator /= (const TrilinosScalar factor);
-
-    /**
-     * Add the given vector to the present one.
-     */
-    VectorBase &operator += (const VectorBase &V);
-
-    /**
-     * Subtract the given vector from the present one.
-     */
-    VectorBase &operator -= (const VectorBase &V);
-
-    /**
-     * Addition of @p s to all components. Note that @p s is a scalar and not
-     * a vector.
-     */
-    void add (const TrilinosScalar s);
-
-    /**
-     * Simple vector addition, equal to the <tt>operator +=</tt>.
-     *
-     * Though, if the second argument <tt>allow_different_maps</tt> is set,
-     * then it is possible to add data from a vector that uses a different
-     * map, i.e., a vector whose elements are split across processors
-     * differently. This may include vectors with ghost elements, for example.
-     * In general, however, adding vectors with a different element-to-
-     * processor map requires communicating data among processors and,
-     * consequently, is a slower operation than when using vectors using the
-     * same map.
-     */
-    void add (const VectorBase &V,
-              const bool        allow_different_maps = false);
-
-    /**
-     * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
-     */
-    void add (const TrilinosScalar  a,
-              const VectorBase     &V);
-
-    /**
-     * Multiple addition of scaled vectors, i.e. <tt>*this += a*V + b*W</tt>.
-     */
-    void add (const TrilinosScalar  a,
-              const VectorBase     &V,
-              const TrilinosScalar  b,
-              const VectorBase     &W);
-
-    /**
-     * Scaling and simple vector addition, i.e.  <tt>*this = s*(*this) +
-     * V</tt>.
-     */
-    void sadd (const TrilinosScalar  s,
-               const VectorBase     &V);
-
-    /**
-     * Scaling and simple addition, i.e.  <tt>*this = s*(*this) + a*V</tt>.
-     */
-    void sadd (const TrilinosScalar  s,
-               const TrilinosScalar  a,
-               const VectorBase     &V);
-
-    /**
-     * Scale each element of this vector by the corresponding element in the
-     * argument. This function is mostly meant to simulate multiplication (and
-     * immediate re-assignment) by a diagonal scaling matrix.
-     */
-    void scale (const VectorBase &scaling_factors);
-
-    /**
-     * Assignment <tt>*this = a*V</tt>.
-     */
-    void equ (const TrilinosScalar  a,
-              const VectorBase     &V);
-    //@}
-
-
-    /**
-     * @name 4: Mixed stuff
-     */
-    //@{
-
-    /**
-     * Return a const reference to the underlying Trilinos Epetra_MultiVector
-     * class.
-     */
-    const Epetra_MultiVector &trilinos_vector () const;
-
-    /**
-     * Return a (modifyable) reference to the underlying Trilinos
-     * Epetra_FEVector class.
-     */
-    Epetra_FEVector &trilinos_vector ();
-
-    /**
-     * Return a const reference to the underlying Trilinos Epetra_Map that
-     * sets the parallel partitioning of the vector.
-     */
-    const Epetra_Map &vector_partitioner () const;
-
-    /**
-     * Print to a stream. @p precision denotes the desired precision with
-     * which values shall be printed, @p scientific whether scientific
-     * notation shall be used. If @p across is @p true then the vector is
-     * printed in a line, while if @p false then the elements are printed on a
-     * separate line each.
-     */
-    void print (std::ostream       &out,
-                const unsigned int  precision  = 3,
-                const bool          scientific = true,
-                const bool          across     = true) const;
-
-    /**
-     * Swap the contents of this vector and the other vector @p v. One could
-     * do this operation with a temporary variable and copying over the data
-     * elements, but this function is significantly more efficient since it
-     * only swaps the pointers to the data of the two vectors and therefore
-     * does not need to allocate temporary storage and move data around. Note
-     * that the vectors need to be of the same size and base on the same map.
-     *
-     * This function is analogous to the the @p swap function of all C++
-     * standard containers. Also, there is a global function
-     * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
-     * analogy to standard functions.
-     */
-    void swap (VectorBase &v);
-
-    /**
-     * Estimate for the memory consumption in bytes.
-     */
-    std::size_t memory_consumption () const;
-
-    /**
-     * Return a reference to the MPI communicator object in use with this
-     * object.
-     */
-    const MPI_Comm &get_mpi_communicator () const;
-    //@}
-
-    /**
-     * Exception
-     */
-    DeclException0 (ExcDifferentParallelPartitioning);
-
-    /**
-     * Exception
-     */
-    DeclException1 (ExcTrilinosError,
-                    int,
-                    << "An error with error number " << arg1
-                    << " occurred while calling a Trilinos function");
-
-    /**
-     * Exception
-     */
-    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:
-    /**
-     * Trilinos doesn't allow to mix additions to matrix entries and
-     * overwriting them (to make synchronisation of parallel computations
-     * simpler). The way we do it is to, for each access operation, store
-     * whether it is an insertion or an addition. If the previous one was of
-     * different type, then we first have to flush the Trilinos buffers;
-     * otherwise, we can simply go on.  Luckily, Trilinos has an object for
-     * this which does already all the parallel communications in such a case,
-     * so we simply use their model, which stores whether the last operation
-     * was an addition or an insertion.
-     */
-    Epetra_CombineMode last_action;
-
-    /**
-     * A boolean variable to hold information on whether the vector is
-     * compressed or not.
-     */
-    bool compressed;
-
-    /**
-     * Whether this vector has ghost elements. This is true on all processors
-     * even if only one of them has any ghost elements.
-     */
-    bool has_ghosts;
-
-    /**
-     * Pointer to the actual Epetra vector object. This may represent a vector
-     * that is in fact distributed among multiple processors. The object
-     * requires an existing Epetra_Map for storing data when setting it up.
-     */
-    std::shared_ptr<Epetra_FEVector> vector;
-
-    /**
-     * A vector object in Trilinos to be used for collecting the non-local
-     * elements if the vector was constructed with an additional IndexSet
-     * describing ghost elements.
-     */
-    std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
-
-    /**
-     * An IndexSet storing the indices this vector owns exclusively.
-     */
-    IndexSet owned_elements;
-
-    /**
-     * Make the reference class a friend.
-     */
-    friend class internal::VectorReference;
-    friend class Vector;
-    friend class MPI::Vector;
-  };
-
-
-
-
-// ------------------- inline and template functions --------------
-
-  /**
-   * Global function swap which overloads the default implementation of the C
-   * standard library which uses a temporary object. The function simply
-   * exchanges the data of the two vectors.
-   *
-   * @relates TrilinosWrappers::VectorBase
-   * @author Martin Kronbichler, Wolfgang Bangerth, 2008
-   */
-  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
-    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::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
-  {
-    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
-  VectorBase::has_ghost_elements() const
-  {
-    return has_ghosts;
-  }
-
-
-
-  inline
-  void
-  VectorBase::update_ghost_values () const
-  {}
-
-
-
-  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() != nullptr,
-            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() != nullptr)
-      nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
-
-    owned_elements = v.owned_elements;
-  }
-
-
-
-  inline
-  VectorBase &
-  VectorBase::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;
-  }
-
-
-
-  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 ::dealii::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)
-      {
-        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
-  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 ::dealii::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)
-          {
-            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
-  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::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::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
-  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
-
+  // All functionality from the former base class VectorBase was moved to
+  // MPI::Vector as this ended up to be the only derived class.
+  // Hence, using the class VectorBase is deprecated and we only provide
+  // an alias for backward compatibility.
+  using VectorBase = MPI::Vector;
 }
-
-/*@}*/
+#endif
 
 DEAL_II_NAMESPACE_CLOSE
 
-#endif // DEAL_II_WITH_TRILINOS
-
-/*----------------------------   trilinos_vector_base.h     ---------------------------*/
-
-#endif
-/*----------------------------   trilinos_vector_base.h     ---------------------------*/
+#endif // dealii__trilinos_vector_base_h
index 5597406a10b19a082971d5fff79b2324ffee9341..3c6dba8d11d0956059e9f765df9dc03081f4bf63 100644 (file)
@@ -121,13 +121,12 @@ IF(DEAL_II_WITH_TRILINOS)
     trilinos_solver.cc
     trilinos_sparse_matrix.cc
     trilinos_sparsity_pattern.cc
-    trilinos_vector_base.cc
     trilinos_vector.cc
   )
   SET(_inst
     ${_inst}
     trilinos_sparse_matrix.inst.in
-    trilinos_vector_base.inst.in
+    trilinos_vector.inst.in
   )
 ENDIF()
 
index 446456765c86bbeeb7ee98b53747995869f9c8c1..acab363c1d9c989c546c35c25a7cea20260939de 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-#include <deal.II/lac/trilinos_index_access.h>
 #include <deal.II/lac/trilinos_vector.h>
 
 #ifdef DEAL_II_WITH_TRILINOS
 
 #  include <deal.II/base/mpi.h>
+#  include <deal.II/lac/trilinos_sparse_matrix.h>
 #  include <deal.II/lac/trilinos_parallel_block_vector.h>
 #  include <deal.II/lac/trilinos_index_access.h>
-#  include <deal.II/lac/trilinos_sparse_matrix.h>
 
 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
 #  include <Epetra_Import.h>
+#  include <Epetra_Export.h>
 #  include <Epetra_Vector.h>
+
+#  include <boost/io/ios_state.hpp>
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
 #  include <cmath>
@@ -35,13 +37,44 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace TrilinosWrappers
 {
-  namespace MPI
+  namespace internal
   {
-    Vector::Vector ()
+    VectorReference::operator TrilinosScalar () const
     {
-      last_action = Zero;
-      vector.reset(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self())));
+      Assert (index < vector.size(),
+              ExcIndexRange (index, 0, vector.size()));
+
+      // Trilinos allows for vectors to be referenced by the [] or ()
+      // operators but only () checks index bounds. We check these bounds by
+      // ourselves, so we can use []. Note that we can only get local values.
+
+      const TrilinosWrappers::types::int_type local_index =
+        vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
+      Assert (local_index >= 0,
+              MPI::Vector::ExcAccessToNonLocalElement (index, vector.local_size(),
+                                                       vector.vector->Map().MinMyGID(),
+                                                       vector.vector->Map().MaxMyGID()));
+
+
+      return (*(vector.vector))[0][local_index];
     }
+  }
+
+  namespace MPI
+  {
+    Vector::Vector ()
+      :
+      last_action (Zero),
+      compressed  (true),
+      has_ghosts  (false),
+#ifdef DEAL_II_WITH_MPI
+      vector(new Epetra_FEVector(
+               Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_SELF))))
+#else
+      vector(new Epetra_FEVector(
+               Epetra_Map(0,0,Epetra_SerialComm())))
+#endif
+    {}
 
 
 
@@ -55,13 +88,13 @@ namespace TrilinosWrappers
 
     Vector::Vector (const Vector &v)
       :
-      VectorBase()
-    {
-      last_action = Zero;
-      vector.reset (new Epetra_FEVector(*v.vector));
-      has_ghosts = v.has_ghosts;
-      owned_elements = v.owned_elements;
-    }
+      Subscriptor(),
+      last_action (Zero),
+      compressed (true),
+      has_ghosts  (v.has_ghosts),
+      vector(new Epetra_FEVector(*v.vector)),
+      owned_elements(v.owned_elements)
+    {}
 
 
 
@@ -78,29 +111,32 @@ namespace TrilinosWrappers
 
 
     Vector::Vector (const IndexSet   &parallel_partitioner,
-                    const VectorBase &v,
+                    const Vector     &v,
                     const MPI_Comm   &communicator)
       :
-      VectorBase()
+      last_action (Zero),
+      compressed  (true),
+      has_ghosts  (false)
     {
       AssertThrow (parallel_partitioner.size() ==
                    static_cast<size_type>(TrilinosWrappers::n_global_elements(v.vector->Map())),
                    ExcDimensionMismatch (parallel_partitioner.size(),
                                          TrilinosWrappers::n_global_elements(v.vector->Map())));
 
-      last_action = Zero;
-
       vector.reset (new Epetra_FEVector
                     (parallel_partitioner.make_trilinos_map(communicator,
                                                             true)));
       reinit (v, false, true);
     }
 
+
+
     Vector::Vector (const IndexSet &local,
                     const IndexSet &ghost,
                     const MPI_Comm &communicator)
       :
-      VectorBase()
+      last_action (Zero),
+      compressed  (true)
     {
       reinit(local, ghost, communicator, false);
     }
@@ -111,6 +147,24 @@ namespace TrilinosWrappers
 
 
 
+    void
+    Vector::clear ()
+    {
+      // When we clear the vector, reset the pointer and generate an empty
+      // vector.
+#ifdef DEAL_II_WITH_MPI
+      Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
+#else
+      Epetra_Map map (0, 0, Epetra_SerialComm());
+#endif
+
+      has_ghosts = false;
+      vector.reset (new Epetra_FEVector(map));
+      last_action = Zero;
+    }
+
+
+
     void
     Vector::reinit (const IndexSet &parallel_partitioner,
                     const MPI_Comm &communicator,
@@ -150,7 +204,7 @@ namespace TrilinosWrappers
 
 
     void
-    Vector::reinit (const VectorBase &v,
+    Vector::reinit (const Vector &v,
                     const bool        omit_zeroing_entries,
                     const bool        allow_different_maps)
     {
@@ -367,6 +421,9 @@ namespace TrilinosWrappers
     Vector &
     Vector::operator = (const Vector &v)
     {
+      Assert (vector.get() != nullptr,
+              ExcMessage("Vector is not constructed properly."));
+
       // check equality for MPI communicators to avoid accessing a possibly
       // invalid MPI_Comm object
 #ifdef DEAL_II_WITH_MPI
@@ -443,19 +500,21 @@ namespace TrilinosWrappers
 
 
 
+    template <typename number>
     Vector &
-    Vector::operator = (const TrilinosWrappers::VectorBase &v)
+    Vector::operator = (const ::dealii::Vector<number> &v)
     {
-      nonlocal_vector.reset();
-
-      Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
-
-      Epetra_Import data_exchange (vector->Map(), v.vector->Map());
-      const int ierr = vector->Import(*v.vector, data_exchange, Insert);
-
-      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      Assert (size() == v.size(),
+              ExcDimensionMismatch(size(), v.size()));
 
-      last_action = Insert;
+      // this is probably not very efficient but works. in particular, we could
+      // do better if we know that number==TrilinosScalar because then we could
+      // elide the copying of elements
+      //
+      // let's hope this isn't a particularly frequent operation
+      std::pair<size_type, size_type> local_range = this->local_range ();
+      for (size_type i=local_range.first; i<local_range.second; ++i)
+        (*vector)[0][i-local_range.first] = v(i);
 
       return *this;
     }
@@ -488,9 +547,298 @@ namespace TrilinosWrappers
       last_action = Insert;
     }
 
+
+
+    void
+    Vector::compress (::dealii::VectorOperation::values given_last_action)
+    {
+      //Select which mode to send to Trilinos. Note that we use last_action if
+      //available and ignore what the user tells us to detect wrongly mixed
+      //operations. Typically given_last_action is only used on machines that do
+      //not execute an operation (because they have no own cells for example).
+      Epetra_CombineMode mode = last_action;
+      if (last_action == Zero)
+        {
+          if (given_last_action==::dealii::VectorOperation::add)
+            mode = Add;
+          else if (given_last_action==::dealii::VectorOperation::insert)
+            mode = Insert;
+        }
+
+#ifdef DEBUG
+#  ifdef DEAL_II_WITH_MPI
+      // check that every process has decided to use the same mode. This will
+      // otherwise result in undefined behaviour in the call to
+      // GlobalAssemble().
+      double double_mode = mode;
+      const Epetra_MpiComm *comm_ptr
+        = dynamic_cast<const Epetra_MpiComm *>(&(vector_partitioner().Comm()));
+      Assert (comm_ptr != nullptr, ExcInternalError());
+      Utilities::MPI::MinMaxAvg result
+        = Utilities::MPI::min_max_avg (double_mode, comm_ptr->GetMpiComm());
+      Assert(result.max-result.min<1e-5,
+             ExcMessage ("Not all processors agree whether the last operation on "
+                         "this vector was an addition or a set operation. This will "
+                         "prevent the compress() operation from succeeding."));
+
+#  endif
+#endif
+
+      // Now pass over the information about what we did last to the vector.
+      int ierr = 0;
+      if (nonlocal_vector.get() == nullptr || mode != Add)
+        ierr = vector->GlobalAssemble(mode);
+      else
+        {
+          Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
+          ierr = vector->Export(*nonlocal_vector, exporter, mode);
+          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+          ierr = nonlocal_vector->PutScalar(0.);
+        }
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      last_action = Zero;
+
+      compressed = true;
+    }
+
+
+
+    TrilinosScalar
+    Vector::operator () (const size_type index) const
+    {
+      // Extract local indices in the vector.
+      TrilinosWrappers::types::int_type trilinos_i =
+        vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
+      TrilinosScalar value = 0.;
+
+      // If the element is not present on the current processor, we can't
+      // continue. This is the main difference to the el() function.
+      if (trilinos_i == -1)
+        {
+          Assert (false, ExcAccessToNonLocalElement(index, local_size(),
+                                                    vector->Map().MinMyGID(),
+                                                    vector->Map().MaxMyGID()));
+        }
+      else
+        value = (*vector)[0][trilinos_i];
+
+      return value;
+    }
+
+
+
+    void
+    Vector::add (const Vector &v,
+                 const bool        allow_different_maps)
+    {
+      if (allow_different_maps == false)
+        *this += v;
+      else
+        {
+          Assert (!has_ghost_elements(), ExcGhostsPresent());
+          AssertThrow (size() == v.size(),
+                       ExcDimensionMismatch (size(), v.size()));
+
+#if DEAL_II_TRILINOS_VERSION_GTE(11,11,0)
+          Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+          int ierr = vector->Import(*v.vector, data_exchange, Epetra_AddLocalAlso);
+          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+          last_action = Add;
+#else
+          // In versions older than 11.11 the Import function is broken for adding
+          // Hence, we provide a workaround in this case
+
+          Epetra_MultiVector dummy(vector->Map(), 1, false);
+          Epetra_Import data_exchange (dummy.Map(), v.vector->Map());
+
+          int ierr = dummy.Import(*v.vector, data_exchange, Insert);
+          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+          ierr = vector->Update (1.0, dummy, 1.0);
+          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+#endif
+        }
+    }
+
+
+
+    bool
+    Vector::operator == (const Vector &v) const
+    {
+      Assert (size() == v.size(),
+              ExcDimensionMismatch(size(), v.size()));
+      if (local_size() != v.local_size())
+        return false;
+
+      size_type i;
+      for (i=0; i<local_size(); i++)
+        if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
+
+      return true;
+    }
+
+
+
+    bool
+    Vector::operator != (const Vector &v) const
+    {
+      Assert (size() == v.size(),
+              ExcDimensionMismatch(size(), v.size()));
+
+      return (!(*this==v));
+    }
+
+
+
+    bool
+    Vector::all_zero () const
+    {
+      // get a representation of the vector and
+      // loop over all the elements
+      TrilinosScalar *start_ptr = (*vector)[0];
+      const TrilinosScalar *ptr  = start_ptr,
+                            *eptr = start_ptr + local_size();
+      unsigned int flag = 0;
+      while (ptr != eptr)
+        {
+          if (*ptr != 0)
+            {
+              flag = 1;
+              break;
+            }
+          ++ptr;
+        }
+
+#ifdef DEAL_II_WITH_MPI
+      // in parallel, check that the vector
+      // is zero on _all_ processors.
+      const Epetra_MpiComm *mpi_comm
+        = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
+      Assert(mpi_comm != nullptr, ExcInternalError());
+      unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
+      return num_nonzero == 0;
+#else
+      return flag == 0;
+#endif
+
+    }
+
+
+
+    bool
+    Vector::is_non_negative () const
+    {
+#ifdef DEAL_II_WITH_MPI
+      // if this vector is a parallel one, then
+      // we need to communicate to determine
+      // the answer to the current
+      // function. this still has to be
+      // implemented
+      AssertThrow(local_size() == size(), ExcNotImplemented());
+#endif
+      // get a representation of the vector and
+      // loop over all the elements
+      TrilinosScalar *start_ptr;
+      int leading_dimension;
+      int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+      // TODO: This
+      // won't work in parallel like
+      // this. Find out a better way to
+      // this in that case.
+      const TrilinosScalar *ptr  = start_ptr,
+                            *eptr = start_ptr + size();
+      bool flag = true;
+      while (ptr != eptr)
+        {
+          if (*ptr < 0.0)
+            {
+              flag = false;
+              break;
+            }
+          ++ptr;
+        }
+
+      return flag;
+    }
+
+
+
+    void
+    Vector::print (std::ostream      &out,
+                   const unsigned int precision,
+                   const bool         scientific,
+                   const bool         across) const
+    {
+      AssertThrow (out, ExcIO());
+      boost::io::ios_flags_saver restore_flags(out);
+
+      // get a representation of the
+      // vector and loop over all
+      // the elements TODO: up to
+      // now only local data printed
+      // out! Find a way to neatly
+      // output distributed data...
+      TrilinosScalar *val;
+      int leading_dimension;
+      int ierr = vector->ExtractView (&val, &leading_dimension);
+
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      out.precision (precision);
+      if (scientific)
+        out.setf (std::ios::scientific, std::ios::floatfield);
+      else
+        out.setf (std::ios::fixed, std::ios::floatfield);
+
+      if (across)
+        for (size_type i=0; i<size(); ++i)
+          out << static_cast<double>(val[i]) << ' ';
+      else
+        for (size_type i=0; i<size(); ++i)
+          out << static_cast<double>(val[i]) << std::endl;
+      out << std::endl;
+
+      // restore the representation
+      // of the vector
+      AssertThrow (out, ExcIO());
+    }
+
+
+
+    void
+    Vector::swap (Vector &v)
+    {
+      std::swap(last_action, v.last_action);
+      std::swap(compressed, v.compressed);
+      std::swap(vector, v.vector);
+    }
+
+
+
+    std::size_t
+    Vector::memory_consumption () const
+    {
+      //TODO[TH]: No accurate memory
+      //consumption for Trilinos vectors
+      //yet. This is a rough approximation with
+      //one index and the value per local
+      //entry.
+      return sizeof(*this)
+             + this->local_size()*( sizeof(double)+
+                                    sizeof(TrilinosWrappers::types::int_type) );
+    }
   } /* end of namespace MPI */
 } /* end of namespace TrilinosWrappers */
 
+namespace TrilinosWrappers
+{
+  namespace MPI
+  {
+#   include "trilinos_vector.inst"
+  }
+}
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_WITH_TRILINOS
similarity index 89%
rename from source/lac/trilinos_vector_base.inst.in
rename to source/lac/trilinos_vector.inst.in
index 2f1c4d3d3c3837677e9b49a46d878aa17b7c7bd5..e94ff81dc975ce85a59c6f14173a8e5b66c9c867 100644 (file)
@@ -18,6 +18,6 @@
 for (SCALAR : REAL_SCALARS)
 {
     template
-    VectorBase &
-    VectorBase::operator = (const ::dealii::Vector<SCALAR> &v);
+    Vector &
+    Vector::operator = (const ::dealii::Vector<SCALAR> &v);
 }
diff --git a/source/lac/trilinos_vector_base.cc b/source/lac/trilinos_vector_base.cc
deleted file mode 100644 (file)
index ec456bf..0000000
+++ /dev/null
@@ -1,462 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2008 - 2017 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/lac/trilinos_vector_base.h>
-
-#ifdef DEAL_II_WITH_TRILINOS
-
-#  include <deal.II/base/mpi.h>
-
-#  include <cmath>
-
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-#  include <Epetra_Import.h>
-#  include <Epetra_Export.h>
-
-#  include <boost/io/ios_state.hpp>
-DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
-
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace TrilinosWrappers
-{
-  namespace internal
-  {
-    VectorReference::operator TrilinosScalar () const
-    {
-      Assert (index < vector.size(),
-              ExcIndexRange (index, 0, vector.size()));
-
-      // Trilinos allows for vectors to be referenced by the [] or ()
-      // operators but only () checks index bounds. We check these bounds by
-      // ourselves, so we can use []. Note that we can only get local values.
-
-      const TrilinosWrappers::types::int_type local_index =
-        vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
-      Assert (local_index >= 0,
-              VectorBase::ExcAccessToNonLocalElement (index, vector.local_size(),
-                                                      vector.vector->Map().MinMyGID(),
-                                                      vector.vector->Map().MaxMyGID()));
-
-
-      return (*(vector.vector))[0][local_index];
-    }
-  }
-
-
-
-  VectorBase::VectorBase ()
-    :
-    last_action (Zero),
-    compressed  (true),
-    has_ghosts  (false),
-#ifdef DEAL_II_WITH_MPI
-    vector(new Epetra_FEVector(
-             Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_SELF))))
-#else
-    vector(new Epetra_FEVector(
-             Epetra_Map(0,0,Epetra_SerialComm())))
-#endif
-  {}
-
-
-
-  VectorBase::VectorBase (const VectorBase &v)
-    :
-    Subscriptor(),
-    last_action (Zero),
-    compressed (true),
-    has_ghosts  (v.has_ghosts),
-    vector(new Epetra_FEVector(*v.vector)),
-    owned_elements(complete_index_set(v.size()))
-  {}
-
-
-
-  VectorBase::~VectorBase ()
-  {}
-
-
-
-  void
-  VectorBase::clear ()
-  {
-    // When we clear the vector, reset the pointer and generate an empty
-    // vector.
-#ifdef DEAL_II_WITH_MPI
-    Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
-#else
-    Epetra_Map map (0, 0, Epetra_SerialComm());
-#endif
-
-    has_ghosts = false;
-    vector.reset (new Epetra_FEVector(map));
-    last_action = Zero;
-  }
-
-
-
-  VectorBase &
-  VectorBase::operator = (const VectorBase &v)
-  {
-    Assert (vector.get() != nullptr,
-            ExcMessage("Vector is not constructed properly."));
-
-    if (local_range() != v.local_range())
-      {
-        last_action = Zero;
-        vector.reset (new Epetra_FEVector(*v.vector));
-        has_ghosts = v.has_ghosts;
-        owned_elements = v.owned_elements;
-      }
-    else
-      {
-        Assert (vector->Map().SameAs(v.vector->Map()) == true,
-                ExcMessage ("The Epetra maps in the assignment operator ="
-                            " do not match, even though the local_range "
-                            " seems to be the same. Check vector setup!"));
-        int ierr;
-        ierr = vector->GlobalAssemble(last_action);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-        ierr = vector->Update(1.0, *v.vector, 0.0);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-        last_action = Zero;
-      }
-
-    return *this;
-  }
-
-
-
-  template <typename number>
-  VectorBase &
-  VectorBase::operator = (const ::dealii::Vector<number> &v)
-  {
-    Assert (size() == v.size(),
-            ExcDimensionMismatch(size(), v.size()));
-
-    // this is probably not very efficient
-    // but works. in particular, we could do
-    // better if we know that
-    // number==TrilinosScalar because then we
-    // could elide the copying of elements
-    //
-    // let's hope this isn't a
-    // particularly frequent operation
-    std::pair<size_type, size_type>
-    local_range = this->local_range ();
-    for (size_type i=local_range.first; i<local_range.second; ++i)
-      (*vector)[0][i-local_range.first] = v(i);
-
-    return *this;
-  }
-
-
-
-  void
-  VectorBase::compress (::dealii::VectorOperation::values given_last_action)
-  {
-    //Select which mode to send to Trilinos. Note that we use last_action if
-    //available and ignore what the user tells us to detect wrongly mixed
-    //operations. Typically given_last_action is only used on machines that do
-    //not execute an operation (because they have no own cells for example).
-    Epetra_CombineMode mode = last_action;
-    if (last_action == Zero)
-      {
-        if (given_last_action==::dealii::VectorOperation::add)
-          mode = Add;
-        else if (given_last_action==::dealii::VectorOperation::insert)
-          mode = Insert;
-      }
-
-#ifdef DEBUG
-#  ifdef DEAL_II_WITH_MPI
-    // check that every process has decided to use the same mode. This will
-    // otherwise result in undefined behaviour in the call to
-    // GlobalAssemble().
-    double double_mode = mode;
-    const Epetra_MpiComm *comm_ptr
-      = dynamic_cast<const Epetra_MpiComm *>(&(vector_partitioner().Comm()));
-    Assert (comm_ptr != nullptr, ExcInternalError());
-    Utilities::MPI::MinMaxAvg result
-      = Utilities::MPI::min_max_avg (double_mode, comm_ptr->GetMpiComm());
-    Assert(result.max-result.min<1e-5,
-           ExcMessage ("Not all processors agree whether the last operation on "
-                       "this vector was an addition or a set operation. This will "
-                       "prevent the compress() operation from succeeding."));
-
-#  endif
-#endif
-
-    // Now pass over the information about what we did last to the vector.
-    int ierr = 0;
-    if (nonlocal_vector.get() == nullptr || mode != Add)
-      ierr = vector->GlobalAssemble(mode);
-    else
-      {
-        Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
-        ierr = vector->Export(*nonlocal_vector, exporter, mode);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-        ierr = nonlocal_vector->PutScalar(0.);
-      }
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-    last_action = Zero;
-
-    compressed = true;
-  }
-
-
-
-  TrilinosScalar
-  VectorBase::operator () (const size_type index) const
-  {
-    // Extract local indices in the vector.
-    TrilinosWrappers::types::int_type trilinos_i =
-      vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
-    TrilinosScalar value = 0.;
-
-    // If the element is not present on the current processor, we can't
-    // continue. This is the main difference to the el() function.
-    if (trilinos_i == -1)
-      {
-        Assert (false, ExcAccessToNonLocalElement(index, local_size(),
-                                                  vector->Map().MinMyGID(),
-                                                  vector->Map().MaxMyGID()));
-      }
-    else
-      value = (*vector)[0][trilinos_i];
-
-    return value;
-  }
-
-
-
-  void
-  VectorBase::add (const VectorBase &v,
-                   const bool        allow_different_maps)
-  {
-    if (allow_different_maps == false)
-      *this += v;
-    else
-      {
-        Assert (!has_ghost_elements(), ExcGhostsPresent());
-        AssertThrow (size() == v.size(),
-                     ExcDimensionMismatch (size(), v.size()));
-
-#if DEAL_II_TRILINOS_VERSION_GTE(11,11,0)
-        Epetra_Import data_exchange (vector->Map(), v.vector->Map());
-        int ierr = vector->Import(*v.vector, data_exchange, Epetra_AddLocalAlso);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-        last_action = Add;
-#else
-        // In versions older than 11.11 the Import function is broken for adding
-        // Hence, we provide a workaround in this case
-
-        Epetra_MultiVector dummy(vector->Map(), 1, false);
-        Epetra_Import data_exchange (dummy.Map(), v.vector->Map());
-
-        int ierr = dummy.Import(*v.vector, data_exchange, Insert);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-        ierr = vector->Update (1.0, dummy, 1.0);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-#endif
-      }
-  }
-
-
-
-  bool
-  VectorBase::operator == (const VectorBase &v) const
-  {
-    Assert (size() == v.size(),
-            ExcDimensionMismatch(size(), v.size()));
-    if (local_size() != v.local_size())
-      return false;
-
-    size_type i;
-    for (i=0; i<local_size(); i++)
-      if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
-
-    return true;
-  }
-
-
-
-  bool
-  VectorBase::operator != (const VectorBase &v) const
-  {
-    Assert (size() == v.size(),
-            ExcDimensionMismatch(size(), v.size()));
-
-    return (!(*this==v));
-  }
-
-
-
-  bool
-  VectorBase::all_zero () const
-  {
-    // get a representation of the vector and
-    // loop over all the elements
-    TrilinosScalar *start_ptr = (*vector)[0];
-    const TrilinosScalar *ptr  = start_ptr,
-                          *eptr = start_ptr + local_size();
-    unsigned int flag = 0;
-    while (ptr != eptr)
-      {
-        if (*ptr != 0)
-          {
-            flag = 1;
-            break;
-          }
-        ++ptr;
-      }
-
-#ifdef DEAL_II_WITH_MPI
-    // in parallel, check that the vector
-    // is zero on _all_ processors.
-    const Epetra_MpiComm *mpi_comm
-      = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
-    Assert(mpi_comm != nullptr, ExcInternalError());
-    unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
-    return num_nonzero == 0;
-#else
-    return flag == 0;
-#endif
-
-  }
-
-
-
-  bool
-  VectorBase::is_non_negative () const
-  {
-#ifdef DEAL_II_WITH_MPI
-    // if this vector is a parallel one, then
-    // we need to communicate to determine
-    // the answer to the current
-    // function. this still has to be
-    // implemented
-    AssertThrow(local_size() == size(), ExcNotImplemented());
-#endif
-    // get a representation of the vector and
-    // loop over all the elements
-    TrilinosScalar *start_ptr;
-    int leading_dimension;
-    int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    // TODO: This
-    // won't work in parallel like
-    // this. Find out a better way to
-    // this in that case.
-    const TrilinosScalar *ptr  = start_ptr,
-                          *eptr = start_ptr + size();
-    bool flag = true;
-    while (ptr != eptr)
-      {
-        if (*ptr < 0.0)
-          {
-            flag = false;
-            break;
-          }
-        ++ptr;
-      }
-
-    return flag;
-  }
-
-
-
-  void
-  VectorBase::print (std::ostream      &out,
-                     const unsigned int precision,
-                     const bool         scientific,
-                     const bool         across) const
-  {
-    AssertThrow (out, ExcIO());
-    boost::io::ios_flags_saver restore_flags(out);
-
-    // get a representation of the
-    // vector and loop over all
-    // the elements TODO: up to
-    // now only local data printed
-    // out! Find a way to neatly
-    // output distributed data...
-    TrilinosScalar *val;
-    int leading_dimension;
-    int ierr = vector->ExtractView (&val, &leading_dimension);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-    out.precision (precision);
-    if (scientific)
-      out.setf (std::ios::scientific, std::ios::floatfield);
-    else
-      out.setf (std::ios::fixed, std::ios::floatfield);
-
-    if (across)
-      for (size_type i=0; i<size(); ++i)
-        out << static_cast<double>(val[i]) << ' ';
-    else
-      for (size_type i=0; i<size(); ++i)
-        out << static_cast<double>(val[i]) << std::endl;
-    out << std::endl;
-
-    // restore the representation
-    // of the vector
-    AssertThrow (out, ExcIO());
-  }
-
-
-
-  void
-  VectorBase::swap (VectorBase &v)
-  {
-    std::swap(last_action, v.last_action);
-    std::swap(compressed, v.compressed);
-    std::swap(vector, v.vector);
-  }
-
-
-
-  std::size_t
-  VectorBase::memory_consumption () const
-  {
-    //TODO[TH]: No accurate memory
-    //consumption for Trilinos vectors
-    //yet. This is a rough approximation with
-    //one index and the value per local
-    //entry.
-    return sizeof(*this)
-           + this->local_size()*( sizeof(double)+
-                                  sizeof(TrilinosWrappers::types::int_type) );
-  }
-
-} /* end of namespace TrilinosWrappers */
-
-
-namespace TrilinosWrappers
-{
-#include "trilinos_vector_base.inst"
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_WITH_TRILINOS

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.