]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a few more matrix-vector products with Trilinos matrices and deal.II vector...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 May 2013 20:11:04 +0000 (20:11 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 May 2013 20:11:04 +0000 (20:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@29483 0785d39b-7218-0410-832d-ea1e28bc413d

19 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
tests/trilinos/block_sparse_matrix_vector_01.cc [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_01/cmp/generic [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_02.cc [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_02/cmp/generic [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_03.cc [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_03/cmp/generic [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_04.cc [new file with mode: 0644]
tests/trilinos/block_sparse_matrix_vector_04/cmp/generic [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_08.cc [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_08/cmp/generic [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_09.cc [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_09/cmp/generic [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_10.cc [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_10/cmp/generic [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_11.cc [new file with mode: 0644]
tests/trilinos/sparse_matrix_vector_11/cmp/generic [new file with mode: 0644]

index 2f39d127f731a16b29e54b21cfa08ed83bbec706..ace3025a0996bed73679f690e9ce89c5190a85fc 100644 (file)
@@ -117,6 +117,14 @@ this function.
 
 <ol>
 
+<li> New: The classes TrilinosWrappers::SparseMatrix and
+TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with
+deal.II's own vector classes Vector<double> and
+parallel::distributed::Vector<double>.
+<br>
+(Martin Kronbichler, 2013/05/08)
+</li>
+
 <li> Improved: The matrix-vector product ChunkSparseMatrix::vmult now runs in
 parallel in shared memory.
 <br>
index f298dad0cc7ef420f0bbee7277ff84bbff0e9684..1550c8a246ddd73aba7eefd4460acd7e77326bf8 100644 (file)
@@ -18,6 +18,7 @@
 #ifdef DEAL_II_WITH_TRILINOS
 
 #  include <deal.II/base/table.h>
+#  include <deal.II/base/template_constraints.h>
 #  include <deal.II/lac/block_matrix_base.h>
 #  include <deal.II/lac/trilinos_sparse_matrix.h>
 #  include <deal.II/lac/trilinos_block_vector.h>
@@ -292,161 +293,23 @@ namespace TrilinosWrappers
     unsigned int n_nonzero_elements () const;
 
     /**
-     * Matrix-vector multiplication:
-     * let $dst = M*src$ with $M$
-     * being this matrix.
+     * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this
+     * matrix. The vector types can be block vectors or non-block vectors
+     * (only if the matrix has only one row or column, respectively), and need
+     * to define TrilinosWrappers::SparseMatrix::vmult.
      */
-    void vmult (MPI::BlockVector       &dst,
-                const MPI::BlockVector &src) const;
-
-
-    /**
-     * Matrix-vector multiplication:
-     * let $dst = M*src$ with $M$
-     * being this matrix, now applied
-     * to localized block vectors
-     * (works only when run on one
-     * processor).
-     */
-    void vmult (BlockVector       &dst,
-                const BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block column.
-     */
-    void vmult (MPI::BlockVector  &dst,
-                const MPI::Vector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block column, now
-     * applied to localized vectors
-     * (works only when run on one
-     * processor).
-     */
-    void vmult (BlockVector  &dst,
-                const Vector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block row.
-     */
-    void vmult (MPI::Vector            &dst,
-                const MPI::BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block row, now
-     * applied to localized vectors
-     * (works only when run on one
-     * processor).
-     */
-    void vmult (Vector            &dst,
-                const BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block.
-     */
-    void vmult (VectorBase       &dst,
-                const VectorBase &src) const;
-
-    /**
-     * Matrix-vector multiplication:
-     * let $dst = M^T*src$ with $M$
-     * being this matrix. This
-     * function does the same as
-     * vmult() but takes the
-     * transposed matrix.
-     */
-    void Tvmult (MPI::BlockVector       &dst,
-                 const MPI::BlockVector &src) const;
-
-    /**
-     * Matrix-vector multiplication:
-     * let $dst = M^T*src$ with $M$
-     * being this matrix. This
-     * function does the same as
-     * vmult() but takes the
-     * transposed matrix, now applied
-     * to localized Trilinos vectors
-     * (works only when run on one
-     * processor).
-     */
-    void Tvmult (BlockVector       &dst,
-                 const BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block row.
-     */
-    void Tvmult (MPI::BlockVector  &dst,
-                 const MPI::Vector &src) const;
+    template <typename VectorType1, typename VectorType2>
+    void vmult (VectorType1       &dst,
+                const VectorType2 &src) const;
 
     /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block row, now
-     * applied to localized Trilinos
-     * vectors (works only when run
-     * on one processor).
-     */
-    void Tvmult (BlockVector  &dst,
-                 const Vector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block column.
-     */
-    void Tvmult (MPI::Vector    &dst,
-                 const MPI::BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block column, now
-     * applied to localized Trilinos
-     * vectors (works only when run
-     * on one processor).
-     */
-    void Tvmult (Vector    &dst,
-                 const BlockVector &src) const;
-
-    /**
-     * Matrix-vector
-     * multiplication. Just like the
-     * previous function, but only
-     * applicable if the matrix has
-     * only one block.
+     * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this
+     * matrix. This function does the same as vmult() but takes the transposed
+     * matrix.
      */
-    void Tvmult (VectorBase       &dst,
-                 const VectorBase &src) const;
+    template <typename VectorType1, typename VectorType2>
+    void Tvmult (VectorType1       &dst,
+                 const VectorType2 &src) const;
 
     /**
      * Compute the residual of an
@@ -623,6 +486,51 @@ namespace TrilinosWrappers
                     << "The blocks [" << arg1 << ',' << arg2 << "] and ["
                     << arg3 << ',' << arg4 << "] have differing column numbers.");
     ///@}
+
+  private:
+    /**
+     * Internal version of (T)vmult with two block vectors
+     */
+    template <typename VectorType1, typename VectorType2>
+    void vmult (VectorType1       &dst,
+                const VectorType2 &src,
+                const bool         transpose,
+                const dealii::internal::bool2type<true>,
+                const dealii::internal::bool2type<true>) const;
+
+    /**
+     * Internal version of (T)vmult where the source vector is a block vector
+     * but the destination vector is a non-block vector
+     */
+    template <typename VectorType1, typename VectorType2>
+    void vmult (VectorType1       &dst,
+                const VectorType2 &src,
+                const bool         transpose,
+                const dealii::internal::bool2type<false>,
+                const dealii::internal::bool2type<true>) const;
+
+    /**
+     * Internal version of (T)vmult where the source vector is a non-block
+     * vector but the destination vector is a block vector
+     */
+    template <typename VectorType1, typename VectorType2>
+    void vmult (VectorType1       &dst,
+                const VectorType2 &src,
+                const bool         transpose,
+                const dealii::internal::bool2type<true>,
+                const dealii::internal::bool2type<false>) const;
+
+    /**
+     * Internal version of (T)vmult where both source vector and the
+     * destination vector are non-block vectors (only defined if the matrix
+     * consists of only one block)
+     */
+    template <typename VectorType1, typename VectorType2>
+    void vmult (VectorType1       &dst,
+                const VectorType2 &src,
+                const bool         transpose,
+                const dealii::internal::bool2type<false>,
+                const dealii::internal::bool2type<false>) const;
   };
 
 
@@ -666,142 +574,97 @@ namespace TrilinosWrappers
 
 
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::vmult (MPI::BlockVector       &dst,
-                            const MPI::BlockVector &src) const
+  BlockSparseMatrix::vmult (VectorType1       &dst,
+                            const VectorType2 &src) const
   {
-    BaseClass::vmult_block_block (dst, src);
+    vmult(dst, src, false,
+          dealii::internal::bool2type<IsBlockVector<VectorType1>::value>(),
+          dealii::internal::bool2type<IsBlockVector<VectorType2>::value>());
   }
 
 
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::vmult (BlockVector       &dst,
-                            const BlockVector &src) const
+  BlockSparseMatrix::Tvmult (VectorType1       &dst,
+                             const VectorType2 &src) const
   {
-    BaseClass::vmult_block_block (dst, src);
+    vmult(dst, src, true,
+          dealii::internal::bool2type<IsBlockVector<VectorType1>::value>(),
+          dealii::internal::bool2type<IsBlockVector<VectorType2>::value>());
   }
 
 
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::vmult (MPI::BlockVector  &dst,
-                            const MPI::Vector &src) const
+  BlockSparseMatrix::vmult (VectorType1       &dst,
+                            const VectorType2 &src,
+                            const bool         transpose,
+                            dealii::internal::bool2type<true>,
+                            dealii::internal::bool2type<true>) const
   {
-    BaseClass::vmult_block_nonblock (dst, src);
+    if (transpose == true)
+      BaseClass::Tvmult_block_block (dst, src);
+    else
+      BaseClass::vmult_block_block (dst, src);
   }
 
 
 
-  inline
-  void
-  BlockSparseMatrix::vmult (BlockVector  &dst,
-                            const Vector &src) const
-  {
-    BaseClass::vmult_block_nonblock (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::vmult (MPI::Vector            &dst,
-                            const MPI::BlockVector &src) const
-  {
-    BaseClass::vmult_nonblock_block (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::vmult (Vector            &dst,
-                            const BlockVector &src) const
-  {
-    BaseClass::vmult_nonblock_block (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::vmult (VectorBase       &dst,
-                            const VectorBase &src) const
-  {
-    BaseClass::vmult_nonblock_nonblock (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::Tvmult (MPI::BlockVector       &dst,
-                             const MPI::BlockVector &src) const
-  {
-    BaseClass::Tvmult_block_block (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::Tvmult (BlockVector       &dst,
-                             const BlockVector &src) const
-  {
-    BaseClass::Tvmult_block_block (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::Tvmult (MPI::BlockVector  &dst,
-                             const MPI::Vector &src) const
-  {
-    BaseClass::Tvmult_block_nonblock (dst, src);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::Tvmult (BlockVector  &dst,
-                             const Vector &src) const
-  {
-    BaseClass::Tvmult_block_nonblock (dst, src);
-  }
-
-
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::Tvmult (MPI::Vector            &dst,
-                             const MPI::BlockVector &src) const
+  BlockSparseMatrix::vmult (VectorType1       &dst,
+                            const VectorType2 &src,
+                            const bool         transpose,
+                            dealii::internal::bool2type<false>,
+                            dealii::internal::bool2type<true>) const
   {
-    BaseClass::Tvmult_nonblock_block (dst, src);
+    if (transpose == true)
+      BaseClass::Tvmult_nonblock_block (dst, src);
+    else
+      BaseClass::vmult_nonblock_block (dst, src);
   }
 
 
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::Tvmult (Vector            &dst,
-                             const BlockVector &src) const
+  BlockSparseMatrix::vmult (VectorType1       &dst,
+                            const VectorType2 &src,
+                            const bool         transpose,
+                            dealii::internal::bool2type<true>,
+                            dealii::internal::bool2type<false>) const
   {
-    BaseClass::Tvmult_nonblock_block (dst, src);
+    if (transpose == true)
+      BaseClass::Tvmult_block_nonblock (dst, src);
+    else
+      BaseClass::vmult_block_nonblock (dst, src);
   }
 
 
 
+  template <typename VectorType1, typename VectorType2>
   inline
   void
-  BlockSparseMatrix::Tvmult (VectorBase       &dst,
-                             const VectorBase &src) const
+  BlockSparseMatrix::vmult (VectorType1       &dst,
+                            const VectorType2 &src,
+                            const bool         transpose,
+                            dealii::internal::bool2type<false>,
+                            dealii::internal::bool2type<false>) const
   {
-    BaseClass::Tvmult_nonblock_nonblock (dst, src);
+    if (transpose == true)
+      BaseClass::Tvmult_nonblock_nonblock (dst, src);
+    else
+      BaseClass::vmult_nonblock_nonblock (dst, src);
   }
 
 }
index daf89a6c57cbd3efd574b9b35bb4dd15795985cf..5acd553a719d3bdf9c0c2f2fcb4554b9911d6c4b 100644 (file)
@@ -1851,171 +1851,123 @@ namespace TrilinosWrappers
 //@{
 
     /**
-     * Matrix-vector multiplication:
-     * let <i>dst = M*src</i> with
-     * <i>M</i> being this matrix.
+     * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i>
+     * being this matrix.
      *
-     * Source and destination must
-     * not be the same vector.
+     * Source and destination must not be the same vector.
      *
-     * 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
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::Vector).
+     * 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
+     * TrilinosWrappers::VectorBase class (or one of the two derived classes
+     * Vector and MPI::Vector).
      *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * In case of a localized Vector, this function will only work when
+     * 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
+     * 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.
+     * Same as before, but working with deal.II's own vector class.
      */
     void vmult (dealii::Vector<TrilinosScalar>       &dst,
                 const dealii::Vector<TrilinosScalar> &src) const;
 
     /**
-     * Matrix-vector multiplication:
-     * let <i>dst =
-     * M<sup>T</sup>*src</i> with
-     * <i>M</i> being this
-     * matrix. This function does the
-     * same as vmult() but takes the
-     * transposed matrix.
+     * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
+     * <i>M</i> being this matrix. This function does the same as vmult() but
+     * takes the transposed matrix.
      *
-     * Source and destination must
-     * not be the same vector.
+     * Source and destination must not be the same vector.
      *
-     * 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
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::Vector).
+     * 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
+     * TrilinosWrappers::VectorBase class (or one of the two derived classes
+     * Vector and MPI::Vector).
      *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * In case of a localized Vector, this function will only work when
+     * 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
+     * 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.
+     * Same as before, but working with deal.II's own vector class.
      */
     void Tvmult (dealii::Vector<TrilinosScalar>       &dst,
                  const dealii::Vector<TrilinosScalar> &src) const;
 
     /**
-     * Adding Matrix-vector
-     * multiplication. Add
-     * <i>M*src</i> on <i>dst</i>
-     * with <i>M</i> being this
-     * matrix.
+     * Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
+     * with <i>M</i> being this matrix.
      *
-     * Source and destination must
-     * not be the same vector.
+     * Source and destination must not be the same vector.
      *
-     * 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
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::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>.
+     *
+     * When using a vector of type TrilinosWrappers::MPI::Vector, both vectors
+     * have to be distributed vectors generated using the same Map as was used
+     * for the matrix rows and columns in case you work on a distributed
+     * memory architecture, using the interface in the
+     * TrilinosWrappers::VectorBase class.
+     *
+     * In case of a localized Vector (i.e., TrilinosWrappers::Vector or
+     * Vector<double>), this function will only work when running on one
+     * processor, since the matrix object is inherently
+     * distributed. Otherwise, and exception will be thrown.
      *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
      */
-    void vmult_add (VectorBase       &dst,
-                    const VectorBase &src) const;
+    template<typename VectorType>
+    void vmult_add (VectorType       &dst,
+                    const VectorType &src) const;
 
     /**
-     * Adding Matrix-vector
-     * multiplication. Add
-     * <i>M<sup>T</sup>*src</i> to
-     * <i>dst</i> with <i>M</i> being
-     * this matrix. This function
-     * does the same as vmult_add()
-     * but takes the transposed
-     * matrix.
+     * Adding matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
+     * <i>dst</i> with <i>M</i> being this matrix. This function does the same
+     * as vmult_add() but takes the transposed matrix.
      *
-     * Source and destination must
-     * not be the same vector.
+     * Source and destination must not be the same vector.
      *
-     * 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
-     * TrilinosWrappers::VectorBase
-     * class (or one of the two
-     * derived classes Vector and
-     * MPI::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>.
      *
-     * In case of a localized Vector,
-     * this function will only work
-     * when running on one processor,
-     * since the matrix object is
-     * inherently
-     * distributed. Otherwise, and
-     * exception will be thrown.
+     * When using a vector of type TrilinosWrappers::MPI::Vector, both vectors
+     * have to be distributed vectors generated using the same Map as was used
+     * for the matrix rows and columns in case you work on a distributed
+     * memory architecture, using the interface in the
+     * TrilinosWrappers::VectorBase class.
+     *
+     * In case of a localized Vector (i.e., TrilinosWrappers::Vector or
+     * Vector<double>), this function will only work when running on one
+     * processor, since the matrix object is inherently
+     * distributed. Otherwise, and exception will be thrown.
      */
-    void Tvmult_add (VectorBase       &dst,
-                     const VectorBase &src) const;
+    template <typename VectorType>
+    void Tvmult_add (VectorType       &dst,
+                     const VectorType &src) const;
 
     /**
      * Return the square of the norm
@@ -2585,33 +2537,6 @@ namespace TrilinosWrappers
      */
     bool compressed;
 
-    /**
-     * An internal Trilinos vector that
-     * is used for accelerating vmult_add
-     * functions (in order not to need to
-     * recreate temporary vectors every
-     * time that function is called).
-     */
-    mutable VectorBase temp_vector;
-
-    /**
-     * An internal array of integer
-     * values that is used to store the
-     * column indices when
-     * adding/inserting local data into
-     * the (large) sparse matrix.
-     */
-    std::vector<unsigned int> column_indices;
-
-    /**
-     * An internal array of double values
-     * that is used to store the column
-     * indices when adding/inserting
-     * local data into the (large) sparse
-     * matrix.
-     */
-    std::vector<TrilinosScalar> column_values;
-
     /**
      *  To allow calling protected
      *  prepare_add() and
@@ -3187,26 +3112,39 @@ namespace TrilinosWrappers
     last_action = Insert;
 
     int *col_index_ptr;
-    TrilinosScalar const *col_value_ptr;
+    TrilinosScalar *col_value_ptr;
     int n_columns;
 
-    // If we don't elide zeros, the pointers
-    // are already available...
+    TrilinosScalar short_val_array[100];
+    int short_index_array[100];
+    std::vector<TrilinosScalar> long_val_array;
+    std::vector<int> long_index_array;
+
+
+    // If we don't elide zeros, the pointers are already available... need to
+    // cast to non-const pointers as that is the format taken by Trilinos (but
+    // we will not modify const data)
     if (elide_zero_values == false)
       {
         col_index_ptr = (int *)col_indices;
-        col_value_ptr = values;
+        col_value_ptr = const_cast<TrilinosScalar*>(values);
         n_columns = n_cols;
       }
     else
       {
-        // Otherwise, extract nonzero values in
-        // each row and get the respective
-        // indices.
-        if (column_indices.size() < n_cols)
+        // Otherwise, extract nonzero values in each row and get the
+        // respective indices.
+        if (n_cols > 100)
           {
-            column_indices.resize(n_cols);
-            column_values.resize(n_cols);
+            long_val_array.resize(n_cols);
+            long_index_array.resize(n_cols);
+            col_index_ptr = &long_index_array[0];
+            col_value_ptr = &long_val_array[0];
+          }
+        else
+          {
+            col_index_ptr = &short_index_array[0];
+            col_value_ptr = &short_val_array[0];
           }
 
         n_columns = 0;
@@ -3216,66 +3154,51 @@ namespace TrilinosWrappers
             Assert (numbers::is_finite(value), ExcNumberNotFinite());
             if (value != 0)
               {
-                column_indices[n_columns] = col_indices[j];
-                column_values[n_columns] = value;
+                col_index_ptr[n_columns] = col_indices[j];
+                col_value_ptr[n_columns] = value;
                 n_columns++;
               }
           }
 
         Assert(n_columns <= (int)n_cols, ExcInternalError());
-
-        col_index_ptr = (int *)&column_indices[0];
-        col_value_ptr = &column_values[0];
       }
 
 
-    // If the calling matrix owns the row to
-    // which we want to insert values, we
-    // can directly call the Epetra_CrsMatrix
-    // input function, which is much faster
-    // than the Epetra_FECrsMatrix
-    // function. We distinguish between two
-    // cases: the first one is when the matrix
-    // is not filled (i.e., it is possible to
-    // add new elements to the sparsity pattern),
-    // and the second one is when the pattern is
-    // already fixed. In the former case, we
-    // add the possibility to insert new values,
-    // and in the second we just replace
+    // If the calling matrix owns the row to which we want to insert values,
+    // we can directly call the Epetra_CrsMatrix input function, which is much
+    // faster than the Epetra_FECrsMatrix function. We distinguish between two
+    // cases: the first one is when the matrix is not filled (i.e., it is
+    // possible to add new elements to the sparsity pattern), and the second
+    // one is when the pattern is already fixed. In the former case, we add
+    // the possibility to insert new values, and in the second we just replace
     // data.
     if (row_partitioner().MyGID(static_cast<int>(row)) == true)
       {
         if (matrix->Filled() == false)
           {
             ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_columns,
-                                                                const_cast<double *>(col_value_ptr),
+                                                                col_value_ptr,
                                                                 col_index_ptr);
 
-            // When inserting elements, we do
-            // not want to create exceptions in
-            // the case when inserting non-local
-            // data (since that's what we want
-            // to do right now).
+            // When inserting elements, we do not want to create exceptions in
+            // the case when inserting non-local data (since that's what we
+            // want to do right now).
             if (ierr > 0)
               ierr = 0;
           }
         else
           ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_columns,
-                                                               const_cast<double *>(col_value_ptr),
+                                                               col_value_ptr,
                                                                col_index_ptr);
       }
     else
       {
-        // When we're at off-processor data, we
-        // have to stick with the standard
-        // Insert/ReplaceGlobalValues
-        // function. Nevertheless, the way we
-        // call it is the fastest one (any other
-        // will lead to repeated allocation and
-        // deallocation of memory in order to
-        // call the function we already use,
-        // which is very unefficient if writing
-        // one element at a time).
+        // When we're at off-processor data, we have to stick with the
+        // standard Insert/ReplaceGlobalValues function. Nevertheless, the way
+        // we call it is the fastest one (any other will lead to repeated
+        // allocation and deallocation of memory in order to call the function
+        // we already use, which is very unefficient if writing one element at
+        // a time).
         compressed = false;
 
         if (matrix->Filled() == false)
@@ -3310,16 +3233,12 @@ namespace TrilinosWrappers
 
     if (value == 0)
       {
-        // we have to do checkings on Insert/Add
-        // in any case
-        // to be consistent with the MPI
-        // communication model (see the comments
-        // in the documentation of
-        // TrilinosWrappers::Vector), but we can
-        // save some work if the addend is
-        // zero. However, these actions are done
-        // in case we pass on to the other
-        // function.
+        // we have to do checkings on Insert/Add in any case to be consistent
+        // with the MPI communication model (see the comments in the
+        // documentation of TrilinosWrappers::Vector), but we can save some
+        // work if the addend is zero. However, these actions are done in case
+        // we pass on to the other function.
+
         // TODO: fix this (do not run compress here, but fail)
         if (last_action == Insert)
           {
@@ -3416,15 +3335,21 @@ namespace TrilinosWrappers
     last_action = Add;
 
     int *col_index_ptr;
-    TrilinosScalar const *col_value_ptr;
+    TrilinosScalar *col_value_ptr;
     int n_columns;
 
-    // If we don't elide zeros, the pointers
-    // are already available...
+    double short_val_array[100];
+    int short_index_array[100];
+    std::vector<TrilinosScalar> long_val_array;
+    std::vector<int> long_index_array;
+
+    // If we don't elide zeros, the pointers are already available... need to
+    // cast to non-const pointers as that is the format taken by Trilinos (but
+    // we will not modify const data)
     if (elide_zero_values == false)
       {
         col_index_ptr = (int *)col_indices;
-        col_value_ptr = values;
+        col_value_ptr = const_cast<TrilinosScalar*>(values);
         n_columns = n_cols;
 #ifdef DEBUG
         for (unsigned int j=0; j<n_cols; ++j)
@@ -3433,56 +3358,56 @@ namespace TrilinosWrappers
       }
     else
       {
-        // Otherwise, extract nonzero values in
-        // each row and the corresponding index.
-        if (column_indices.size() < n_cols)
+        // Otherwise, extract nonzero values in each row and the corresponding
+        // index.
+        if (n_cols > 100)
+          {
+            long_val_array.resize(n_cols);
+            long_index_array.resize(n_cols);
+            col_index_ptr = &long_index_array[0];
+            col_value_ptr = &long_val_array[0];
+          }
+        else
           {
-            column_indices.resize(n_cols);
-            column_values.resize(n_cols);
+            col_index_ptr = &short_index_array[0];
+            col_value_ptr = &short_val_array[0];
           }
 
         n_columns = 0;
         for (unsigned int j=0; j<n_cols; ++j)
           {
             const double value = values[j];
+
             Assert (numbers::is_finite(value), ExcNumberNotFinite());
             if (value != 0)
               {
-                column_indices[n_columns] = col_indices[j];
-                column_values[n_columns] = value;
+                col_index_ptr[n_columns] = col_indices[j];
+                col_value_ptr[n_columns] = value;
                 n_columns++;
               }
           }
 
         Assert(n_columns <= (int)n_cols, ExcInternalError());
 
-        col_index_ptr = (int *)&column_indices[0];
-        col_value_ptr = &column_values[0];
       }
 
-    // If the calling matrix owns the row to
-    // which we want to add values, we
-    // can directly call the Epetra_CrsMatrix
-    // input function, which is much faster
-    // than the Epetra_FECrsMatrix function.
+    // If the calling matrix owns the row to which we want to add values, we
+    // can directly call the Epetra_CrsMatrix input function, which is much
+    // faster than the Epetra_FECrsMatrix function.
     if (row_partitioner().MyGID(static_cast<int>(row)) == true)
       {
         ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns,
-                                                             const_cast<double *>(col_value_ptr),
+                                                             col_value_ptr,
                                                              col_index_ptr);
       }
     else
       {
-        // When we're at off-processor data, we
-        // have to stick with the standard
-        // SumIntoGlobalValues
-        // function. Nevertheless, the way we
-        // call it is the fastest one (any other
-        // will lead to repeated allocation and
-        // deallocation of memory in order to
-        // call the function we already use,
-        // which is very unefficient if writing
-        // one element at a time).
+        // When we're at off-processor data, we have to stick with the
+        // standard SumIntoGlobalValues function. Nevertheless, the way we
+        // call it is the fastest one (any other will lead to repeated
+        // allocation and deallocation of memory in order to call the function
+        // we already use, which is very inefficient if writing one element at
+        // a time).
         compressed = false;
 
         ierr = matrix->SumIntoGlobalValues (1, (int *)&row, n_columns,
@@ -3831,18 +3756,18 @@ namespace TrilinosWrappers
 
 
 
+  template <typename VectorType>
   inline
   void
-  SparseMatrix::vmult_add (VectorBase       &dst,
-                           const VectorBase &src) const
+  SparseMatrix::vmult_add (VectorType       &dst,
+                           const VectorType &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
 
-    // Choose to reinit the vector with fast
-    // argument set, which does not overwrite
-    // the content -- this is what we need
-    // since we're going to overwrite that
-    // anyway in the vmult operation.
+    // 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);
@@ -3851,13 +3776,14 @@ namespace TrilinosWrappers
 
 
 
+  template <typename VectorType>
   inline
   void
-  SparseMatrix::Tvmult_add (VectorBase       &dst,
-                            const VectorBase &src) const
+  SparseMatrix::Tvmult_add (VectorType       &dst,
+                            const VectorType &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
-
+    VectorType temp_vector;
     temp_vector.reinit(dst, true);
 
     Tvmult (temp_vector, src);
@@ -3873,7 +3799,8 @@ namespace TrilinosWrappers
     Assert (row_partitioner().SameAs(domain_partitioner()),
             ExcNotQuadratic());
 
-    temp_vector.reinit(v);
+    VectorBase temp_vector;
+    temp_vector.reinit(v, true);
 
     vmult (temp_vector, v);
     return temp_vector*v;
@@ -3889,7 +3816,8 @@ namespace TrilinosWrappers
     Assert (row_partitioner().SameAs(domain_partitioner()),
             ExcNotQuadratic());
 
-    temp_vector.reinit(v);
+    VectorBase temp_vector;
+    temp_vector.reinit(v, true);
 
     vmult (temp_vector, v);
     return u*temp_vector;
diff --git a/tests/trilinos/block_sparse_matrix_vector_01.cc b/tests/trilinos/block_sparse_matrix_vector_01.cc
new file mode 100644 (file)
index 0000000..e9948ba
--- /dev/null
@@ -0,0 +1,126 @@
+//----------------------------  trilinos_block_sparse_matrix_vector_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 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_block_sparse_matrix_vector_01.cc  ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II block vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+  TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+  m.compress (VectorOperation::insert);
+
+  TrilinosWrappers::BlockSparseMatrix m_block;
+  m_block.reinit(2,2);
+  m_block.block(0,0).copy_from(m);
+  m_block.block(1,0).copy_from(m);
+  m_block.block(0,1).copy_from(m);
+  m_block.block(1,1).copy_from(m);
+  m_block.collect_sizes();
+
+  BlockVector<double> v(2), w(2);
+  v.block(0).reinit(size_2);
+  v.block(1).reinit(size_2);
+  v.collect_sizes();
+  w.block(0).reinit(size_1);
+  w.block(1).reinit(size_1);
+  w.collect_sizes();
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+
+                                   // w:=Mv
+  m_block.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*(j+m.n());
+      Assert (w(i) == result, ExcInternalError());
+      Assert (w(i+m.m()) == result, ExcInternalError());
+    }
+
+  for (unsigned int i=0; i<w.size(); ++i)
+    w(i) = i;
+
+  m_block.Tvmult (v, w);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*(j+m.m());
+      Assert (v(i) == result, ExcInternalError());
+      Assert (v(i+m.n()) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("block_sparse_matrix_vector_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      test (50,47);
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/block_sparse_matrix_vector_01/cmp/generic b/tests/trilinos/block_sparse_matrix_vector_01/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/block_sparse_matrix_vector_02.cc b/tests/trilinos/block_sparse_matrix_vector_02.cc
new file mode 100644 (file)
index 0000000..2effc6a
--- /dev/null
@@ -0,0 +1,120 @@
+//----------------------------  trilinos_block_sparse_matrix_vector_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 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_block_sparse_matrix_vector_02.cc  ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II block vector/vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+  TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+  m.compress (VectorOperation::insert);
+
+  TrilinosWrappers::BlockSparseMatrix m_block;
+  m_block.reinit(2,1);
+  m_block.block(0,0).copy_from(m);
+  m_block.block(1,0).copy_from(m);
+  m_block.collect_sizes();
+
+  BlockVector<double> w(2);
+  w.block(0).reinit(size_1);
+  w.block(1).reinit(size_1);
+  w.collect_sizes();
+  Vector<double>v(size_2);
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+
+                                   // w:=Mv
+  m_block.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result, ExcInternalError());
+      Assert (w(i+m.m()) == result, ExcInternalError());
+    }
+
+  for (unsigned int i=0; i<w.size(); ++i)
+    w(i) = i;
+
+  m_block.Tvmult (v, w);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*(j+m.m());
+      Assert (v(i) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("block_sparse_matrix_vector_02/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      test (50,47);
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/block_sparse_matrix_vector_02/cmp/generic b/tests/trilinos/block_sparse_matrix_vector_02/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/block_sparse_matrix_vector_03.cc b/tests/trilinos/block_sparse_matrix_vector_03.cc
new file mode 100644 (file)
index 0000000..9f7a24e
--- /dev/null
@@ -0,0 +1,120 @@
+//----------------------------  trilinos_block_sparse_matrix_vector_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 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_block_sparse_matrix_vector_03.cc  ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II vector/block vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+  TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+  m.compress (VectorOperation::insert);
+
+  TrilinosWrappers::BlockSparseMatrix m_block;
+  m_block.reinit(1,2);
+  m_block.block(0,0).copy_from(m);
+  m_block.block(0,1).copy_from(m);
+  m_block.collect_sizes();
+
+  BlockVector<double> v(2);
+  v.block(0).reinit(size_2);
+  v.block(1).reinit(size_2);
+  v.collect_sizes();
+  Vector<double>w(size_1);
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+
+                                   // w:=Mv
+  m_block.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*(j+m.n());
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  for (unsigned int i=0; i<w.size(); ++i)
+    w(i) = i;
+
+  m_block.Tvmult (v, w);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (v(i) == result, ExcInternalError());
+      Assert (v(i+m.n()) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("block_sparse_matrix_vector_03/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      test (50,47);
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/block_sparse_matrix_vector_03/cmp/generic b/tests/trilinos/block_sparse_matrix_vector_03/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/block_sparse_matrix_vector_04.cc b/tests/trilinos/block_sparse_matrix_vector_04.cc
new file mode 100644 (file)
index 0000000..f1183eb
--- /dev/null
@@ -0,0 +1,112 @@
+//----------------------------  trilinos_block_sparse_matrix_vector_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 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_block_sparse_matrix_vector_04.cc  ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II vector/vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+  TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+  m.compress (VectorOperation::insert);
+
+  TrilinosWrappers::BlockSparseMatrix m_block;
+  m_block.reinit(1,1);
+  m_block.block(0,0).copy_from(m);
+  m_block.collect_sizes();
+
+  Vector<double> w(size_1), v(size_2);
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+
+                                   // w:=Mv
+  m_block.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  for (unsigned int i=0; i<w.size(); ++i)
+    w(i) = i;
+
+  m_block.Tvmult (v, w);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (v(i) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("block_sparse_matrix_vector_04/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      test (50,47);
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/block_sparse_matrix_vector_04/cmp/generic b/tests/trilinos/block_sparse_matrix_vector_04/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/sparse_matrix_vector_08.cc b/tests/trilinos/sparse_matrix_vector_08.cc
new file mode 100644 (file)
index 0000000..0030444
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  trilinos_sparse_matrix_vector_08.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 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_sparse_matrix_vector_08.cc  ---------------------------
+
+
+// check SparseMatrix::vmult, vmult_add with deal.II vector
+
+#include "../tests.h" 
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (Vector<double> &v,
+           Vector<double> &w)
+{
+  TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size());
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  m.compress (VectorOperation::insert);
+
+                                   // w:=Mv
+  m.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  m.vmult_add (w, v);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("sparse_matrix_vector_08/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10); 
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      {
+        Vector<double> v (100);
+        Vector<double> w (95);
+        test (v,w);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/sparse_matrix_vector_08/cmp/generic b/tests/trilinos/sparse_matrix_vector_08/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/sparse_matrix_vector_09.cc b/tests/trilinos/sparse_matrix_vector_09.cc
new file mode 100644 (file)
index 0000000..a600928
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  trilinos_sparse_matrix_vector_09.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 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_sparse_matrix_vector_09.cc  ---------------------------
+
+
+// check SparseMatrix::Tvmult, Tvmult_add with deal.II vector
+
+#include "../tests.h" 
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (Vector<double> &v,
+           Vector<double> &w)
+{
+  TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size());
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  m.compress (VectorOperation::insert);
+
+                                   // w:=Mv
+  m.Tvmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  m.Tvmult_add (w, v);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == result+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("sparse_matrix_vector_09/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10); 
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      {
+        Vector<double> v (95);
+        Vector<double> w (100);
+        test (v,w);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/sparse_matrix_vector_09/cmp/generic b/tests/trilinos/sparse_matrix_vector_09/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/sparse_matrix_vector_10.cc b/tests/trilinos/sparse_matrix_vector_10.cc
new file mode 100644 (file)
index 0000000..693b567
--- /dev/null
@@ -0,0 +1,109 @@
+//----------------------------  trilinos_sparse_matrix_vector_10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 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_sparse_matrix_vector_10.cc  ---------------------------
+
+
+// check SparseMatrix::vmult, vmult_add with distributed deal.II vector (but
+// without using distributed things)
+
+#include "../tests.h" 
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (parallel::distributed::Vector<double> &v,
+           parallel::distributed::Vector<double> &w)
+{
+  TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size());
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  m.compress (VectorOperation::insert);
+
+                                   // w:=Mv
+  m.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  m.vmult_add (w, v);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("sparse_matrix_vector_10/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10); 
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      {
+        parallel::distributed::Vector<double> v (100);
+        parallel::distributed::Vector<double> w (95);
+        test (v,w);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/sparse_matrix_vector_10/cmp/generic b/tests/trilinos/sparse_matrix_vector_10/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/sparse_matrix_vector_11.cc b/tests/trilinos/sparse_matrix_vector_11.cc
new file mode 100644 (file)
index 0000000..450e8f0
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  trilinos_sparse_matrix_vector_11.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 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_sparse_matrix_vector_11.cc  ---------------------------
+
+
+// check SparseMatrix::Tvmult, Tvmult_add with distributed deal.II vector
+
+#include "../tests.h" 
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (parallel::distributed::Vector<double> &v,
+           parallel::distributed::Vector<double> &w)
+{
+  TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size());
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  m.compress (VectorOperation::insert);
+
+                                   // w:=Mv
+  m.Tvmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  m.Tvmult_add (w, v);
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<m.n(); ++i)
+    {
+      double result = 0;
+      for (unsigned int j=0; j<m.m(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == result+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("sparse_matrix_vector_11/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10); 
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+  try
+    {
+      {
+        parallel::distributed::Vector<double> v (95);
+        parallel::distributed::Vector<double> w (100);
+        test (v,w);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/sparse_matrix_vector_11/cmp/generic b/tests/trilinos/sparse_matrix_vector_11/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::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.