From f609c91d539c2d8713e2475c2d29e89d5b4bd9d1 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 18 Oct 2013 12:19:18 +0000 Subject: [PATCH] Improve handling of ghost elements in parallel vector: Distinguish a state with read access to ghost elements and one with write access (assembly mode). git-svn-id: https://svn.dealii.org/trunk@31301 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 25 ++ .../deal.II/lac/parallel_block_vector.h | 50 ++- deal.II/include/deal.II/lac/parallel_vector.h | 353 ++++++++++++------ .../deal.II/lac/parallel_vector.templates.h | 168 +++++---- .../include/deal.II/matrix_free/matrix_free.h | 129 +++++-- tests/mpi/matrix_free_02.cc | 2 - tests/mpi/parallel_block_vector_01.cc | 3 +- tests/mpi/parallel_block_vector_02.cc | 162 ++++++++ .../ncpu_1/cmp/generic | 3 + .../ncpu_10/cmp/generic | 3 + .../ncpu_4/cmp/generic | 3 + 11 files changed, 672 insertions(+), 229 deletions(-) create mode 100644 tests/mpi/parallel_block_vector_02.cc create mode 100644 tests/mpi/parallel_block_vector_02/ncpu_1/cmp/generic create mode 100644 tests/mpi/parallel_block_vector_02/ncpu_10/cmp/generic create mode 100644 tests/mpi/parallel_block_vector_02/ncpu_4/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 54d0a2412b..c399798fed 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -24,6 +24,23 @@ inconvenience this causes.

    +
  1. + Changed: The ghost handling of the parallel::distributed::Vector class has + been reworked: The vector now carries a global state that stores whether + ghost elements have been updated or not. If a vector has ghost elements, it + does not allow calls to compress() any more. Instead, a compress operation + can now only be done when the ghost entries have been cleared before by + calling zero_out_ghosts() or operator=0. The state can be queried by the new + method has_ghost_elements(). This change avoids spurious entries to be + inserted with compress(), but requires some change in user codes. The + behavior of a ghosted vector is now very similar to ghosted PETSc and + Trilinos vectors. The only difference is that the same vector can + also be used as a non-ghosted vector which is designed for use in assembly + routines. +
    + (Martin Kronbichler, 2013/10/18) +
    2. +
  2. Removed: GridTools::collect_periodic_face_pairs. This function is superseded by GridTools::collect_periodic_faces which exports an @@ -145,6 +162,14 @@ inconvenience this causes.

    Specific improvements

      +
    1. + New: parallel::distributed::BlockVector has now methods update_ghost_values, + compress, set_out_ghosts, and has_ghost_elements that do the respective + operation on each block of parallel::distributed::Vector. +
      + (Martin Kronbichler, 2013/10/18) +
      2. +
    2. Fixed: When deriving from DataOut to filter the cells where output is generated, there were two different bugs that result in segmentation faults or wrong cells written (example, step-18).
      diff --git a/deal.II/include/deal.II/lac/parallel_block_vector.h b/deal.II/include/deal.II/lac/parallel_block_vector.h index 3be362a635..7e3f1bb190 100644 --- a/deal.II/include/deal.II/lac/parallel_block_vector.h +++ b/deal.II/include/deal.II/lac/parallel_block_vector.h @@ -261,6 +261,21 @@ namespace parallel */ void update_ghost_values () const; + /** + * This method zeros the entries on ghost dofs, but does not touch + * locally owned DoFs. + * + * After calling this method, read access to ghost elements of the + * vector is forbidden and an exception is thrown. Only write access to + * ghost elements is allowed in this state. + */ + void zero_out_ghosts (); + + /** + * Returns if this Vector contains ghost elements. + */ + bool has_ghost_elements() const; + /** * Return whether the vector contains only elements with value * zero. This function is mainly for internal consistency checks and @@ -584,8 +599,12 @@ namespace parallel void BlockVector::compress (::dealii::VectorOperation::values operation) { + // start all requests for all blocks before finishing the transfers as + // this saves repeated synchronizations for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).compress(operation); + this->block(block).compress_start(block*10 + 8273, operation); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).compress_finish(operation); } @@ -596,7 +615,34 @@ namespace parallel BlockVector::update_ghost_values () const { for (unsigned int block=0; blockn_blocks(); ++block) - this->block(block).update_ghost_values(); + this->block(block).update_ghost_values_start(block*10 + 9923); + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).update_ghost_values_finish(); + } + + + + template + inline + void + BlockVector::zero_out_ghosts () + { + for (unsigned int block=0; blockn_blocks(); ++block) + this->block(block).zero_out_ghosts(); + } + + + + template + inline + bool + BlockVector::has_ghost_elements () const + { + bool has_ghost_elements = false; + for (unsigned int block=0; blockn_blocks(); ++block) + if (this->block(block).has_ghost_elements() == true) + has_ghost_elements = true; + return has_ghost_elements; } diff --git a/deal.II/include/deal.II/lac/parallel_vector.h b/deal.II/include/deal.II/lac/parallel_vector.h index 7e1d853e52..c568a4aba7 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.h +++ b/deal.II/include/deal.II/lac/parallel_vector.h @@ -83,6 +83,22 @@ namespace parallel * - Of course, reduction operations (like norms) make use of collective * all-to-all MPI communications. * + * This vector can take two different states with respect to ghost + * elements: + * - After creation and whenever zero_out_ghosts() is called (or + * operator = (0.)), the vector does only allow writing + * into ghost elements but not reading from ghost elements. + * - After a call to update_ghost_values(), the vector does not allow + * writing into ghost elements but only reading from them. This is in + * order to avoid undesired ghost data artifacts when calling compress() + * after modifying some vector entries. + * The current statues of the ghost entries (read mode or write mode) can + * be queried by the method has_ghost_elements(), which returns + * true exactly when ghost elements have been updated and + * false otherwise, irrespective of the actual number of + * ghost entries in the vector layout (for that information, use + * n_ghost_entries() instead). + * * @author Katharina Kormann, Martin Kronbichler, 2010, 2011 */ template @@ -307,6 +323,16 @@ namespace parallel * ghost data is changed. This is needed to allow functions with a @p * const vector to perform the data exchange without creating * temporaries. + * + * After calling this method, write access to ghost elements of the + * vector is forbidden and an exception is thrown. Only read access to + * ghost elements is allowed in this state. Note that all subsequent + * operations on this vector, like global vector addition, etc., will + * also update the ghost values by a call to this method after the + * operation. However, global reduction operations like norms or the + * inner product will always ignore ghost elements in order to avoid + * counting the ghost data more than once. To allow writing to ghost + * elements again, call zero_out_ghosts(). */ void update_ghost_values () const; @@ -332,7 +358,14 @@ namespace parallel * the communication to finish. Once it is finished, add or set the data * (depending on the flag operation) to the respective positions in the * owning processor, and clear the contents in the ghost data - * fields. The meaning of this argument is the same as in compress(). + * fields. The meaning of this argument is the same as in + * compress(). + * + * This function should be called exactly once per vector after calling + * compress_start, otherwise the result is undefined. In particular, it + * is not well-defined to call compress_start on the same vector again + * before compress_finished has been called. However, there is no + * warning to prevent this situation. * * Must follow a call to the @p compress_start function. */ @@ -372,9 +405,23 @@ namespace parallel /** * This method zeros the entries on ghost dofs, but does not touch * locally owned DoFs. + * + * After calling this method, read access to ghost elements of the + * vector is forbidden and an exception is thrown. Only write access to + * ghost elements is allowed in this state. */ void zero_out_ghosts (); + /** + * Returns whether the vector currently is in a state where ghost values + * can be read or not. This is the same functionality as other parallel + * vectors have. If this method returns false, this only means that + * read-access to ghost elements is prohibited whereas write access is + * still possible (to those entries specified as ghosts during + * initialization), not that there are no ghost elements at all. + */ + bool has_ghost_elements() const; + /** * Return whether the vector contains only elements with value * zero. This function is mainly for internal consistency checks and @@ -909,6 +956,15 @@ namespace parallel */ mutable Number *import_data; + /** + * Stores whether the vector currently allows for reading ghost elements + * or not. Note that this is to ensure consistent ghost data and does + * not indicate whether the vector actually can store ghost elements. In + * particular, when assembling a vector we do not allow reading + * elements, only writing them. + */ + mutable bool vector_is_ghosted; + /** * Provide this class with all functionality of ::dealii::Vector by * creating a VectorView object. @@ -977,6 +1033,7 @@ namespace parallel allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) {} @@ -990,6 +1047,7 @@ namespace parallel allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) { reinit (v, true); @@ -1007,6 +1065,7 @@ namespace parallel allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) { reinit (local_range, ghost_indices, communicator); @@ -1017,11 +1076,12 @@ namespace parallel template inline Vector::Vector (const IndexSet &local_range, - const MPI_Comm communicator) + const MPI_Comm communicator) : allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) { IndexSet ghost_indices(local_range.size()); @@ -1037,6 +1097,7 @@ namespace parallel allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) { reinit (size, false); @@ -1052,6 +1113,7 @@ namespace parallel allocated_size (0), val (0), import_data (0), + vector_is_ghosted (false), vector_view (0, static_cast(0)) { reinit (partitioner); @@ -1083,11 +1145,14 @@ namespace parallel { Assert (c.partitioner.get() != 0, ExcNotInitialized()); - // check whether the two vectors use the same - // parallel partitioner. if not, check if all - // local ranges are the same (that way, we can - // exchange data between different parallel - // layouts) + // we update ghost values whenever one of the input or output vector + // already held ghost values or when we import data from a vector with + // the same local range but different ghost layout + bool must_update_ghost_values = true; + + // check whether the two vectors use the same parallel partitioner. if + // not, check if all local ranges are the same (that way, we can + // exchange data between different parallel layouts) if (partitioner.get() == 0) reinit (c, true); else if (partitioner.get() != c.partitioner.get()) @@ -1101,8 +1166,12 @@ namespace parallel local_ranges_different_loc) reinit (c, true); } + else + must_update_ghost_values = vector_is_ghosted || c.vector_is_ghosted; + vector_view = c.vector_view; - update_ghost_values(); + if (must_update_ghost_values) + update_ghost_values(); return *this; } @@ -1135,8 +1204,12 @@ namespace parallel reinit (c, true); } vector_view.reinit (partitioner->local_size(), val); - if (partitioner->local_size() > 0) + + if (partitioner->local_size()) vector_view.equ (1., c.vector_view); + + if (vector_is_ghosted || c.vector_is_ghosted) + update_ghost_values(); return *this; } @@ -1195,6 +1268,17 @@ namespace parallel std::fill_n (&val[partitioner->local_size()], partitioner->n_ghost_indices(), Number()); + vector_is_ghosted = false; + } + + + + template + inline + bool + Vector::has_ghost_elements () const + { + return vector_is_ghosted; } @@ -1360,7 +1444,7 @@ namespace parallel Vector::mean_value_local () const { Assert (partitioner->size()!=0, ExcEmptyObject()); - return (partitioner->local_size()>0 ? + return (partitioner->local_size() ? vector_view.mean_value() : Number()); } @@ -1389,7 +1473,7 @@ namespace parallel typename Vector::real_type Vector::l1_norm_local () const { - return partitioner->local_size()>0 ? vector_view.l1_norm() : real_type(); + return partitioner->local_size() ? vector_view.l1_norm() : real_type(); } @@ -1424,7 +1508,7 @@ namespace parallel typename Vector::real_type Vector::lp_norm_local (const real_type p) const { - return partitioner->local_size()>0 ? vector_view.lp_norm(p) : real_type(); + return partitioner->local_size() ? vector_view.lp_norm(p) : real_type(); } @@ -1450,7 +1534,7 @@ namespace parallel typename Vector::real_type Vector::linfty_norm_local () const { - return partitioner->local_size()>0 ? vector_view.linfty_norm() : real_type(); + return partitioner->local_size() ? vector_view.linfty_norm() : real_type(); } @@ -1519,8 +1603,7 @@ namespace parallel { IndexSet is (size()); - const std::pair x = local_range(); - is.add_range (x.first, x.second); + is.add_range (local_range().first, local_range().second); return is; } @@ -1602,6 +1685,10 @@ namespace parallel Number Vector::operator() (const size_type global_index) const { + // do not allow reading a vector which is not in ghost mode + Assert (in_local_range (global_index) || vector_is_ghosted == true, + ExcMessage("You tried to read a ghost element of this vector, " + "but it has not imported its ghost values.")); return val[partitioner->global_to_local(global_index)]; } @@ -1612,6 +1699,12 @@ namespace parallel Number & Vector::operator() (const size_type global_index) { + // we would like to prevent reading ghosts from a vector that does not + // have them imported, but this is not possible because we might be in a + // part of the code where the vector has enabled ghosts but is non-const + // (then, the compiler picks this method according to the C++ rule book + // even if a human would pick the const method when this subsequent use + // is just a read) return val[partitioner->global_to_local (global_index)]; } @@ -1656,10 +1749,11 @@ namespace parallel const ForwardIterator indices_end, OutputIterator values_begin) const { - while (indices_begin != indices_end) { - *values_begin = operator()(*indices_begin); - indices_begin++; values_begin++; - } + while (indices_begin != indices_end) + { + *values_begin = operator()(*indices_begin); + indices_begin++; values_begin++; + } } @@ -1672,6 +1766,10 @@ namespace parallel AssertIndexRange (local_index, partitioner->local_size()+ partitioner->n_ghost_indices()); + // do not allow reading a vector which is not in ghost mode + Assert (local_index < local_size() || vector_is_ghosted == true, + ExcMessage("You tried to read a ghost element of this vector, " + "but it has not imported its ghost values.")); return val[local_index]; } @@ -1695,8 +1793,8 @@ namespace parallel Vector & Vector::operator = (const Number s) { - // if we call Vector::operator=0, we want to - // zero out all the entries plus ghosts. + // if we call Vector::operator=0, we want to zero out all the entries + // plus ghosts. if (partitioner->local_size() > 0) vector_view.dealii::template Vector::operator= (s); if (s==Number()) @@ -1713,11 +1811,15 @@ namespace parallel Vector::operator += (const Vector &v) { AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation + + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation if (local_size()>0) vector_view += v.vector_view; + + if (vector_is_ghosted) + update_ghost_values(); + return *this; } @@ -1729,11 +1831,15 @@ namespace parallel Vector::operator -= (const Vector &v) { AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation + + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation if (local_size()>0) vector_view -= v.vector_view; + + if (vector_is_ghosted) + update_ghost_values(); + return *this; } @@ -1771,7 +1877,7 @@ namespace parallel void Vector::add (const size_type n_indices, const size_type *indices, - const OtherNumber *values) + const OtherNumber *values) { for (size_type i=0; i::add (const Number a) { - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.add (a); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1802,12 +1910,13 @@ namespace parallel void Vector::add (const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.add (v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1818,12 +1927,13 @@ namespace parallel Vector::add (const Number a, const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.add (a, v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1836,13 +1946,13 @@ namespace parallel const Number b, const Vector &w) { - AssertDimension (local_size(), v.local_size()); - AssertDimension (local_size(), w.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.add (a, v.vector_view, b, w.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1853,12 +1963,13 @@ namespace parallel Vector::sadd (const Number x, const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.sadd (x, v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1870,12 +1981,13 @@ namespace parallel const Number a, const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.sadd (x, a, v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1889,13 +2001,13 @@ namespace parallel const Number b, const Vector &w) { - AssertDimension (local_size(), v.local_size()); - AssertDimension (local_size(), w.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.sadd (x, a, v.vector_view, b, w.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1911,15 +2023,14 @@ namespace parallel const Number c, const Vector &x) { - AssertDimension (local_size(), v.local_size()); - AssertDimension (local_size(), w.local_size()); - AssertDimension (local_size(), x.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.sadd (s, a, v.vector_view, b, w.vector_view, c, x.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1929,11 +2040,7 @@ namespace parallel void Vector::scale (const Number factor) { - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) - vector_view *= factor; + operator *=(factor); } @@ -1943,11 +2050,14 @@ namespace parallel Vector & Vector::operator *= (const Number factor) { - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) - vector_view.operator *= (factor); + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) + vector_view *= factor; + + if (vector_is_ghosted) + update_ghost_values(); + return *this; } @@ -1958,11 +2068,7 @@ namespace parallel Vector & Vector::operator /= (const Number factor) { - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) - vector_view.operator /= (factor); + operator *= (1./factor); return *this; } @@ -1973,11 +2079,13 @@ namespace parallel void Vector::scale (const Vector &scaling_factors) { - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.scale (scaling_factors.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1988,7 +2096,11 @@ namespace parallel void Vector::scale (const Vector &scaling_factors) { - vector_view.template scale (scaling_factors.vector_view); + if (local_size()) + vector_view.template scale (scaling_factors.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -1999,12 +2111,13 @@ namespace parallel Vector::equ (const Number a, const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.equ (a, v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -2016,12 +2129,13 @@ namespace parallel Vector::equ (const Number a, const Vector &v) { - AssertDimension (local_size(), v.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.equ (a, v.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -2034,13 +2148,13 @@ namespace parallel const Number b, const Vector &w) { - AssertDimension (local_size(), v.local_size()); - AssertDimension (local_size(), w.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.equ (a, v.vector_view, b, w.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -2055,15 +2169,14 @@ namespace parallel const Number c, const Vector &x) { - AssertDimension (local_size(), v.local_size()); - AssertDimension (local_size(), w.local_size()); - AssertDimension (local_size(), w.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.equ (a, v.vector_view, b, w.vector_view, c, x.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -2074,13 +2187,13 @@ namespace parallel Vector::ratio (const Vector &a, const Vector &b) { - AssertDimension (local_size(), a.local_size()); - AssertDimension (local_size(), b.local_size()); - // dealii::Vector does not allow empty fields - // but this might happen on some processors - // for parallel implementation - if (local_size()>0) + // dealii::Vector does not allow empty fields but this might happen on + // some processors for parallel implementation + if (local_size()) vector_view.ratio (a.vector_view, b.vector_view); + + if (vector_is_ghosted) + update_ghost_values(); } diff --git a/deal.II/include/deal.II/lac/parallel_vector.templates.h b/deal.II/include/deal.II/lac/parallel_vector.templates.h index 334966c5f0..624168b2e7 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.templates.h +++ b/deal.II/include/deal.II/lac/parallel_vector.templates.h @@ -92,6 +92,8 @@ namespace parallel // set entries to zero if so requested if (fast == false) this->operator = (Number()); + + vector_is_ghosted = false; } @@ -134,6 +136,8 @@ namespace parallel // call these methods and hence do not need to have the storage. import_data = 0; } + + vector_is_ghosted = false; } @@ -194,6 +198,8 @@ namespace parallel // call these methods and hence do not need to have the storage. import_data = 0; } + + vector_is_ghosted = false; } @@ -209,6 +215,8 @@ namespace parallel vector_view = c.vector_view; if (call_update_ghost_values == true) update_ghost_values(); + else + vector_is_ghosted = false; } @@ -218,10 +226,12 @@ namespace parallel Vector::compress_start (const unsigned int counter, ::dealii::VectorOperation::values operation) { -#ifdef DEAL_II_WITH_MPI + Assert (vector_is_ghosted == false, + ExcMessage ("Cannot call compress() on a ghosted vector")); +#ifdef DEAL_II_WITH_MPI // nothing to do for insert (only need to zero ghost entries in - // compress_finish(). in debug mode we still want to check consistency + // compress_finish()). in debug mode we want to check consistency // of the inserted data, therefore the communication is still // initialized. Having different code in debug and optimized mode is // somewhat dangerous, but it really saves communication so it seems @@ -244,20 +254,17 @@ namespace parallel const size_type n_import_targets = part.import_targets().size(); const size_type n_ghost_targets = part.ghost_targets().size(); - // Need to send and receive the data. Use - // non-blocking communication, where it is - // generally less overhead to first initiate - // the receive and then actually send the data + // Need to send and receive the data. Use non-blocking communication, + // where it is generally less overhead to first initiate the receive and + // then actually send the data if (compress_requests.size() == 0) { - // set channels in different range from - // update_ghost_values channels + // set channels in different range from update_ghost_values channels const unsigned int channel = counter + 400; unsigned int current_index_start = 0; compress_requests.resize (n_import_targets + n_ghost_targets); - // allocate import_data in case it is not set - // up yet + // allocate import_data in case it is not set up yet if (import_data == 0) import_data = new Number[part.n_import_indices()]; for (size_type i=0; i 0) { - int ierr; - ierr = MPI_Startall(compress_requests.size(),&compress_requests[0]); + int ierr = MPI_Startall(compress_requests.size(),&compress_requests[0]); Assert (ierr == MPI_SUCCESS, ExcInternalError()); } @@ -305,7 +311,6 @@ namespace parallel (void)counter; (void)operation; #endif - } @@ -328,8 +333,7 @@ namespace parallel const Utilities::MPI::Partitioner &part = *partitioner; - // nothing to do when we neither have import - // nor ghost indices. + // nothing to do when we neither have import nor ghost indices. if (part.n_ghost_indices()==0 && part.n_import_indices()==0) return; @@ -346,9 +350,8 @@ namespace parallel // first wait for the receive to complete if (compress_requests.size() > 0 && n_import_targets > 0) { - int ierr; - ierr = MPI_Waitall (n_import_targets, &compress_requests[0], - MPI_STATUSES_IGNORE); + int ierr = MPI_Waitall (n_import_targets, &compress_requests[0], + MPI_STATUSES_IGNORE); Assert (ierr == MPI_SUCCESS, ExcInternalError()); Number *read_position = import_data; @@ -377,10 +380,9 @@ namespace parallel if (compress_requests.size() > 0 && n_ghost_targets > 0) { - int ierr; - ierr = MPI_Waitall (n_ghost_targets, - &compress_requests[n_import_targets], - MPI_STATUSES_IGNORE); + int ierr = MPI_Waitall (n_ghost_targets, + &compress_requests[n_import_targets], + MPI_STATUSES_IGNORE); Assert (ierr == MPI_SUCCESS, ExcInternalError()); } else @@ -401,8 +403,7 @@ namespace parallel #ifdef DEAL_II_WITH_MPI const Utilities::MPI::Partitioner &part = *partitioner; - // nothing to do when we neither have import - // nor ghost indices. + // nothing to do when we neither have import nor ghost indices. if (part.n_ghost_indices()==0 && part.n_import_indices()==0) return; @@ -412,10 +413,9 @@ namespace parallel const size_type n_import_targets = part.import_targets().size(); const size_type n_ghost_targets = part.ghost_targets().size(); - // Need to send and receive the data. Use - // non-blocking communication, where it is - // generally less overhead to first initiate - // the receive and then actually send the data + // Need to send and receive the data. Use non-blocking communication, + // where it is generally less overhead to first initiate the receive and + // then actually send the data if (update_ghost_values_requests.size() == 0) { Assert (part.local_size() == vector_view.size(), @@ -424,8 +424,8 @@ namespace parallel update_ghost_values_requests.resize (n_import_targets+n_ghost_targets); for (size_type i=0; i(&val[current_index_start]), part.ghost_targets()[i].second*sizeof(Number), MPI_BYTE, @@ -439,8 +439,7 @@ namespace parallel AssertDimension (current_index_start, part.local_size()+part.n_ghost_indices()); - // allocate import_data in case it is not set - // up yet + // allocate import_data in case it is not set up yet if (import_data == 0 && part.n_import_indices() > 0) import_data = new Number[part.n_import_indices()]; current_index_start = 0; @@ -458,8 +457,7 @@ namespace parallel AssertDimension (current_index_start, part.n_import_indices()); } - // copy the data that is actually to be send - // to the import_data field + // copy the data that is actually to be send to the import_data field if (part.n_import_indices() > 0) { Assert (import_data != 0, ExcInternalError()); @@ -475,9 +473,8 @@ namespace parallel update_ghost_values_requests.size()); if (update_ghost_values_requests.size() > 0) { - int ierr; - ierr = MPI_Startall(update_ghost_values_requests.size(), - &update_ghost_values_requests[0]); + int ierr = MPI_Startall(update_ghost_values_requests.size(), + &update_ghost_values_requests[0]); Assert (ierr == MPI_SUCCESS, ExcInternalError()); } #else @@ -492,10 +489,8 @@ namespace parallel Vector::update_ghost_values_finish () const { #ifdef DEAL_II_WITH_MPI - // wait for both sends and receives to - // complete, even though only receives are - // really necessary. this gives (much) better - // performance + // wait for both sends and receives to complete, even though only + // receives are really necessary. this gives (much) better performance AssertDimension (partitioner->ghost_targets().size() + partitioner->import_targets().size(), update_ghost_values_requests.size()); @@ -504,13 +499,13 @@ namespace parallel // make this function thread safe Threads::Mutex::ScopedLock lock (mutex); - int ierr; - ierr = MPI_Waitall (update_ghost_values_requests.size(), - &update_ghost_values_requests[0], - MPI_STATUSES_IGNORE); + int ierr = MPI_Waitall (update_ghost_values_requests.size(), + &update_ghost_values_requests[0], + MPI_STATUSES_IGNORE); Assert (ierr == MPI_SUCCESS, ExcInternalError()); } #endif + vector_is_ghosted = true; } @@ -520,28 +515,38 @@ namespace parallel Vector::swap (Vector &v) { #ifdef DEAL_II_WITH_MPI - // introduce a Barrier over all MPI processes - // to make sure that the compress request are - // no longer used before changing the owner - if (v.partitioner->n_mpi_processes() > 1) - MPI_Barrier (v.partitioner->get_communicator()); - if (partitioner->n_mpi_processes() > 1 && - v.partitioner->n_mpi_processes() != - partitioner->n_mpi_processes()) - MPI_Barrier (partitioner->get_communicator()); + +#ifdef DEBUG + // make sure that there are not outstanding requests from updating ghost + // values or compress + int flag = 1; + int ierr = MPI_Testall (update_ghost_values_requests.size(), + &update_ghost_values_requests[0], + &flag, MPI_STATUSES_IGNORE); + Assert (ierr == MPI_SUCCESS, ExcInternalError()); + Assert (flag == 1, + ExcMessage("MPI found unfinished update_ghost_values() requests" + "when calling swap, which is not allowed")); + ierr = MPI_Testall (compress_requests.size(), &compress_requests[0], + &flag, MPI_STATUSES_IGNORE); + Assert (ierr == MPI_SUCCESS, ExcInternalError()); + Assert (flag == 1, + ExcMessage("MPI found unfinished compress() requests " + "when calling swap, which is not allowed")); +#endif std::swap (compress_requests, v.compress_requests); std::swap (update_ghost_values_requests, v.update_ghost_values_requests); #endif - std::swap (partitioner, v.partitioner); - std::swap (allocated_size, v.allocated_size); - std::swap (val, v.val); - std::swap (import_data, v.import_data); + std::swap (partitioner, v.partitioner); + std::swap (allocated_size, v.allocated_size); + std::swap (val, v.val); + std::swap (import_data, v.import_data); + std::swap (vector_is_ghosted, v.vector_is_ghosted); - // vector view cannot be swapped so reset it - // manually (without touching the vector - // elements) + // vector view cannot be swapped so reset it manually (without touching + // the vector elements) vector_view.reinit (partitioner->local_size(), val); v.vector_view.reinit (v.partitioner->local_size(), v.val); } @@ -555,10 +560,9 @@ namespace parallel std::size_t memory = sizeof(*this); memory += sizeof (Number) * static_cast(allocated_size); - // if the partitioner is shared between more - // processors, just count a fraction of that - // memory, since we're not actually using more - // memory for it. + // if the partitioner is shared between more processors, just count a + // fraction of that memory, since we're not actually using more memory + // for it. if (partitioner.use_count() > 0) memory += partitioner->memory_consumption()/partitioner.use_count()+1; if (import_data != 0) @@ -587,10 +591,9 @@ namespace parallel else out.setf (std::ios::fixed, std::ios::floatfield); - // to make the vector write out all the - // information in order, use as many barriers - // as there are processors and start writing - // when it's our turn + // to make the vector write out all the information in order, use as + // many barriers as there are processors and start writing when it's our + // turn #ifdef DEAL_II_WITH_MPI if (partitioner->n_mpi_processes() > 1) for (unsigned int i=0; ithis_mpi_process(); i++) @@ -609,17 +612,22 @@ namespace parallel for (size_type i=0; ilocal_size(); ++i) out << local_element(i) << std::endl; out << std::endl; - out << "Ghost entries (global index / value):" << std::endl; - if (across) - for (size_type i=0; in_ghost_indices(); ++i) - out << '(' << partitioner->ghost_indices().nth_index_in_set(i) - << '/' << local_element(partitioner->local_size()+i) << ") "; - else - for (size_type i=0; in_ghost_indices(); ++i) - out << '(' << partitioner->ghost_indices().nth_index_in_set(i) - << '/' << local_element(partitioner->local_size()+i) << ")" - << std::endl; - out << std::endl << std::flush; + + if (vector_is_ghosted) + { + out << "Ghost entries (global index / value):" << std::endl; + if (across) + for (size_type i=0; in_ghost_indices(); ++i) + out << '(' << partitioner->ghost_indices().nth_index_in_set(i) + << '/' << local_element(partitioner->local_size()+i) << ") "; + else + for (size_type i=0; in_ghost_indices(); ++i) + out << '(' << partitioner->ghost_indices().nth_index_in_set(i) + << '/' << local_element(partitioner->local_size()+i) << ")" + << std::endl; + out << std::endl; + } + out << std::flush; #ifdef DEAL_II_WITH_MPI if (partitioner->n_mpi_processes() > 1) diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h index 8c8b0a83f6..59b9ae536f 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -1643,17 +1643,22 @@ reinit(const Mapping &mapping, // internal helper functions that define how to call MPI data exchange // functions: for generic vectors, do nothing at all. For distributed vectors, -// can call update_ghost_values_start function and so on. If we have -// collections of vectors, just do the individual functions of the -// components. this is a bit messy for block vectors, which use some -// additional helper functions to select the blocks +// call update_ghost_values_start function and so on. If we have collections +// of vectors, just do the individual functions of the components. In order to +// keep ghost values consistent (whether we are in read or write mode). the whole situation is a bit complicated by the fact +// that we need to treat block vectors differently, which use some additional +// helper functions to select the blocks and template magic. namespace internal { template - void update_ghost_values_start_block (const VectorStruct &vec, + bool update_ghost_values_start_block (const VectorStruct &vec, const unsigned int channel, internal::bool2type); template + void reset_ghost_values_block (const VectorStruct &vec, + const bool zero_out_ghosts, + internal::bool2type); + template void update_ghost_values_finish_block (const VectorStruct &vec, internal::bool2type); template @@ -1665,9 +1670,16 @@ namespace internal internal::bool2type); template - void update_ghost_values_start_block (const VectorStruct &, + bool update_ghost_values_start_block (const VectorStruct &, const unsigned int, internal::bool2type) + { + return false; + } + template + void reset_ghost_values_block (const VectorStruct &, + const bool, + internal::bool2type) {} template void update_ghost_values_finish_block (const VectorStruct &, @@ -1685,55 +1697,125 @@ namespace internal + // returns true if the vector was in a state without ghost values before, + // i.e., we need to zero out ghosts in the very end template inline - void update_ghost_values_start (const VectorStruct &vec, + bool update_ghost_values_start (const VectorStruct &vec, const unsigned int channel = 0) { - update_ghost_values_start_block(vec, channel, - internal::bool2type::value>()); + return + update_ghost_values_start_block(vec, channel, + internal::bool2type::value>()); } template inline - void update_ghost_values_start (const parallel::distributed::Vector &vec, + bool update_ghost_values_start (const parallel::distributed::Vector &vec, const unsigned int channel = 0) { + bool return_value = !vec.has_ghost_elements(); vec.update_ghost_values_start(channel); + return return_value; } template inline - void update_ghost_values_start (const std::vector &vec) + bool update_ghost_values_start (const std::vector &vec) { + bool return_value = false; for (unsigned int comp=0; comp inline - void update_ghost_values_start (const std::vector &vec) + bool update_ghost_values_start (const std::vector &vec) { + bool return_value = false; for (unsigned int comp=0; comp inline - void update_ghost_values_start_block (const VectorStruct &vec, + bool update_ghost_values_start_block (const VectorStruct &vec, const unsigned int channel, internal::bool2type) + { + bool return_value = false; + for (unsigned int i=0; i + inline + void reset_ghost_values (const VectorStruct &vec, + const bool zero_out_ghosts) + { + reset_ghost_values_block(vec, zero_out_ghosts, + internal::bool2type::value>()); + } + + + + template + inline + void reset_ghost_values (const parallel::distributed::Vector &vec, + const bool zero_out_ghosts) + { + if (zero_out_ghosts) + const_cast&>(vec).zero_out_ghosts(); + } + + + + template + inline + void reset_ghost_values (const std::vector &vec, + const bool zero_out_ghosts) + { + for (unsigned int comp=0; comp + inline + void reset_ghost_values (const std::vector &vec, + const bool zero_out_ghosts) + { + for (unsigned int comp=0; comp + inline + void reset_ghost_values_block (const VectorStruct &vec, + const bool zero_out_ghosts, + internal::bool2type) { for (unsigned int i=0; i::cell_loop OutVector &dst, const InVector &src) const { + // in any case, need to start the ghost import at the beginning + bool ghosts_were_not_set = internal::update_ghost_values_start (src); + #ifdef DEAL_II_WITH_THREADS // Use multithreading if so requested and if there is enough work to do in @@ -2181,7 +2266,6 @@ MatrixFree::cell_loop if (task_info.use_partition_partition == true) { - internal::update_ghost_values_start(src); tbb::empty_task *root = new( tbb::task::allocate_root() ) tbb::empty_task; unsigned int evens = task_info.evens; @@ -2252,11 +2336,9 @@ MatrixFree::cell_loop root->wait_for_all(); root->destroy(*root); - internal::compress_finish(dst); } else // end of partition-partition, start of partition-color { - internal::update_ghost_values_start(src); unsigned int evens = task_info.evens; unsigned int odds = task_info.odds; @@ -2387,7 +2469,6 @@ MatrixFree::cell_loop internal::compress_start(dst); } - internal::compress_finish(dst); } } else @@ -2396,8 +2477,6 @@ MatrixFree::cell_loop { std::pair cell_range; - internal::update_ghost_values_start (src); - // First operate on cells where no ghost data is needed (inner cells) { cell_range.first = 0; @@ -2426,9 +2505,11 @@ MatrixFree::cell_loop cell_range.second = size_info.n_macro_cells; cell_operation (*this, dst, src, cell_range); } - - internal::compress_finish(dst); } + + // In every case, we need to finish transfers at the very end + internal::compress_finish(dst); + internal::reset_ghost_values(src, ghosts_were_not_set); } diff --git a/tests/mpi/matrix_free_02.cc b/tests/mpi/matrix_free_02.cc index a9b26a0512..56fad9a8eb 100644 --- a/tests/mpi/matrix_free_02.cc +++ b/tests/mpi/matrix_free_02.cc @@ -188,8 +188,6 @@ void test () for (unsigned int parallel_option = 0; parallel_option < 3; ++parallel_option) { - out -= ref; - const double diff_norm = out.linfty_norm(); const QGauss<1> quad (fe_degree+1); typename MatrixFree::AdditionalData data; data.mpi_communicator = MPI_COMM_WORLD; diff --git a/tests/mpi/parallel_block_vector_01.cc b/tests/mpi/parallel_block_vector_01.cc index 6efa357d4b..e518efcd3b 100644 --- a/tests/mpi/parallel_block_vector_01.cc +++ b/tests/mpi/parallel_block_vector_01.cc @@ -47,7 +47,8 @@ void test () local_relevant = local_owned; local_relevant.add_range(1,2); - parallel::distributed::Vector v(local_owned, local_owned, MPI_COMM_WORLD); + parallel::distributed::Vector v(local_owned, local_relevant, + MPI_COMM_WORLD); // set local values if (myid < 8) diff --git a/tests/mpi/parallel_block_vector_02.cc b/tests/mpi/parallel_block_vector_02.cc new file mode 100644 index 0000000000..28d08f0df7 --- /dev/null +++ b/tests/mpi/parallel_block_vector_02.cc @@ -0,0 +1,162 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2011 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check ghost handling on parallel block vectors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + + +void test () +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + if (myid==0) deallog << "numproc=" << numproc << std::endl; + + + // each processor from processor 1 to 8 owns 2 indices (the other processors + // do not own any dof), and all processors are ghosting element 1 + IndexSet local_owned(std::min(16U,numproc*2)); + if (myid < 8) + local_owned.add_range(myid*2,myid*2+2); + IndexSet local_relevant(numproc*2); + local_relevant = local_owned; + local_relevant.add_range(1,2); + + parallel::distributed::Vector v(local_owned, local_relevant, + MPI_COMM_WORLD); + + // set local values + if (myid < 8) + { + v(myid*2)=myid*2.0; + v(myid*2+1)=myid*2.0+1.0; + } + v.compress(VectorOperation::insert); + + parallel::distributed::BlockVector w(3); + for (unsigned int i=0; i<3; ++i) + { + w.block(i) = v; + w.block(i) *= (i+1); + } + w.collect_sizes(); + + // see if we can access the ghosted entry 1 in each block with the correct + // value: the initialization should also update the ghost values, so all + // processors should have i+1 + for (unsigned int i=0; i<3; ++i) + AssertDimension(i+1,(unsigned int)w.block(i)(1)); + + // import ghost values, all processors should still have i+1 + w.update_ghost_values(); + for (unsigned int i=0; i<3; ++i) + AssertDimension(i+1,(unsigned int)w.block(i)(1)); + + // zero out ghosts, now all processors except processor 1 should have 0. + w.zero_out_ghosts(); + for (unsigned int i=0; i<3; ++i) + if (myid == 0) + { + AssertDimension(i+1,(unsigned int)w.block(i)(1)); + } + else + { + AssertDimension(0,(unsigned int)w.block(i)(1)); + } + + // create a vector copy that gets the entries from w. First, it should have + // updated the ghosts because it is created from an empty state. + parallel::distributed::BlockVector x(w); + Assert(x.has_ghost_elements() == true, ExcInternalError()); + for (unsigned int i=0; i<3; ++i) + AssertDimension(i+1,(unsigned int)x.block(i)(1)); + + // now we zero the vector, which should disable ghost elements + x = 0; + Assert(x.has_ghost_elements() == false, ExcInternalError()); + + // we copy from w (i.e., the same vector but one that does not have ghosts + // enabled) -> should not have ghosts enabled + x = w; + Assert(x.has_ghost_elements() == false, ExcInternalError()); + for (unsigned int i=0; i<3; ++i) + if (myid == 0) + { + AssertDimension(i+1,(unsigned int)x.block(i)(1)); + } + else + { + AssertDimension(0,(unsigned int)x.block(i)(1)); + } + + x.update_ghost_values(); + Assert(x.has_ghost_elements() == true, ExcInternalError()); + + + // add something to entry 1 on all processors + w(1) += myid+1; + w.compress(VectorOperation::add); + if (myid == 0) + AssertDimension((unsigned int)w(1), 1+(numproc*(numproc+1))/2); + + // add again and check if everything is still correct + w(1+v.size()) += myid+1; + w.compress(VectorOperation::add); + if (myid == 0) + AssertDimension((unsigned int)w(1), 1+(numproc*(numproc+1))/2); + if (myid == 0) + AssertDimension((unsigned int)w(v.size()+1), 2+(numproc*(numproc+1))/2); + + w.update_ghost_values(); + AssertDimension((unsigned int)w(1), 1+(numproc*(numproc+1))/2); + AssertDimension((unsigned int)w(v.size()+1), 2+(numproc*(numproc+1))/2); + + if (myid == 0) + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile(output_file_for_mpi("parallel_block_vector_02").c_str()); + deallog.attach(logfile); + deallog << std::setprecision(4); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test(); + } + else + test(); + +} diff --git a/tests/mpi/parallel_block_vector_02/ncpu_1/cmp/generic b/tests/mpi/parallel_block_vector_02/ncpu_1/cmp/generic new file mode 100644 index 0000000000..bb1593da60 --- /dev/null +++ b/tests/mpi/parallel_block_vector_02/ncpu_1/cmp/generic @@ -0,0 +1,3 @@ + +DEAL:0::numproc=1 +DEAL:0::OK diff --git a/tests/mpi/parallel_block_vector_02/ncpu_10/cmp/generic b/tests/mpi/parallel_block_vector_02/ncpu_10/cmp/generic new file mode 100644 index 0000000000..999baf1a98 --- /dev/null +++ b/tests/mpi/parallel_block_vector_02/ncpu_10/cmp/generic @@ -0,0 +1,3 @@ + +DEAL:0::numproc=10 +DEAL:0::OK diff --git a/tests/mpi/parallel_block_vector_02/ncpu_4/cmp/generic b/tests/mpi/parallel_block_vector_02/ncpu_4/cmp/generic new file mode 100644 index 0000000000..6f5d2775d9 --- /dev/null +++ b/tests/mpi/parallel_block_vector_02/ncpu_4/cmp/generic @@ -0,0 +1,3 @@ + +DEAL:0::numproc=4 +DEAL:0::OK -- 2.39.5