From: kronbichler Date: Thu, 9 May 2013 13:13:18 +0000 (+0000) Subject: Simplify TrilinosWrappers::SparseMatrix::vmult a bit more by using the same code... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8cdaba8f7724d193321da7d76015fb08554872d3;p=dealii-svn.git Simplify TrilinosWrappers::SparseMatrix::vmult a bit more by using the same code for all different vector types. git-svn-id: https://svn.dealii.org/trunk@29487 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/block_matrix_base.h b/deal.II/include/deal.II/lac/block_matrix_base.h index 1f796a01fc..f98b791278 100644 --- a/deal.II/include/deal.II/lac/block_matrix_base.h +++ b/deal.II/include/deal.II/lac/block_matrix_base.h @@ -2521,13 +2521,9 @@ BlockMatrixBase::vmult_add (BlockVectorType &dst, ExcDimensionMismatch(src.n_blocks(), n_block_cols())); for (unsigned int row=0; row::Tvmult_add (BlockVectorType &dst, ExcDimensionMismatch(src.n_blocks(), n_block_rows())); for (unsigned int row=0; rowvector<> class of the C++ - * standard library by returning - * iterators to the start and end of the - * locally owned elements of this vector. + * Make the @p Vector class a bit like the vector<> class of + * the C++ standard library by returning iterators to the start and end + * of the locally owned elements of this vector. + * + * It holds that end() - begin() == local_size(). */ iterator begin (); /** - * Return constant iterator to the start of - * the vector. + * Return constant iterator to the start of the locally owned elements + * of the vector. */ const_iterator begin () const; /** - * Return an iterator pointing to the - * element past the end of the array of - * locally owned entries. + * Return an iterator pointing to the element past the end of the array + * of locally owned entries. */ iterator end (); /** - * Return a constant iterator pointing to - * the element past the end of the array - * of the locally owned entries. + * Return a constant iterator pointing to the element past the end of + * the array of the locally owned entries. */ const_iterator end () const; //@} 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 5acd553a71..ff56ed065c 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h @@ -22,8 +22,8 @@ # include # include # include -# include -# include +# include +# include # include # include @@ -52,7 +52,6 @@ class SparsityPattern; namespace TrilinosWrappers { // forward declarations - class VectorBase; class SparseMatrix; class SparsityPattern; @@ -1856,6 +1855,11 @@ namespace TrilinosWrappers * * Source and destination must not be the same vector. * + * This function can be called with several different vector objects, + * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well + * as deal.II's own vector classes Vector and + * parallel::distributed::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 @@ -1866,21 +1870,9 @@ namespace TrilinosWrappers * running on one processor, since the matrix object is inherently * distributed. Otherwise, and exception will be thrown. */ - void vmult (VectorBase &dst, - const VectorBase &src) const; - - /** - * Same as before, but working with deal.II's own distributed vector - * class. - */ - void vmult (parallel::distributed::Vector &dst, - const parallel::distributed::Vector &src) const; - - /** - * Same as before, but working with deal.II's own vector class. - */ - void vmult (dealii::Vector &dst, - const dealii::Vector &src) const; + template + void vmult (VectorType &dst, + const VectorType &src) const; /** * Matrix-vector multiplication: let dst = MT*src with @@ -1889,6 +1881,11 @@ namespace TrilinosWrappers * * Source and destination must not be the same vector. * + * This function can be called with several different vector objects, + * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well + * as deal.II's own vector classes Vector and + * parallel::distributed::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 @@ -1899,21 +1896,9 @@ namespace TrilinosWrappers * running on one processor, since the matrix object is inherently * distributed. Otherwise, and exception will be thrown. */ - void Tvmult (VectorBase &dst, - const VectorBase &src) const; - - /** - * Same as before, but working with deal.II's own distributed vector - * class. - */ - void Tvmult (parallel::distributed::Vector &dst, - const parallel::distributed::Vector &src) const; - - /** - * Same as before, but working with deal.II's own vector class. - */ - void Tvmult (dealii::Vector &dst, - const dealii::Vector &src) const; + template + void Tvmult (VectorType &dst, + const VectorType &src) const; /** * Adding matrix-vector multiplication. Add M*src on dst @@ -3622,65 +3607,50 @@ namespace TrilinosWrappers - inline - void - SparseMatrix::vmult (VectorBase &dst, - const VectorBase &src) const + namespace internal { - Assert (&src != &dst, ExcSourceEqualsDestination()); - Assert (matrix->Filled(), ExcMatrixNotCompressed()); - - Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true, - ExcMessage ("Column map of matrix does not fit with vector map!")); - Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true, - ExcMessage ("Row map of matrix does not fit with vector map!")); - - const int ierr = matrix->Multiply (false, src.trilinos_vector(), - dst.trilinos_vector()); - Assert (ierr == 0, ExcTrilinosError(ierr)); - (void)ierr; // removes -Wunused-variable in optimized mode - } - - - - inline - void - SparseMatrix::vmult (parallel::distributed::Vector &dst, - const parallel::distributed::Vector &src) const - { - Assert (&src != &dst, ExcSourceEqualsDestination()); - Assert (matrix->Filled(), ExcMatrixNotCompressed()); - - AssertDimension (dst.local_size(), static_cast(matrix->RangeMap().NumMyElements())); - AssertDimension (src.local_size(), static_cast(matrix->DomainMap().NumMyElements())); - - Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin()); - Epetra_Vector tril_src (View, matrix->DomainMap(), - const_cast(src.begin())); + namespace SparseMatrix + { + template + inline + void check_vector_map_equality(const Epetra_CrsMatrix &, + const VectorType &, + const VectorType &) + { + } - const int ierr = matrix->Multiply (false, tril_src, tril_dst); - Assert (ierr == 0, ExcTrilinosError(ierr)); - (void)ierr; // removes -Wunused-variable in optimized mode + inline + void check_vector_map_equality(const Epetra_CrsMatrix &m, + const TrilinosWrappers::MPI::Vector &in, + const TrilinosWrappers::MPI::Vector &out) + { + Assert (in.vector_partitioner().SameAs(m.DomainMap()) == true, + ExcMessage ("Column map of matrix does not fit with vector map!")); + Assert (out.vector_partitioner().SameAs(m.RangeMap()) == true, + ExcMessage ("Row map of matrix does not fit with vector map!")); + } + } } - + template inline void - SparseMatrix::vmult (dealii::Vector &dst, - const dealii::Vector &src) const + SparseMatrix::vmult (VectorType &dst, + const VectorType &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); Assert (matrix->Filled(), ExcMatrixNotCompressed()); - AssertDimension (static_cast(matrix->DomainMap().NumMyElements()), - static_cast(matrix->DomainMap().NumGlobalElements())); - AssertDimension (dst.size(), static_cast(matrix->RangeMap().NumMyElements())); - AssertDimension (src.size(), static_cast(matrix->DomainMap().NumMyElements())); + internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst); + const int dst_local_size = dst.end() - dst.begin(); + AssertDimension (dst_local_size, matrix->RangeMap().NumMyElements()); + const int src_local_size = src.end() - src.begin(); + AssertDimension (src_local_size, matrix->DomainMap().NumMyElements()); Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin()); Epetra_Vector tril_src (View, matrix->DomainMap(), - const_cast(src.begin())); + const_cast(src.begin())); const int ierr = matrix->Multiply (false, tril_src, tril_dst); Assert (ierr == 0, ExcTrilinosError(ierr)); @@ -3688,62 +3658,21 @@ namespace TrilinosWrappers } - - inline - void - SparseMatrix::Tvmult (VectorBase &dst, - const VectorBase &src) const - { - Assert (&src != &dst, ExcSourceEqualsDestination()); - Assert (matrix->Filled(), ExcMatrixNotCompressed()); - - Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true, - ExcMessage ("Column map of matrix does not fit with vector map!")); - Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true, - ExcMessage ("Row map of matrix does not fit with vector map!")); - - const int ierr = matrix->Multiply (true, src.trilinos_vector(), - dst.trilinos_vector()); - Assert (ierr == 0, ExcTrilinosError(ierr)); - (void)ierr; // removes -Wunused-variable in optimized mode - } - - - - inline - void - SparseMatrix::Tvmult (parallel::distributed::Vector &dst, - const parallel::distributed::Vector &src) const - { - Assert (&src != &dst, ExcSourceEqualsDestination()); - Assert (matrix->Filled(), ExcMatrixNotCompressed()); - - AssertDimension (dst.local_size(), static_cast(matrix->DomainMap().NumMyElements())); - AssertDimension (src.local_size(), static_cast(matrix->RangeMap().NumMyElements())); - - Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin()); - Epetra_Vector tril_src (View, matrix->RangeMap(), - const_cast(src.begin())); - - const int ierr = matrix->Multiply (true, tril_src, tril_dst); - Assert (ierr == 0, ExcTrilinosError(ierr)); - (void)ierr; // removes -Wunused-variable in optimized mode - } - - - + + template inline void - SparseMatrix::Tvmult (dealii::Vector &dst, - const dealii::Vector &src) const + SparseMatrix::Tvmult (VectorType &dst, + const VectorType &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); Assert (matrix->Filled(), ExcMatrixNotCompressed()); - AssertDimension (static_cast(matrix->DomainMap().NumMyElements()), - static_cast(matrix->DomainMap().NumGlobalElements())); - AssertDimension (dst.size(), static_cast(matrix->DomainMap().NumMyElements())); - AssertDimension (src.size(), static_cast(matrix->RangeMap().NumMyElements())); + internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src); + const int dst_local_size = dst.end() - dst.begin(); + AssertDimension (dst_local_size, matrix->DomainMap().NumMyElements()); + const int src_local_size = src.end() - src.begin(); + AssertDimension (src_local_size, matrix->RangeMap().NumMyElements()); Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin()); Epetra_Vector tril_src (View, matrix->RangeMap(), @@ -3764,14 +3693,18 @@ namespace TrilinosWrappers { Assert (&src != &dst, ExcSourceEqualsDestination()); - // Choose to reinit the vector with fast argument set, which does not - // overwrite the content -- this is what we needs since we're going to - // overwrite that anyway in the vmult operation. - VectorType temp_vector; - temp_vector.reinit(dst, true); - - vmult (temp_vector, src); - dst += temp_vector; + // Reinit a temporary vector with fast argument set, which does not + // overwrite the content (to save time). However, the + // TrilinosWrappers::Vector classes do not support this, so create a + // deal.II local vector that has this fast setting. It will be accepted in + // vmult because it only checks the local size. + dealii::Vector temp_vector; + temp_vector.reinit(dst.end()-dst.begin(), true); + dealii::VectorView src_view(src.end()-src.begin(), src.begin()); + dealii::VectorView dst_view(dst.end()-dst.begin(), dst.begin()); + vmult (temp_vector, static_cast&>(src_view)); + if (dst_view.size() > 0) + dst_view += temp_vector; } @@ -3783,11 +3716,19 @@ namespace TrilinosWrappers const VectorType &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); - VectorType temp_vector; - temp_vector.reinit(dst, true); - Tvmult (temp_vector, src); - dst += temp_vector; + // Reinit a temporary vector with fast argument set, which does not + // overwrite the content (to save time). However, the + // TrilinosWrappers::Vector classes do not support this, so create a + // deal.II local vector that has this fast setting. It will be accepted in + // vmult because it only checks the local size. + dealii::Vector temp_vector; + temp_vector.reinit(dst.end()-dst.begin(), true); + dealii::VectorView src_view(src.end()-src.begin(), src.begin()); + dealii::VectorView dst_view(dst.end()-dst.begin(), dst.begin()); + Tvmult (temp_vector, static_cast&>(src_view)); + if (dst_view.size() > 0) + dst_view += temp_vector; } diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index 134bc47064..f9c2145a90 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -267,10 +267,12 @@ namespace TrilinosWrappers * C standard libraries * vector<...> class. */ - typedef TrilinosScalar value_type; - typedef TrilinosScalar real_type; - typedef std::size_t size_type; - typedef internal::VectorReference reference; + typedef TrilinosScalar value_type; + typedef TrilinosScalar real_type; + typedef std::size_t size_type; + typedef value_type *iterator; + typedef const value_type *const_iterator; + typedef internal::VectorReference reference; typedef const internal::VectorReference const_reference; /** @@ -659,6 +661,33 @@ namespace TrilinosWrappers */ TrilinosScalar el (const unsigned int index) const; + /** + * Make the Vector class a bit like the vector<> class of + * the C++ standard library by returning iterators to the start and end + * of the locally owned elements of this vector. The ordering of local elements corresponds to the one given + * + * It holds that end() - begin() == local_size(). + */ + iterator begin (); + + /** + * Return constant iterator to the start of the locally owned elements + * of the vector. + */ + const_iterator begin () const; + + /** + * Return an iterator pointing to the element past the end of the array + * of locally owned entries. + */ + iterator end (); + + /** + * Return a constant iterator pointing to the element past the end of + * the array of the locally owned entries. + */ + const_iterator end () const; + /** * A collective set operation: * instead of setting individual @@ -1232,6 +1261,42 @@ namespace TrilinosWrappers + inline + typename VectorBase::iterator + VectorBase::begin() + { + return (*vector)[0]; + } + + + + inline + typename VectorBase::iterator + VectorBase::end() + { + return (*vector)[0]+local_size(); + } + + + + inline + typename VectorBase::const_iterator + VectorBase::begin() const + { + return (*vector)[0]; + } + + + + inline + typename VectorBase::const_iterator + VectorBase::end() const + { + return (*vector)[0]+local_size(); + } + + + inline void VectorBase::reinit (const VectorBase &v, diff --git a/tests/mpi/trilinos_matvec_01.cc b/tests/mpi/trilinos_matvec_01.cc index f1305de745..eb84b5b91b 100644 --- a/tests/mpi/trilinos_matvec_01.cc +++ b/tests/mpi/trilinos_matvec_01.cc @@ -108,7 +108,7 @@ void test () // (should be no roundoff difference) for (unsigned int i=0; i +#include +#include +#include +#include +#include +#include +#include +#include + + +void test () +{ + const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + const unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + const unsigned int n_rows = 3; + const unsigned int n_cols = 4; + + IndexSet row_partitioning (n_rows); + IndexSet col_partitioning (n_cols); + + if (n_procs == 1) + { + row_partitioning.add_range(0, n_rows); + col_partitioning.add_range(0, n_cols); + } + else if (n_procs == 2) + { + // row_partitioning should be { [0, 2), [2, n_rows) } + // col_partitioning should be { [0, 2), [2, n_cols) } + // col_relevant_set should be { [0, 3), [1, n_cols) } + if (my_id == 0) + { + row_partitioning.add_range(0, n_rows); + col_partitioning.add_range(0, n_cols); + } + } + else + Assert (false, ExcNotImplemented()); + + TrilinosWrappers::SparsityPattern sp (row_partitioning, + col_partitioning, MPI_COMM_WORLD); + if (my_id == 0) + { + sp.add (0, 0); + sp.add (0, 2); + sp.add (2, 3); + } + sp.compress(); + + TrilinosWrappers::SparseMatrix A; + A.clear (); + A.reinit (sp); + if (my_id == 0) + { + A.add (0, 0, 1); + A.add (0, 2, 1); + A.add (2, 3, 2.0); + } + A.compress(VectorOperation::add); + + TrilinosWrappers::MPI::Vector x, y; + x.reinit (col_partitioning, MPI_COMM_WORLD); + y.reinit (row_partitioning, MPI_COMM_WORLD); + + parallel::distributed::Vector + dx (col_partitioning, col_partitioning, MPI_COMM_WORLD), + dy (row_partitioning, row_partitioning, MPI_COMM_WORLD); + + for (unsigned int i=0; i