</p>
<ol>
+ <li>
+ 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 <i>same</i> vector can
+ also be used as a non-ghosted vector which is designed for use in assembly
+ routines.
+ <br>
+ (Martin Kronbichler, 2013/10/18)
+ </li>
+
<li>
Removed: GridTools::collect_periodic_face_pairs. This function is superseded
by GridTools::collect_periodic_faces which exports an
<h3>Specific improvements</h3>
<ol>
+ <li>
+ 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.
+ <br>
+ (Martin Kronbichler, 2013/10/18)
+ </li>
+
<li>
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).
<br>
*/
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
void
BlockVector<Number>::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; block<this->n_blocks(); ++block)
- this->block(block).compress(operation);
+ this->block(block).compress_start(block*10 + 8273, operation);
+ for (unsigned int block=0; block<this->n_blocks(); ++block)
+ this->block(block).compress_finish(operation);
}
BlockVector<Number>::update_ghost_values () const
{
for (unsigned int block=0; block<this->n_blocks(); ++block)
- this->block(block).update_ghost_values();
+ this->block(block).update_ghost_values_start(block*10 + 9923);
+ for (unsigned int block=0; block<this->n_blocks(); ++block)
+ this->block(block).update_ghost_values_finish();
+ }
+
+
+
+ template <typename Number>
+ inline
+ void
+ BlockVector<Number>::zero_out_ghosts ()
+ {
+ for (unsigned int block=0; block<this->n_blocks(); ++block)
+ this->block(block).zero_out_ghosts();
+ }
+
+
+
+ template <typename Number>
+ inline
+ bool
+ BlockVector<Number>::has_ghost_elements () const
+ {
+ bool has_ghost_elements = false;
+ for (unsigned int block=0; block<this->n_blocks(); ++block)
+ if (this->block(block).has_ghost_elements() == true)
+ has_ghost_elements = true;
+ return has_ghost_elements;
}
* - 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
+ * <code>operator = (0.)</code>), 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
+ * <code>true</code> exactly when ghost elements have been updated and
+ * <code>false</code> 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 <typename Number>
* 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;
* 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.
*/
/**
* 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
*/
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.
allocated_size (0),
val (0),
import_data (0),
+ vector_is_ghosted (false),
vector_view (0, static_cast<Number *>(0))
{}
allocated_size (0),
val (0),
import_data (0),
+ vector_is_ghosted (false),
vector_view (0, static_cast<Number *>(0))
{
reinit (v, true);
allocated_size (0),
val (0),
import_data (0),
+ vector_is_ghosted (false),
vector_view (0, static_cast<Number *>(0))
{
reinit (local_range, ghost_indices, communicator);
template <typename Number>
inline
Vector<Number>::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<Number *>(0))
{
IndexSet ghost_indices(local_range.size());
allocated_size (0),
val (0),
import_data (0),
+ vector_is_ghosted (false),
vector_view (0, static_cast<Number *>(0))
{
reinit (size, false);
allocated_size (0),
val (0),
import_data (0),
+ vector_is_ghosted (false),
vector_view (0, static_cast<Number *>(0))
{
reinit (partitioner);
{
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())
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;
}
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;
}
std::fill_n (&val[partitioner->local_size()],
partitioner->n_ghost_indices(),
Number());
+ vector_is_ghosted = false;
+ }
+
+
+
+ template <typename Number>
+ inline
+ bool
+ Vector<Number>::has_ghost_elements () const
+ {
+ return vector_is_ghosted;
}
Vector<Number>::mean_value_local () const
{
Assert (partitioner->size()!=0, ExcEmptyObject());
- return (partitioner->local_size()>0 ?
+ return (partitioner->local_size() ?
vector_view.mean_value()
: Number());
}
typename Vector<Number>::real_type
Vector<Number>::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();
}
typename Vector<Number>::real_type
Vector<Number>::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();
}
typename Vector<Number>::real_type
Vector<Number>::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();
}
{
IndexSet is (size());
- const std::pair<types::global_dof_index,types::global_dof_index> x = local_range();
- is.add_range (x.first, x.second);
+ is.add_range (local_range().first, local_range().second);
return is;
}
Number
Vector<Number>::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)];
}
Number &
Vector<Number>::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)];
}
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++;
+ }
}
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];
}
Vector<Number> &
Vector<Number>::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<Number>::operator= (s);
if (s==Number())
Vector<Number>::operator += (const Vector<Number> &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;
}
Vector<Number>::operator -= (const Vector<Number> &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;
}
void
Vector<Number>::add (const size_type n_indices,
const size_type *indices,
- const OtherNumber *values)
+ const OtherNumber *values)
{
for (size_type i=0; i<n_indices; ++i)
{
void
Vector<Number>::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();
}
void
Vector<Number>::add (const Vector<Number> &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();
}
Vector<Number>::add (const Number a,
const Vector<Number> &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();
}
const Number b,
const Vector<Number> &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();
}
Vector<Number>::sadd (const Number x,
const Vector<Number> &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();
}
const Number a,
const Vector<Number> &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();
}
const Number b,
const Vector<Number> &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();
}
const Number c,
const Vector<Number> &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();
}
void
Vector<Number>::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);
}
Vector<Number> &
Vector<Number>::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;
}
Vector<Number> &
Vector<Number>::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;
}
void
Vector<Number>::scale (const Vector<Number> &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();
}
void
Vector<Number>::scale (const Vector<Number2> &scaling_factors)
{
- vector_view.template scale<Number2> (scaling_factors.vector_view);
+ if (local_size())
+ vector_view.template scale<Number2> (scaling_factors.vector_view);
+
+ if (vector_is_ghosted)
+ update_ghost_values();
}
Vector<Number>::equ (const Number a,
const Vector<Number> &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();
}
Vector<Number>::equ (const Number a,
const Vector<Number2> &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();
}
const Number b,
const Vector<Number> &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();
}
const Number c,
const Vector<Number> &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();
}
Vector<Number>::ratio (const Vector<Number> &a,
const Vector<Number> &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();
}
// set entries to zero if so requested
if (fast == false)
this->operator = (Number());
+
+ vector_is_ghosted = false;
}
// call these methods and hence do not need to have the storage.
import_data = 0;
}
+
+ vector_is_ghosted = false;
}
// call these methods and hence do not need to have the storage.
import_data = 0;
}
+
+ vector_is_ghosted = false;
}
vector_view = c.vector_view;
if (call_update_ghost_values == true)
update_ghost_values();
+ else
+ vector_is_ghosted = false;
}
Vector<Number>::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
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<n_import_targets; i++)
compress_requests.size());
if (compress_requests.size() > 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());
}
(void)counter;
(void)operation;
#endif
-
}
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;
// 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;
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
#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;
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(),
update_ghost_values_requests.resize (n_import_targets+n_ghost_targets);
for (size_type i=0; i<n_ghost_targets; i++)
{
- // allow writing into ghost indices even
- // though we are in a const function
+ // allow writing into ghost indices even though we are in a
+ // const function
MPI_Recv_init (const_cast<Number *>(&val[current_index_start]),
part.ghost_targets()[i].second*sizeof(Number),
MPI_BYTE,
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;
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());
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
Vector<Number>::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());
// 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;
}
Vector<Number>::swap (Vector<Number> &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);
}
std::size_t memory = sizeof(*this);
memory += sizeof (Number) * static_cast<std::size_t>(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)
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; i<partitioner->this_mpi_process(); i++)
for (size_type i=0; i<partitioner->local_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; i<partitioner->n_ghost_indices(); ++i)
- out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
- << '/' << local_element(partitioner->local_size()+i) << ") ";
- else
- for (size_type i=0; i<partitioner->n_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; i<partitioner->n_ghost_indices(); ++i)
+ out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
+ << '/' << local_element(partitioner->local_size()+i) << ") ";
+ else
+ for (size_type i=0; i<partitioner->n_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)
// 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<typename VectorStruct>
- void update_ghost_values_start_block (const VectorStruct &vec,
+ bool update_ghost_values_start_block (const VectorStruct &vec,
const unsigned int channel,
internal::bool2type<true>);
template<typename VectorStruct>
+ void reset_ghost_values_block (const VectorStruct &vec,
+ const bool zero_out_ghosts,
+ internal::bool2type<true>);
+ template<typename VectorStruct>
void update_ghost_values_finish_block (const VectorStruct &vec,
internal::bool2type<true>);
template<typename VectorStruct>
internal::bool2type<true>);
template<typename VectorStruct>
- void update_ghost_values_start_block (const VectorStruct &,
+ bool update_ghost_values_start_block (const VectorStruct &,
const unsigned int,
internal::bool2type<false>)
+ {
+ return false;
+ }
+ template<typename VectorStruct>
+ void reset_ghost_values_block (const VectorStruct &,
+ const bool,
+ internal::bool2type<false>)
{}
template<typename VectorStruct>
void update_ghost_values_finish_block (const VectorStruct &,
+ // 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<typename VectorStruct>
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<IsBlockVector<VectorStruct>::value>());
+ return
+ update_ghost_values_start_block(vec, channel,
+ internal::bool2type<IsBlockVector<VectorStruct>::value>());
}
template<typename Number>
inline
- void update_ghost_values_start (const parallel::distributed::Vector<Number> &vec,
+ bool update_ghost_values_start (const parallel::distributed::Vector<Number> &vec,
const unsigned int channel = 0)
{
+ bool return_value = !vec.has_ghost_elements();
vec.update_ghost_values_start(channel);
+ return return_value;
}
template <typename VectorStruct>
inline
- void update_ghost_values_start (const std::vector<VectorStruct> &vec)
+ bool update_ghost_values_start (const std::vector<VectorStruct> &vec)
{
+ bool return_value = false;
for (unsigned int comp=0; comp<vec.size(); comp++)
- update_ghost_values_start(vec[comp], comp);
+ return_value = update_ghost_values_start(vec[comp], comp);
+ return return_value;
}
template <typename VectorStruct>
inline
- void update_ghost_values_start (const std::vector<VectorStruct *> &vec)
+ bool update_ghost_values_start (const std::vector<VectorStruct *> &vec)
{
+ bool return_value = false;
for (unsigned int comp=0; comp<vec.size(); comp++)
- update_ghost_values_start(*vec[comp], comp);
+ return_value = update_ghost_values_start(*vec[comp], comp);
+ return return_value;
}
template<typename VectorStruct>
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<true>)
+ {
+ bool return_value = false;
+ for (unsigned int i=0; i<vec.n_blocks(); ++i)
+ return_value = update_ghost_values_start(vec.block(i), channel+509*i);
+ return return_value;
+ }
+
+
+
+ // if the input vector did not have ghosts imported, clear them here again
+ // in order to avoid subsequent operations e.g. in linear solvers to work
+ // with ghosts all the time
+ template<typename VectorStruct>
+ inline
+ void reset_ghost_values (const VectorStruct &vec,
+ const bool zero_out_ghosts)
+ {
+ reset_ghost_values_block(vec, zero_out_ghosts,
+ internal::bool2type<IsBlockVector<VectorStruct>::value>());
+ }
+
+
+
+ template<typename Number>
+ inline
+ void reset_ghost_values (const parallel::distributed::Vector<Number> &vec,
+ const bool zero_out_ghosts)
+ {
+ if (zero_out_ghosts)
+ const_cast<parallel::distributed::Vector<Number>&>(vec).zero_out_ghosts();
+ }
+
+
+
+ template <typename VectorStruct>
+ inline
+ void reset_ghost_values (const std::vector<VectorStruct> &vec,
+ const bool zero_out_ghosts)
+ {
+ for (unsigned int comp=0; comp<vec.size(); comp++)
+ reset_ghost_values(vec[comp], zero_out_ghosts);
+ }
+
+
+
+ template <typename VectorStruct>
+ inline
+ void reset_ghost_values (const std::vector<VectorStruct *> &vec,
+ const bool zero_out_ghosts)
+ {
+ for (unsigned int comp=0; comp<vec.size(); comp++)
+ reset_ghost_values(*vec[comp], zero_out_ghosts);
+ }
+
+
+
+ template<typename VectorStruct>
+ inline
+ void reset_ghost_values_block (const VectorStruct &vec,
+ const bool zero_out_ghosts,
+ internal::bool2type<true>)
{
for (unsigned int i=0; i<vec.n_blocks(); ++i)
- update_ghost_values_start(vec.block(i), channel+500*i);
+ reset_ghost_values(vec.block(i), zero_out_ghosts);
}
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
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;
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;
internal::compress_start(dst);
}
- internal::compress_finish(dst);
}
}
else
{
std::pair<unsigned int,unsigned int> cell_range;
- internal::update_ghost_values_start (src);
-
// First operate on cells where no ghost data is needed (inner cells)
{
cell_range.first = 0;
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);
}
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<dim,number>::AdditionalData data;
data.mpi_communicator = MPI_COMM_WORLD;
local_relevant = local_owned;
local_relevant.add_range(1,2);
- parallel::distributed::Vector<double> v(local_owned, local_owned, MPI_COMM_WORLD);
+ parallel::distributed::Vector<double> v(local_owned, local_relevant,
+ MPI_COMM_WORLD);
// set local values
if (myid < 8)
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/parallel_block_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+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<double> 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<double> 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<double> 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();
+
+}
--- /dev/null
+
+DEAL:0::numproc=1
+DEAL:0::OK
--- /dev/null
+
+DEAL:0::numproc=10
+DEAL:0::OK
--- /dev/null
+
+DEAL:0::numproc=4
+DEAL:0::OK
projects / dealii.git / commitdiff