From 4cbb72f66fbeb9ae3e72e25261a1bdaf63758b46 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 8 May 2013 20:11:04 +0000 Subject: [PATCH] Implement a few more matrix-vector products with Trilinos matrices and deal.II vectors. Avoid usage of temporary arrays in Trilinos sparse matrix add/set functoin to make it really thread safe. git-svn-id: https://svn.dealii.org/trunk@29483 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 8 + .../lac/trilinos_block_sparse_matrix.h | 359 +++++---------- .../deal.II/lac/trilinos_sparse_matrix.h | 408 ++++++++---------- .../trilinos/block_sparse_matrix_vector_01.cc | 126 ++++++ .../block_sparse_matrix_vector_01/cmp/generic | 2 + .../trilinos/block_sparse_matrix_vector_02.cc | 120 ++++++ .../block_sparse_matrix_vector_02/cmp/generic | 2 + .../trilinos/block_sparse_matrix_vector_03.cc | 120 ++++++ .../block_sparse_matrix_vector_03/cmp/generic | 2 + .../trilinos/block_sparse_matrix_vector_04.cc | 112 +++++ .../block_sparse_matrix_vector_04/cmp/generic | 2 + tests/trilinos/sparse_matrix_vector_08.cc | 108 +++++ .../sparse_matrix_vector_08/cmp/generic | 2 + tests/trilinos/sparse_matrix_vector_09.cc | 108 +++++ .../sparse_matrix_vector_09/cmp/generic | 2 + tests/trilinos/sparse_matrix_vector_10.cc | 109 +++++ .../sparse_matrix_vector_10/cmp/generic | 2 + tests/trilinos/sparse_matrix_vector_11.cc | 108 +++++ .../sparse_matrix_vector_11/cmp/generic | 2 + 19 files changed, 1214 insertions(+), 488 deletions(-) create mode 100644 tests/trilinos/block_sparse_matrix_vector_01.cc create mode 100644 tests/trilinos/block_sparse_matrix_vector_01/cmp/generic create mode 100644 tests/trilinos/block_sparse_matrix_vector_02.cc create mode 100644 tests/trilinos/block_sparse_matrix_vector_02/cmp/generic create mode 100644 tests/trilinos/block_sparse_matrix_vector_03.cc create mode 100644 tests/trilinos/block_sparse_matrix_vector_03/cmp/generic create mode 100644 tests/trilinos/block_sparse_matrix_vector_04.cc create mode 100644 tests/trilinos/block_sparse_matrix_vector_04/cmp/generic create mode 100644 tests/trilinos/sparse_matrix_vector_08.cc create mode 100644 tests/trilinos/sparse_matrix_vector_08/cmp/generic create mode 100644 tests/trilinos/sparse_matrix_vector_09.cc create mode 100644 tests/trilinos/sparse_matrix_vector_09/cmp/generic create mode 100644 tests/trilinos/sparse_matrix_vector_10.cc create mode 100644 tests/trilinos/sparse_matrix_vector_10/cmp/generic create mode 100644 tests/trilinos/sparse_matrix_vector_11.cc create mode 100644 tests/trilinos/sparse_matrix_vector_11/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 2f39d127f7..ace3025a09 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -117,6 +117,14 @@ this function.
    +
  1. New: The classes TrilinosWrappers::SparseMatrix and +TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with +deal.II's own vector classes Vector and +parallel::distributed::Vector. +
    +(Martin Kronbichler, 2013/05/08) +
    2. +
  2. Improved: The matrix-vector product ChunkSparseMatrix::vmult now runs in parallel in shared memory.
    diff --git a/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h index f298dad0cc..1550c8a246 100644 --- a/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h @@ -18,6 +18,7 @@ #ifdef DEAL_II_WITH_TRILINOS # include +# include # include # include # include @@ -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 + 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 + 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 + void vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const dealii::internal::bool2type, + const dealii::internal::bool2type) const; + + /** + * Internal version of (T)vmult where the source vector is a block vector + * but the destination vector is a non-block vector + */ + template + void vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const dealii::internal::bool2type, + const dealii::internal::bool2type) const; + + /** + * Internal version of (T)vmult where the source vector is a non-block + * vector but the destination vector is a block vector + */ + template + void vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const dealii::internal::bool2type, + const dealii::internal::bool2type) 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 + void vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + const dealii::internal::bool2type, + const dealii::internal::bool2type) const; }; @@ -666,142 +574,97 @@ namespace TrilinosWrappers + template 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::value>(), + dealii::internal::bool2type::value>()); } + template 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::value>(), + dealii::internal::bool2type::value>()); } + template 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, + dealii::internal::bool2type) 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 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, + dealii::internal::bool2type) const { - BaseClass::Tvmult_nonblock_block (dst, src); + if (transpose == true) + BaseClass::Tvmult_nonblock_block (dst, src); + else + BaseClass::vmult_nonblock_block (dst, src); } + template inline void - BlockSparseMatrix::Tvmult (Vector &dst, - const BlockVector &src) const + BlockSparseMatrix::vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + dealii::internal::bool2type, + dealii::internal::bool2type) const { - BaseClass::Tvmult_nonblock_block (dst, src); + if (transpose == true) + BaseClass::Tvmult_block_nonblock (dst, src); + else + BaseClass::vmult_block_nonblock (dst, src); } + template inline void - BlockSparseMatrix::Tvmult (VectorBase &dst, - const VectorBase &src) const + BlockSparseMatrix::vmult (VectorType1 &dst, + const VectorType2 &src, + const bool transpose, + dealii::internal::bool2type, + dealii::internal::bool2type) const { - BaseClass::Tvmult_nonblock_nonblock (dst, src); + if (transpose == true) + BaseClass::Tvmult_nonblock_nonblock (dst, src); + else + BaseClass::vmult_nonblock_nonblock (dst, src); } } diff --git a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h index daf89a6c57..5acd553a71 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h @@ -1851,171 +1851,123 @@ namespace TrilinosWrappers //@{ /** - * Matrix-vector multiplication: - * let dst = M*src with - * M being this matrix. + * Matrix-vector multiplication: let dst = M*src with M + * 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 &dst, const parallel::distributed::Vector &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 &dst, const dealii::Vector &src) const; /** - * Matrix-vector multiplication: - * let dst = - * MT*src with - * M being this - * matrix. This function does the - * same as vmult() but takes the - * transposed matrix. + * Matrix-vector multiplication: let dst = MT*src with + * M 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 &dst, const parallel::distributed::Vector &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 &dst, const dealii::Vector &src) const; /** - * Adding Matrix-vector - * multiplication. Add - * M*src on dst - * with M being this - * matrix. + * Adding matrix-vector multiplication. Add M*src on dst + * with M 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 and + * parallel::distributed::Vector. + * + * 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), 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 + void vmult_add (VectorType &dst, + const VectorType &src) const; /** - * Adding Matrix-vector - * multiplication. Add - * MT*src to - * dst with M being - * this matrix. This function - * does the same as vmult_add() - * but takes the transposed - * matrix. + * Adding matrix-vector multiplication. Add MT*src to + * dst with M 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 and + * parallel::distributed::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. + * 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), 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 + 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 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 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 long_val_array; + std::vector 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(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(row)) == true) { if (matrix->Filled() == false) { ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_columns, - const_cast(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(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 long_val_array; + std::vector 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(values); n_columns = n_cols; #ifdef DEBUG for (unsigned int j=0; j 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(row)) == true) { ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns, - const_cast(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 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 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 index 0000000000..e9948ba514 --- /dev/null +++ b/tests/trilinos/block_sparse_matrix_vector_01.cc @@ -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 +#include +#include +#include +#include +#include + + +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 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 +#include +#include +#include +#include +#include + + +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 w(2); + w.block(0).reinit(size_1); + w.block(1).reinit(size_1); + w.collect_sizes(); + Vectorv(size_2); + for (unsigned int i=0; i +#include +#include +#include +#include +#include + + +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 v(2); + v.block(0).reinit(size_2); + v.block(1).reinit(size_2); + v.collect_sizes(); + Vectorw(size_1); + for (unsigned int i=0; i +#include +#include +#include +#include +#include + + +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 w(size_1), v(size_2); + for (unsigned int i=0; i +#include +#include +#include +#include +#include + + +void test (Vector &v, + Vector &w) +{ + TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size()); + for (unsigned int i=0; i v (100); + Vector 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_08/cmp/generic @@ -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 index 0000000000..a600928b8a --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_09.cc @@ -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 +#include +#include +#include +#include +#include + + +void test (Vector &v, + Vector &w) +{ + TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size()); + for (unsigned int i=0; i v (95); + Vector 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_09/cmp/generic @@ -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 index 0000000000..693b567b3a --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_10.cc @@ -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 +#include +#include +#include +#include +#include + + +void test (parallel::distributed::Vector &v, + parallel::distributed::Vector &w) +{ + TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size()); + for (unsigned int i=0; i v (100); + parallel::distributed::Vector 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_10/cmp/generic @@ -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 index 0000000000..450e8f05c9 --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_11.cc @@ -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 +#include +#include +#include +#include +#include + + +void test (parallel::distributed::Vector &v, + parallel::distributed::Vector &w) +{ + TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size()); + for (unsigned int i=0; i v (95); + parallel::distributed::Vector 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/sparse_matrix_vector_11/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5