]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify TrilinosWrappers::SparseMatrix::vmult a bit more by using the same code...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 9 May 2013 13:13:18 +0000 (13:13 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 9 May 2013 13:13:18 +0000 (13:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@29487 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/block_matrix_base.h
deal.II/include/deal.II/lac/parallel_vector.h
deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
tests/mpi/trilinos_matvec_01.cc
tests/mpi/trilinos_matvec_03.cc [new file with mode: 0644]
tests/mpi/trilinos_matvec_03/ncpu_2/cmp/generic [new file with mode: 0644]

index 1f796a01fce3ed59dbe84da9e342d47b9b61264c..f98b791278622de43985f1d2aaef010766cfd5f5 100644 (file)
@@ -2521,13 +2521,9 @@ BlockMatrixBase<MatrixType>::vmult_add (BlockVectorType       &dst,
           ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
 
   for (unsigned int row=0; row<n_block_rows(); ++row)
-    {
-      block(row,0).vmult_add (dst.block(row),
-                              src.block(0));
-      for (unsigned int col=1; col<n_block_cols(); ++col)
-        block(row,col).vmult_add (dst.block(row),
-                                  src.block(col));
-    };
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      block(row,col).vmult_add (dst.block(row),
+                                src.block(col));
 }
 
 
@@ -2628,11 +2624,9 @@ BlockMatrixBase<MatrixType>::Tvmult_add (BlockVectorType &dst,
           ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
 
   for (unsigned int row=0; row<n_block_rows(); ++row)
-    {
-      for (unsigned int col=0; col<n_block_cols(); ++col)
-        block(row,col).Tvmult_add (dst.block(col),
-                                   src.block(row));
-    };
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      block(row,col).Tvmult_add (dst.block(col),
+                                 src.block(row));
 }
 
 
index 088d0528e040e526cde2eb423cb1e9a8f712dbcc..64cade9c6f41794625785b995c861e4f9d4f7372 100644 (file)
@@ -526,31 +526,29 @@ namespace parallel
       bool is_ghost_entry (const types::global_dof_index global_index) const;
 
       /**
-       * Make the @p 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.
+       * Make the @p 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 <i>locally owned</i> elements of this vector.
+       *
+       * It holds that end() - begin() == local_size().
        */
       iterator begin ();
 
       /**
-       * Return constant iterator to the start of
-       * the vector.
+       * 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.
+       * 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.
+       * Return a constant iterator pointing to the element past the end of
+       * the array of the locally owned entries.
        */
       const_iterator end () const;
       //@}
index 5acd553a719d3bdf9c0c2f2fcb4554b9911d6c4b..ff56ed065c51f33478317558ffbe7eddfdef1daa 100644 (file)
@@ -22,8 +22,8 @@
 #  include <deal.II/base/index_set.h>
 #  include <deal.II/lac/full_matrix.h>
 #  include <deal.II/lac/exceptions.h>
-#  include <deal.II/lac/trilinos_vector_base.h>
-#  include <deal.II/lac/parallel_vector.h>
+#  include <deal.II/lac/trilinos_vector.h>
+#  include <deal.II/lac/vector_view.h>
 
 #  include <vector>
 #  include <cmath>
@@ -52,7 +52,6 @@ class SparsityPattern;
 namespace TrilinosWrappers
 {
   // forward declarations
-  class VectorBase;
   class SparseMatrix;
   class SparsityPattern;
 
@@ -1856,6 +1855,11 @@ namespace TrilinosWrappers
      *
      * Source and destination must not be the same vector.
      *
+     * This function can be called with several different vector objects,
+     * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well
+     * as deal.II's own vector classes Vector<double> and
+     * parallel::distributed::Vector<double>.
+     *
      * Note that both vectors have to be distributed vectors generated using
      * the same Map as was used for the matrix in case you work on a
      * distributed memory architecture, using the interface in the
@@ -1866,21 +1870,9 @@ namespace TrilinosWrappers
      * running on one processor, since the matrix object is inherently
      * distributed. Otherwise, and exception will be thrown.
      */
-    void vmult (VectorBase       &dst,
-                const VectorBase &src) const;
-
-    /**
-     * Same as before, but working with deal.II's own distributed vector
-     * class.
-     */
-    void vmult (parallel::distributed::Vector<TrilinosScalar>       &dst,
-                const parallel::distributed::Vector<TrilinosScalar> &src) const;
-
-    /**
-     * Same as before, but working with deal.II's own vector class.
-     */
-    void vmult (dealii::Vector<TrilinosScalar>       &dst,
-                const dealii::Vector<TrilinosScalar> &src) const;
+    template<typename VectorType>
+    void vmult (VectorType       &dst,
+                const VectorType &src) const;
 
     /**
      * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
@@ -1889,6 +1881,11 @@ namespace TrilinosWrappers
      *
      * Source and destination must not be the same vector.
      *
+     * This function can be called with several different vector objects,
+     * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well
+     * as deal.II's own vector classes Vector<double> and
+     * parallel::distributed::Vector<double>.
+     *
      * Note that both vectors have to be distributed vectors generated using
      * the same Map as was used for the matrix in case you work on a
      * distributed memory architecture, using the interface in the
@@ -1899,21 +1896,9 @@ namespace TrilinosWrappers
      * running on one processor, since the matrix object is inherently
      * distributed. Otherwise, and exception will be thrown.
      */
-    void Tvmult (VectorBase       &dst,
-                 const VectorBase &src) const;
-
-    /**
-     * Same as before, but working with deal.II's own distributed vector
-     * class.
-     */
-    void Tvmult (parallel::distributed::Vector<TrilinosScalar>       &dst,
-                 const parallel::distributed::Vector<TrilinosScalar> &src) const;
-
-    /**
-     * Same as before, but working with deal.II's own vector class.
-     */
-    void Tvmult (dealii::Vector<TrilinosScalar>       &dst,
-                 const dealii::Vector<TrilinosScalar> &src) const;
+    template <typename VectorType>
+    void Tvmult (VectorType       &dst,
+                 const VectorType &src) const;
 
     /**
      * Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
@@ -3622,65 +3607,50 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  void
-  SparseMatrix::vmult (VectorBase       &dst,
-                       const VectorBase &src) const
+  namespace internal
   {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-    Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
-    Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
-            ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
-            ExcMessage ("Row map of matrix does not fit with vector map!"));
-
-    const int ierr = matrix->Multiply (false, src.trilinos_vector(),
-                                       dst.trilinos_vector());
-    Assert (ierr == 0, ExcTrilinosError(ierr));
-    (void)ierr; // removes -Wunused-variable in optimized mode
-  }
-
-
-
-  inline
-  void
-  SparseMatrix::vmult (parallel::distributed::Vector<TrilinosScalar>       &dst,
-                       const parallel::distributed::Vector<TrilinosScalar> &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-    Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
-    AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
-    AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
-
-    Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
-    Epetra_Vector tril_src (View, matrix->DomainMap(),
-                            const_cast<double *>(src.begin()));
+    namespace SparseMatrix
+    {
+      template <typename VectorType>
+      inline
+      void check_vector_map_equality(const Epetra_CrsMatrix &,
+                                     const VectorType &,
+                                     const VectorType &)
+      {
+      }
 
-    const int ierr = matrix->Multiply (false, tril_src, tril_dst);
-    Assert (ierr == 0, ExcTrilinosError(ierr));
-    (void)ierr; // removes -Wunused-variable in optimized mode
+      inline
+      void check_vector_map_equality(const Epetra_CrsMatrix              &m,
+                                     const TrilinosWrappers::MPI::Vector &in,
+                                     const TrilinosWrappers::MPI::Vector &out)
+      {
+        Assert (in.vector_partitioner().SameAs(m.DomainMap()) == true,
+                ExcMessage ("Column map of matrix does not fit with vector map!"));
+        Assert (out.vector_partitioner().SameAs(m.RangeMap()) == true,
+                ExcMessage ("Row map of matrix does not fit with vector map!"));
+      }
+    }
   }
 
 
-
+  template <typename VectorType>
   inline
   void
-  SparseMatrix::vmult (dealii::Vector<TrilinosScalar>       &dst,
-                       const dealii::Vector<TrilinosScalar> &src) const
+  SparseMatrix::vmult (VectorType       &dst,
+                       const VectorType &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
     Assert (matrix->Filled(), ExcMatrixNotCompressed());
 
-    AssertDimension (static_cast<unsigned int>(matrix->DomainMap().NumMyElements()),
-                     static_cast<unsigned int>(matrix->DomainMap().NumGlobalElements()));
-    AssertDimension (dst.size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
-    AssertDimension (src.size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
+    internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst);
+    const int dst_local_size = dst.end() - dst.begin();
+    AssertDimension (dst_local_size, matrix->RangeMap().NumMyElements());
+    const int src_local_size = src.end() - src.begin();
+    AssertDimension (src_local_size, matrix->DomainMap().NumMyElements());
 
     Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
     Epetra_Vector tril_src (View, matrix->DomainMap(),
-                            const_cast<double *>(src.begin()));
+                            const_cast<TrilinosScalar *>(src.begin()));
 
     const int ierr = matrix->Multiply (false, tril_src, tril_dst);
     Assert (ierr == 0, ExcTrilinosError(ierr));
@@ -3688,62 +3658,21 @@ namespace TrilinosWrappers
   }
 
 
-
-  inline
-  void
-  SparseMatrix::Tvmult (VectorBase       &dst,
-                        const VectorBase &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-    Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
-    Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
-            ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
-            ExcMessage ("Row map of matrix does not fit with vector map!"));
-
-    const int ierr = matrix->Multiply (true, src.trilinos_vector(),
-                                       dst.trilinos_vector());
-    Assert (ierr == 0, ExcTrilinosError(ierr));
-    (void)ierr; // removes -Wunused-variable in optimized mode
-  }
-
-
-
-  inline
-  void
-  SparseMatrix::Tvmult (parallel::distributed::Vector<TrilinosScalar>      &dst,
-                        const parallel::distributed::Vector<TrilinosScalar> &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-    Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
-    AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
-    AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
-
-    Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
-    Epetra_Vector tril_src (View, matrix->RangeMap(),
-                            const_cast<double *>(src.begin()));
-
-    const int ierr = matrix->Multiply (true, tril_src, tril_dst);
-    Assert (ierr == 0, ExcTrilinosError(ierr));
-    (void)ierr; // removes -Wunused-variable in optimized mode
-  }
-
-
-
+  template <typename VectorType>
   inline
   void
-  SparseMatrix::Tvmult (dealii::Vector<TrilinosScalar>      &dst,
-                        const dealii::Vector<TrilinosScalar> &src) const
+  SparseMatrix::Tvmult (VectorType       &dst,
+                        const VectorType &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
     Assert (matrix->Filled(), ExcMatrixNotCompressed());
 
-    AssertDimension (static_cast<unsigned int>(matrix->DomainMap().NumMyElements()),
-                     static_cast<unsigned int>(matrix->DomainMap().NumGlobalElements()));
-    AssertDimension (dst.size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
-    AssertDimension (src.size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
+    internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src);
+    const int dst_local_size = dst.end() - dst.begin();
+    AssertDimension (dst_local_size, matrix->DomainMap().NumMyElements());
+    const int src_local_size = src.end() - src.begin();
+    AssertDimension (src_local_size, matrix->RangeMap().NumMyElements());
 
     Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
     Epetra_Vector tril_src (View, matrix->RangeMap(),
@@ -3764,14 +3693,18 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
 
-    // Choose to reinit the vector with fast argument set, which does not
-    // overwrite the content -- this is what we needs since we're going to
-    // overwrite that anyway in the vmult operation.
-    VectorType temp_vector;
-    temp_vector.reinit(dst, true);
-
-    vmult (temp_vector, src);
-    dst += temp_vector;
+    // Reinit a temporary vector with fast argument set, which does not
+    // overwrite the content (to save time). However, the
+    // TrilinosWrappers::Vector classes do not support this, so create a
+    // deal.II local vector that has this fast setting. It will be accepted in
+    // vmult because it only checks the local size.
+    dealii::Vector<TrilinosScalar> temp_vector;
+    temp_vector.reinit(dst.end()-dst.begin(), true);
+    dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
+    dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+    vmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
+    if (dst_view.size() > 0)
+      dst_view += temp_vector;
   }
 
 
@@ -3783,11 +3716,19 @@ namespace TrilinosWrappers
                             const VectorType &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
-    VectorType temp_vector;
-    temp_vector.reinit(dst, true);
 
-    Tvmult (temp_vector, src);
-    dst += temp_vector;
+    // Reinit a temporary vector with fast argument set, which does not
+    // overwrite the content (to save time). However, the
+    // TrilinosWrappers::Vector classes do not support this, so create a
+    // deal.II local vector that has this fast setting. It will be accepted in
+    // vmult because it only checks the local size.
+    dealii::Vector<TrilinosScalar> temp_vector;
+    temp_vector.reinit(dst.end()-dst.begin(), true);
+    dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
+    dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+    Tvmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
+    if (dst_view.size() > 0)
+      dst_view += temp_vector;
   }
 
 
index 134bc47064ad90670160120b6f3be424cf97b5ac..f9c2145a9012700a451fd47508c4c09c39340761 100644 (file)
@@ -267,10 +267,12 @@ namespace TrilinosWrappers
      * <tt>C</tt> standard libraries
      * <tt>vector<...></tt> class.
      */
-    typedef TrilinosScalar            value_type;
-    typedef TrilinosScalar            real_type;
-    typedef std::size_t               size_type;
-    typedef internal::VectorReference reference;
+    typedef TrilinosScalar                  value_type;
+    typedef TrilinosScalar                  real_type;
+    typedef std::size_t                     size_type;
+    typedef value_type                     *iterator;
+    typedef const value_type               *const_iterator;
+    typedef internal::VectorReference       reference;
     typedef const internal::VectorReference const_reference;
 
     /**
@@ -659,6 +661,33 @@ namespace TrilinosWrappers
      */
     TrilinosScalar el (const unsigned int index) 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 
+     *
+     * 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;
+
     /**
      * A collective set operation:
      * instead of setting individual
@@ -1232,6 +1261,42 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  typename VectorBase::iterator
+  VectorBase::begin()
+  {
+    return (*vector)[0];
+  }
+
+
+
+  inline
+  typename VectorBase::iterator
+  VectorBase::end()
+  {
+    return (*vector)[0]+local_size();
+  }
+
+
+
+  inline
+  typename VectorBase::const_iterator
+  VectorBase::begin() const
+  {
+    return (*vector)[0];
+  }
+
+
+
+  inline
+  typename VectorBase::const_iterator
+  VectorBase::end() const
+  {
+    return (*vector)[0]+local_size();
+  }
+
+
+
   inline
   void
   VectorBase::reinit (const VectorBase &v,
index f1305de745f61782b2f5ce4222eebd875c8203ef..eb84b5b91b99b47c1500bc1f54664abe398b1542 100644 (file)
@@ -108,7 +108,7 @@ void test ()
                                // (should be no roundoff difference)
   for (unsigned int i=0; i<row_partitioning.n_elements(); ++i)
     {
-      const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+      const unsigned int global_index = row_partitioning.nth_index_in_set(i);
       Assert (dy(global_index) == y(global_index), ExcInternalError());
     }
   if (my_id == 0) deallog << "OK" << std::endl;
diff --git a/tests/mpi/trilinos_matvec_03.cc b/tests/mpi/trilinos_matvec_03.cc
new file mode 100644 (file)
index 0000000..94460fc
--- /dev/null
@@ -0,0 +1,162 @@
+//----------------------  trilinos_matvec_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2011, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------  trilinos_matvec_03.cc  ---------------------------
+
+
+// Test whether TrilinosWrappers::SparseMatrix::(T)vmult(_add) gives same
+// result with Trilinos vector and parallel distributed vector when one
+// processor does not have any rows
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test ()
+{
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  const unsigned int n_rows = 3;
+  const unsigned int n_cols = 4;
+
+  IndexSet row_partitioning (n_rows);
+  IndexSet col_partitioning (n_cols);
+
+  if (n_procs == 1)
+    {
+      row_partitioning.add_range(0, n_rows);
+      col_partitioning.add_range(0, n_cols);
+    }
+  else if (n_procs == 2)
+    {
+      // row_partitioning should be { [0, 2), [2, n_rows) }
+      // col_partitioning should be { [0, 2), [2, n_cols) }
+      // col_relevant_set should be { [0, 3), [1, n_cols) }
+      if (my_id == 0)
+       {
+         row_partitioning.add_range(0, n_rows);
+         col_partitioning.add_range(0, n_cols);
+       }
+    }
+  else
+    Assert (false, ExcNotImplemented());
+
+  TrilinosWrappers::SparsityPattern sp (row_partitioning,
+                                       col_partitioning, MPI_COMM_WORLD);
+  if (my_id == 0)
+    {
+      sp.add (0, 0);
+      sp.add (0, 2);
+      sp.add (2, 3);
+    }
+  sp.compress();
+
+  TrilinosWrappers::SparseMatrix A;
+  A.clear ();
+  A.reinit (sp);
+  if (my_id == 0)
+    {
+      A.add (0, 0, 1);
+      A.add (0, 2, 1);
+      A.add (2, 3, 2.0);
+    }
+  A.compress(VectorOperation::add);
+
+  TrilinosWrappers::MPI::Vector x, y;
+  x.reinit (col_partitioning, MPI_COMM_WORLD);
+  y.reinit (row_partitioning, MPI_COMM_WORLD);
+
+  parallel::distributed::Vector<double>
+    dx (col_partitioning, col_partitioning, MPI_COMM_WORLD),
+    dy (row_partitioning, row_partitioning, MPI_COMM_WORLD);
+
+  for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+    {
+      const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+      dx(global_index) = (double)rand()/RAND_MAX;
+      x(global_index)  = dx(global_index);
+    }
+  dy = 1.;
+
+  A.vmult (y, x);
+  A.vmult (dy, dx);
+
+                               // compare whether we got the same result
+                               // (should be no roundoff difference)
+  for (unsigned int i=0; i<row_partitioning.n_elements(); ++i)
+    {
+      const unsigned int global_index = row_partitioning.nth_index_in_set(i);
+      Assert (dy(global_index) == y(global_index), ExcInternalError());
+    }
+
+  A.vmult_add (y, x);
+  A.vmult_add (dy, dx);
+
+                               // compare whether we got the same result
+                               // (should be no roundoff difference)
+  for (unsigned int i=0; i<row_partitioning.n_elements(); ++i)
+    {
+      const unsigned int global_index = row_partitioning.nth_index_in_set(i);
+      Assert (dy(global_index) == y(global_index), ExcInternalError());
+    }
+
+  A.Tvmult (x, y);
+  A.Tvmult (dx, dy);
+  for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+    {
+      const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+      Assert (dx(global_index) == x(global_index), ExcInternalError());
+    }
+
+  A.Tvmult_add (x, y);
+  A.Tvmult_add (dx, dy);
+  for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+    {
+      const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+      Assert (dx(global_index) == x(global_index), ExcInternalError());
+    }
+
+  if (my_id == 0) deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile(output_file_for_mpi("trilinos_matvec_03").c_str());
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      test();
+    }
+  else
+    test();
+
+}
diff --git a/tests/mpi/trilinos_matvec_03/ncpu_2/cmp/generic b/tests/mpi/trilinos_matvec_03/ncpu_2/cmp/generic
new file mode 100644 (file)
index 0000000..be8d055
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL:0::OK

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.