]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Improve handling of ghost elements in parallel vector: Distinguish a state with read...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Oct 2013 12:19:18 +0000 (12:19 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Oct 2013 12:19:18 +0000 (12:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@31301 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/parallel_block_vector.h
deal.II/include/deal.II/lac/parallel_vector.h
deal.II/include/deal.II/lac/parallel_vector.templates.h
deal.II/include/deal.II/matrix_free/matrix_free.h
tests/mpi/matrix_free_02.cc
tests/mpi/parallel_block_vector_01.cc
tests/mpi/parallel_block_vector_02.cc [new file with mode: 0644]
tests/mpi/parallel_block_vector_02/ncpu_1/cmp/generic [new file with mode: 0644]
tests/mpi/parallel_block_vector_02/ncpu_10/cmp/generic [new file with mode: 0644]
tests/mpi/parallel_block_vector_02/ncpu_4/cmp/generic [new file with mode: 0644]

index 54d0a2412bf08d629557c48f4ba06edb117c3cbf..c399798fed1b1d39196cbf42372ed7eedaefb702 100644 (file)
@@ -24,6 +24,23 @@ inconvenience this causes.
 </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
@@ -145,6 +162,14 @@ inconvenience this causes.
 <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>
index 3be362a6359e73908308ccb266f51928be8ffe3e..7e3f1bb1905bd343f1050b8aab06e45aeea336e4 100644 (file)
@@ -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<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);
     }
 
 
@@ -596,7 +615,34 @@ namespace parallel
     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;
     }
 
 
index 7e1d853e52f5b31128ce007133c9b546fea25d7b..c568a4aba7a804f9bef4ea65baab170b2372ea8e 100644 (file)
@@ -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
+     *   <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>
@@ -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<Number *>(0))
     {}
 
@@ -990,6 +1047,7 @@ namespace parallel
       allocated_size (0),
       val (0),
       import_data (0),
+      vector_is_ghosted (false),
       vector_view (0, static_cast<Number *>(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<Number *>(0))
     {
       reinit (local_range, ghost_indices, communicator);
@@ -1017,11 +1076,12 @@ namespace parallel
     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());
@@ -1037,6 +1097,7 @@ namespace parallel
       allocated_size (0),
       val (0),
       import_data (0),
+      vector_is_ghosted (false),
       vector_view (0, static_cast<Number *>(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<Number *>(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 <typename Number>
+    inline
+    bool
+    Vector<Number>::has_ghost_elements () const
+    {
+      return vector_is_ghosted;
     }
 
 
@@ -1360,7 +1444,7 @@ namespace parallel
     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());
     }
@@ -1389,7 +1473,7 @@ namespace parallel
     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();
     }
 
 
@@ -1424,7 +1508,7 @@ namespace parallel
     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();
     }
 
 
@@ -1450,7 +1534,7 @@ namespace parallel
     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();
     }
 
 
@@ -1519,8 +1603,7 @@ namespace parallel
     {
       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;
     }
@@ -1602,6 +1685,10 @@ namespace parallel
     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)];
     }
 
@@ -1612,6 +1699,12 @@ namespace parallel
     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)];
     }
 
@@ -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<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())
@@ -1713,11 +1811,15 @@ namespace parallel
     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;
     }
 
@@ -1729,11 +1831,15 @@ namespace parallel
     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;
     }
 
@@ -1771,7 +1877,7 @@ namespace parallel
     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)
         {
@@ -1788,11 +1894,13 @@ namespace parallel
     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();
     }
 
 
@@ -1802,12 +1910,13 @@ namespace parallel
     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();
     }
 
 
@@ -1818,12 +1927,13 @@ namespace parallel
     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();
     }
 
 
@@ -1836,13 +1946,13 @@ namespace parallel
                          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();
     }
 
 
@@ -1853,12 +1963,13 @@ namespace parallel
     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();
     }
 
 
@@ -1870,12 +1981,13 @@ namespace parallel
                           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();
     }
 
 
@@ -1889,13 +2001,13 @@ namespace parallel
                           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();
     }
 
 
@@ -1911,15 +2023,14 @@ namespace parallel
                           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();
     }
 
 
@@ -1929,11 +2040,7 @@ namespace parallel
     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);
     }
 
 
@@ -1943,11 +2050,14 @@ namespace parallel
     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;
     }
 
@@ -1958,11 +2068,7 @@ namespace parallel
     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;
     }
 
@@ -1973,11 +2079,13 @@ namespace parallel
     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();
     }
 
 
@@ -1988,7 +2096,11 @@ namespace parallel
     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();
     }
 
 
@@ -1999,12 +2111,13 @@ namespace parallel
     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();
     }
 
 
@@ -2016,12 +2129,13 @@ namespace parallel
     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();
     }
 
 
@@ -2034,13 +2148,13 @@ namespace parallel
                          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();
     }
 
 
@@ -2055,15 +2169,14 @@ namespace parallel
                          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();
     }
 
 
@@ -2074,13 +2187,13 @@ namespace parallel
     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();
     }
 
 
index 334966c5f099f57c4c842d5bc689c229ee984886..624168b2e757e497fe45a7e9c43e31b1bca8af54 100644 (file)
@@ -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<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
@@ -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<n_import_targets; i++)
@@ -296,8 +303,7 @@ namespace parallel
                       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());
         }
 
@@ -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<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,
@@ -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<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());
@@ -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<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);
     }
@@ -555,10 +560,9 @@ namespace parallel
       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)
@@ -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; i<partitioner->this_mpi_process(); i++)
@@ -609,17 +612,22 @@ namespace parallel
         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)
index 8c8b0a83f66385ecdb23e319b1d66a1541fbcc1e..59b9ae536f57139f7b664fbf66c55d7d6e4d987c 100644 (file)
@@ -1643,17 +1643,22 @@ reinit(const Mapping<dim>                         &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<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>
@@ -1665,9 +1670,16 @@ namespace internal
                               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 &,
@@ -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<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);
   }
 
 
@@ -2161,6 +2243,9 @@ MatrixFree<dim, Number>::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<dim, Number>::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<dim, Number>::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<dim, Number>::cell_loop
 
               internal::compress_start(dst);
             }
-          internal::compress_finish(dst);
         }
     }
   else
@@ -2396,8 +2477,6 @@ MatrixFree<dim, Number>::cell_loop
     {
       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;
@@ -2426,9 +2505,11 @@ MatrixFree<dim, Number>::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);
 }
 
 
index a9b26a05124031d9d9002d121fb3eae165be7a3a..56fad9a8eb9e32111af99a12d60612b3be532733 100644 (file)
@@ -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<dim,number>::AdditionalData data;
       data.mpi_communicator = MPI_COMM_WORLD;
index 6efa357d4b460f371275f99235d13f371813bd29..e518efcd3b68e3c4112c958d1c49a7b8111956d8 100644 (file)
@@ -47,7 +47,8 @@ void test ()
   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)
diff --git a/tests/mpi/parallel_block_vector_02.cc b/tests/mpi/parallel_block_vector_02.cc
new file mode 100644 (file)
index 0000000..28d08f0
--- /dev/null
@@ -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 <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();
+
+}
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 (file)
index 0000000..bb1593d
--- /dev/null
@@ -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 (file)
index 0000000..999baf1
--- /dev/null
@@ -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 (file)
index 0000000..6f5d277
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0::numproc=4
+DEAL:0::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.