]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow multigrid transfer to use external vector partitioners
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 18 Sep 2019 20:45:46 +0000 (22:45 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 12 Nov 2019 11:07:26 +0000 (12:07 +0100)
include/deal.II/multigrid/mg_transfer_internal.h
include/deal.II/multigrid/mg_transfer_matrix_free.h
source/multigrid/mg_level_global_transfer.cc
source/multigrid/mg_transfer_internal.cc
source/multigrid/mg_transfer_internal.inst.in
source/multigrid/mg_transfer_matrix_free.cc

index e90e7ca3da9a83b01c2837a58d9b6432a08f0fc8..64e310ac23cedc1d79c88c40b8921d7de399e430 100644 (file)
@@ -120,8 +120,10 @@ namespace internal
     template <int dim, typename Number>
     void
     setup_transfer(
-      const DoFHandler<dim> &                 mg_dof,
-      const MGConstrainedDoFs *               mg_constrained_dofs,
+      const DoFHandler<dim> &  mg_dof,
+      const MGConstrainedDoFs *mg_constrained_dofs,
+      const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+        &                                     external_partitioners,
       ElementInfo<Number> &                   elem_info,
       std::vector<std::vector<unsigned int>> &level_dof_indices,
       std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
@@ -130,8 +132,8 @@ namespace internal
       std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
       std::vector<std::vector<Number>> &                     weights_on_refined,
       std::vector<Table<2, unsigned int>> &copy_indices_global_mine,
-      MGLevelObject<LinearAlgebra::distributed::Vector<Number>>
-        &ghosted_level_vector);
+      MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>>
+        &vector_partitioners);
 
   } // namespace MGTransfer
 } // namespace internal
index 111ca464bf1de5d861ca8da2da802959a081eac6..a3d02ee30c05d5dde36edc594985326b28c98772 100644 (file)
@@ -89,9 +89,24 @@ public:
 
   /**
    * Actually build the information for the prolongation for each level.
+   *
+   * The optional second argument of external partitioners allows the user to
+   * suggest vector partitioning on the levels. In case the partitioners
+   * are found to contain all ghost unknowns that are visited through the
+   * transfer, the given partitioners are chosen. This ensures compatibility
+   * of vectors during prolongate and restrict with external partitioners as
+   * given by the user, which in turn saves some copy operations. However, in
+   * case there are unknowns missing -- and this is typically the case at some
+   * point during h-coarsening since processors will need to drop out and
+   * thus children's unknowns on some processor will be needed as ghosts to a
+   * parent cell on another processor -- the provided external partitioners are
+   * ignored and internal variants are used instead.
    */
   void
-  build(const DoFHandler<dim, dim> &mg_dof);
+  build(const DoFHandler<dim, dim> &mg_dof,
+        const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+          &external_partitioners =
+            std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>());
 
   /**
    * Prolongate a vector from level <tt>to_level-1</tt> to level
@@ -245,6 +260,15 @@ private:
    */
   std::vector<std::vector<std::vector<unsigned short>>> dirichlet_indices;
 
+  /**
+   * A vector that holds shared pointers to the partitioners of the
+   * transfer. These partitioners might be shared with what was passed in from
+   * the outside through build() or be shared with the level vectors inherited
+   * from MGLevelGlobalTransfer.
+   */
+  MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>>
+    vector_partitioners;
+
   /**
    * Perform the prolongation operation.
    */
index 633b4ff98011c51da36760aae3d7547401437511..964f1a9d13efd4bb74c643eb287eb8531cd29989 100644 (file)
@@ -300,15 +300,39 @@ MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>::
                 ghosted_global_vector,
                 ghosted_level_vector);
 
-  fill_internal(mg_dof,
-                mg_constrained_dofs,
-                mpi_communicator,
-                true,
-                this->solution_copy_indices,
-                this->solution_copy_indices_global_mine,
-                this->solution_copy_indices_level_mine,
-                solution_ghosted_global_vector,
-                solution_ghosted_level_vector);
+  // in case we have hanging nodes which imply different ownership of
+  // global and level dof indices, we must also fill the solution indices
+  // which contain additional indices, otherwise they can re-use the
+  // indices of the standard transfer.
+  int have_refinement_edge_dofs = 0;
+  if (mg_constrained_dofs != nullptr)
+    for (unsigned int level = 0;
+         level < mg_dof.get_triangulation().n_global_levels();
+         ++level)
+      if (mg_constrained_dofs->get_refinement_edge_indices(level).n_elements() >
+          0)
+        {
+          have_refinement_edge_dofs = 1;
+          break;
+        }
+  if (Utilities::MPI::max(have_refinement_edge_dofs, mpi_communicator) == 1)
+    fill_internal(mg_dof,
+                  mg_constrained_dofs,
+                  mpi_communicator,
+                  true,
+                  this->solution_copy_indices,
+                  this->solution_copy_indices_global_mine,
+                  this->solution_copy_indices_level_mine,
+                  solution_ghosted_global_vector,
+                  solution_ghosted_level_vector);
+  else
+    {
+      this->solution_copy_indices             = this->copy_indices;
+      this->solution_copy_indices_global_mine = this->copy_indices_global_mine;
+      this->solution_copy_indices_level_mine  = this->copy_indices_level_mine;
+      solution_ghosted_global_vector          = ghosted_global_vector;
+      solution_ghosted_level_vector           = ghosted_level_vector;
+    }
 
   bool my_perform_renumbered_plain_copy =
     (this->copy_indices.back().n_cols() ==
index a11aa5b66cbec27f040384084b25a5925b6970a8..57135a56990e921d829eae95461256ef23073717 100644 (file)
@@ -222,7 +222,8 @@ namespace internal
           "We should only be sending information with a parallel Triangulation!"));
 
 #ifdef DEAL_II_WITH_MPI
-      if (tria)
+      if (tria && Utilities::MPI::sum(send_data_temp.size(),
+                                      tria->get_communicator()) > 0)
         {
           const std::set<types::subdomain_id> &neighbors =
             tria->level_ghost_owners();
@@ -419,14 +420,15 @@ namespace internal
 
     // initialize the vectors needed for the transfer (and merge with the
     // content in copy_indices_global_mine)
-    template <typename Number>
     void
-    reinit_ghosted_vector(
-      const IndexSet &                            locally_owned,
-      std::vector<types::global_dof_index> &      ghosted_level_dofs,
-      const MPI_Comm &                            communicator,
-      LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
-      Table<2, unsigned int> &                    copy_indices_global_mine)
+    reinit_level_partitioner(
+      const IndexSet &                      locally_owned,
+      std::vector<types::global_dof_index> &ghosted_level_dofs,
+      const std::shared_ptr<const Utilities::MPI::Partitioner>
+        &                                                 external_partitioner,
+      const MPI_Comm &                                    communicator,
+      std::shared_ptr<const Utilities::MPI::Partitioner> &target_partitioner,
+      Table<2, unsigned int> &copy_indices_global_mine)
     {
       std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
       IndexSet ghosted_dofs(locally_owned.size());
@@ -436,21 +438,52 @@ namespace internal
       ghosted_dofs.compress();
 
       // Add possible ghosts from the previous content in the vector
-      if (ghosted_level_vector.size() == locally_owned.size())
+      if (target_partitioner.get() != nullptr &&
+          target_partitioner->size() == locally_owned.size())
+        {
+          ghosted_dofs.add_indices(target_partitioner->ghost_indices());
+        }
+
+      // check if the given partitioner's ghosts represent a superset of the
+      // ghosts we require in this function
+      const int ghosts_locally_contained =
+        (external_partitioner.get() != nullptr &&
+         (external_partitioner->ghost_indices() & ghosted_dofs) ==
+           ghosted_dofs) ?
+          1 :
+          0;
+      if (external_partitioner.get() != nullptr &&
+          Utilities::MPI::min(ghosts_locally_contained, communicator) == 1)
         {
           // shift the local number of the copy indices according to the new
-          // partitioner that we are going to use for the vector
-          const auto &part = ghosted_level_vector.get_partitioner();
-          ghosted_dofs.add_indices(part->ghost_indices());
-          for (unsigned int i = 0; i < copy_indices_global_mine.n_cols(); ++i)
-            copy_indices_global_mine(1, i) =
-              locally_owned.n_elements() +
-              ghosted_dofs.index_within_set(
-                part->local_to_global(copy_indices_global_mine(1, i)));
+          // partitioner that we are going to use during the access to the
+          // entries
+          if (target_partitioner.get() != nullptr &&
+              target_partitioner->size() == locally_owned.size())
+            for (unsigned int i = 0; i < copy_indices_global_mine.n_cols(); ++i)
+              copy_indices_global_mine(1, i) =
+                external_partitioner->global_to_local(
+                  target_partitioner->local_to_global(
+                    copy_indices_global_mine(1, i)));
+          target_partitioner = external_partitioner;
+        }
+      else
+        {
+          if (target_partitioner.get() != nullptr &&
+              target_partitioner->size() == locally_owned.size())
+            for (unsigned int i = 0; i < copy_indices_global_mine.n_cols(); ++i)
+              copy_indices_global_mine(1, i) =
+                locally_owned.n_elements() +
+                ghosted_dofs.index_within_set(
+                  target_partitioner->local_to_global(
+                    copy_indices_global_mine(1, i)));
+          target_partitioner.reset(new Utilities::MPI::Partitioner(
+            locally_owned, ghosted_dofs, communicator));
         }
-      ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator);
     }
 
+
+
     // Transform the ghost indices to local index space for the vector
     inline void
     copy_indices_to_mpi_local_numbers(
@@ -470,6 +503,8 @@ namespace internal
           localized_indices[i + mine.size()] = part.global_to_local(remote[i]);
     }
 
+
+
     // given the collection of child cells in lexicographic ordering as seen
     // from the parent, compute the first index of the given child
     template <int dim>
@@ -498,6 +533,8 @@ namespace internal
       return shift;
     }
 
+
+
     // puts the indices on the given child cell in lexicographic ordering with
     // respect to the collection of all child cells as seen from the parent
     template <int dim>
@@ -535,6 +572,8 @@ namespace internal
               }
     }
 
+
+
     template <int dim, typename Number>
     void
     setup_element_info(ElementInfo<Number> &          elem_info,
@@ -596,6 +635,8 @@ namespace internal
               fe.get_prolongation_matrix(c)(renumbering[j], renumbering[i]);
     }
 
+
+
     namespace
     {
       /**
@@ -627,13 +668,17 @@ namespace internal
       }
     } // namespace
 
+
+
     // Sets up most of the internal data structures of the MGTransferMatrixFree
     // class
     template <int dim, typename Number>
     void
     setup_transfer(
-      const dealii::DoFHandler<dim> &         mg_dof,
-      const MGConstrainedDoFs *               mg_constrained_dofs,
+      const dealii::DoFHandler<dim> &mg_dof,
+      const MGConstrainedDoFs *      mg_constrained_dofs,
+      const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+        &                                     external_partitioners,
       ElementInfo<Number> &                   elem_info,
       std::vector<std::vector<unsigned int>> &level_dof_indices,
       std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
@@ -642,8 +687,8 @@ namespace internal
       std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
       std::vector<std::vector<Number>> &                     weights_on_refined,
       std::vector<Table<2, unsigned int>> &copy_indices_global_mine,
-      MGLevelObject<LinearAlgebra::distributed::Vector<Number>>
-        &ghosted_level_vector)
+      MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>>
+        &target_partitioners)
     {
       level_dof_indices.clear();
       parent_child_connect.clear();
@@ -691,11 +736,10 @@ namespace internal
         mg_dof.get_fe().dofs_per_cell);
       dirichlet_indices.resize(n_levels - 1);
 
-      // We use the vectors stored ghosted_level_vector in the base class for
-      // keeping ghosted transfer indices. To avoid keeping two very similar
-      // vectors, we merge them here.
-      if (ghosted_level_vector.max_level() != n_levels - 1)
-        ghosted_level_vector.resize(0, n_levels - 1);
+      AssertDimension(target_partitioners.max_level(), n_levels - 1);
+      Assert(external_partitioners.empty() ||
+               external_partitioners.size() == n_levels,
+             ExcDimensionMismatch(external_partitioners.size(), n_levels));
 
       for (unsigned int level = n_levels - 1; level > 0; --level)
         {
@@ -892,38 +936,48 @@ namespace internal
                   i->first += counter;
                 }
 
-          // step 2.7: Initialize the ghosted vector
+          // step 2.7: Initialize the partitioner for the ghosted vector
+          //
+          // We use a vector based on the target partitioner handed in also in
+          // the base class for keeping ghosted transfer indices. To avoid
+          // keeping two very similar vectors, we keep one single ghosted
+          // vector that is augmented/filled here.
           const parallel::TriangulationBase<dim, dim> *ptria =
             (dynamic_cast<const parallel::TriangulationBase<dim, dim> *>(
               &tria));
           const MPI_Comm communicator =
             ptria != nullptr ? ptria->get_communicator() : MPI_COMM_SELF;
 
-          reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level),
-                                ghosted_level_dofs,
-                                communicator,
-                                ghosted_level_vector[level],
-                                copy_indices_global_mine[level]);
-
-          copy_indices_to_mpi_local_numbers(
-            *ghosted_level_vector[level].get_partitioner(),
-            global_level_dof_indices,
-            global_level_dof_indices_remote,
-            level_dof_indices[level]);
+          reinit_level_partitioner(mg_dof.locally_owned_mg_dofs(level),
+                                   ghosted_level_dofs,
+                                   external_partitioners.empty() ?
+                                     nullptr :
+                                     external_partitioners[level],
+                                   communicator,
+                                   target_partitioners[level],
+                                   copy_indices_global_mine[level]);
+
+          copy_indices_to_mpi_local_numbers(*target_partitioners[level],
+                                            global_level_dof_indices,
+                                            global_level_dof_indices_remote,
+                                            level_dof_indices[level]);
           // step 2.8: Initialize the ghosted vector for level 0
           if (level == 1)
             {
               for (unsigned int i = 0; i < parent_child_connect[0].size(); ++i)
                 parent_child_connect[0][i] = std::make_pair(i, 0U);
 
-              reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
-                                    ghosted_level_dofs_l0,
-                                    communicator,
-                                    ghosted_level_vector[0],
-                                    copy_indices_global_mine[0]);
+              reinit_level_partitioner(mg_dof.locally_owned_mg_dofs(0),
+                                       ghosted_level_dofs_l0,
+                                       external_partitioners.empty() ?
+                                         nullptr :
+                                         external_partitioners[0],
+                                       communicator,
+                                       target_partitioners[0],
+                                       copy_indices_global_mine[0]);
 
               copy_indices_to_mpi_local_numbers(
-                *ghosted_level_vector[0].get_partitioner(),
+                *target_partitioners[0],
                 global_level_dof_indices_l0,
                 std::vector<types::global_dof_index>(),
                 level_dof_indices[0]);
@@ -940,14 +994,15 @@ namespace internal
       weights_on_refined.resize(n_levels - 1);
       for (unsigned int level = 1; level < n_levels; ++level)
         {
-          ghosted_level_vector[level] = 0;
+          LinearAlgebra::distributed::Vector<Number> touch_count(
+            target_partitioners[level]);
           for (unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
             for (unsigned int j = 0; j < elem_info.n_child_cell_dofs; ++j)
-              ghosted_level_vector[level].local_element(
+              touch_count.local_element(
                 level_dof_indices[level][elem_info.n_child_cell_dofs * c +
                                          j]) += Number(1.);
-          ghosted_level_vector[level].compress(VectorOperation::add);
-          ghosted_level_vector[level].update_ghost_values();
+          touch_count.compress(VectorOperation::add);
+          touch_count.update_ghost_values();
 
           std::vector<unsigned int> degree_to_3(n_child_dofs_1d);
           degree_to_3[0] = 0;
@@ -970,7 +1025,7 @@ namespace internal
                                        1][c * Utilities::fixed_power<dim>(3) +
                                           shift + degree_to_3[i]] =
                       Number(1.) /
-                      ghosted_level_vector[level].local_element(
+                      touch_count.local_element(
                         level_dof_indices[level]
                                          [elem_info.n_child_cell_dofs * c + m]);
                 }
index 204503303a92896c23901c674e83521a60ec64cb..0c8841afbff8e9c7ac66cca5f00f47408efeb757 100644 (file)
@@ -79,6 +79,8 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS)
         setup_transfer<deal_II_dimension>(
           const dealii::DoFHandler<deal_II_dimension> &,
           const MGConstrainedDoFs *,
+          const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+            &,
           ElementInfo<S> &,
           std::vector<std::vector<unsigned int>> &,
           std::vector<std::vector<std::pair<unsigned int, unsigned int>>> &,
@@ -86,7 +88,7 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS)
           std::vector<std::vector<std::vector<unsigned short>>> &,
           std::vector<std::vector<S>> &,
           std::vector<Table<2, unsigned int>> &,
-          MGLevelObject<LinearAlgebra::distributed::Vector<S>> &);
+          MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>> &);
       \}
     \}
   }
index 769c532e8c60f7c179c1ea6c742f83a7a5a94201..e7e087e2b13cc0160ab2d8371e71304db27264c9 100644 (file)
@@ -94,10 +94,20 @@ MGTransferMatrixFree<dim, Number>::clear()
 
 template <int dim, typename Number>
 void
-MGTransferMatrixFree<dim, Number>::build(const DoFHandler<dim, dim> &mg_dof)
+MGTransferMatrixFree<dim, Number>::build(
+  const DoFHandler<dim, dim> &mg_dof,
+  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+    &external_partitioners)
 {
   this->fill_and_communicate_copy_indices(mg_dof);
 
+  vector_partitioners.resize(0,
+                             mg_dof.get_triangulation().n_global_levels() - 1);
+  for (unsigned int level = 0; level <= this->ghosted_level_vector.max_level();
+       ++level)
+    vector_partitioners[level] =
+      this->ghosted_level_vector[level].get_partitioner();
+
   std::vector<std::vector<Number>> weights_unvectorized;
 
   internal::MGTransfer::ElementInfo<Number> elem_info;
@@ -105,6 +115,7 @@ MGTransferMatrixFree<dim, Number>::build(const DoFHandler<dim, dim> &mg_dof)
   internal::MGTransfer::setup_transfer<dim, Number>(
     mg_dof,
     this->mg_constrained_dofs,
+    external_partitioners,
     elem_info,
     level_dof_indices,
     parent_child_connect,
@@ -112,7 +123,19 @@ MGTransferMatrixFree<dim, Number>::build(const DoFHandler<dim, dim> &mg_dof)
     dirichlet_indices,
     weights_unvectorized,
     this->copy_indices_global_mine,
-    this->ghosted_level_vector);
+    vector_partitioners);
+
+  // reinit the ghosted level vector in case its partitioner differs from the
+  // one the user intends to use externally
+  this->ghosted_level_vector.resize(0, vector_partitioners.max_level());
+  for (unsigned int level = 0; level <= vector_partitioners.max_level();
+       ++level)
+    if (external_partitioners.size() == vector_partitioners.max_level() + 1 &&
+        external_partitioners[level].get() == vector_partitioners[level].get())
+      this->ghosted_level_vector[level].reinit(0);
+    else
+      this->ghosted_level_vector[level].reinit(vector_partitioners[level]);
+
   // unpack element info data
   fe_degree             = elem_info.fe_degree;
   element_is_continuous = elem_info.element_is_continuous;
@@ -163,69 +186,73 @@ MGTransferMatrixFree<dim, Number>::prolongate(
   Assert((to_level >= 1) && (to_level <= level_dof_indices.size()),
          ExcIndexRange(to_level, 1, level_dof_indices.size() + 1));
 
-  AssertDimension(this->ghosted_level_vector[to_level].local_size(),
-                  dst.local_size());
-  AssertDimension(this->ghosted_level_vector[to_level - 1].local_size(),
-                  src.local_size());
+  const bool src_inplace = src.get_partitioner().get() ==
+                           this->vector_partitioners[to_level - 1].get();
+  if (src_inplace == false)
+    {
+      if (this->ghosted_level_vector[to_level - 1].get_partitioner().get() !=
+          this->vector_partitioners[to_level - 1].get())
+        this->ghosted_level_vector[to_level - 1].reinit(
+          this->vector_partitioners[to_level - 1]);
+      this->ghosted_level_vector[to_level - 1].copy_locally_owned_data_from(
+        src);
+    }
+
+  const bool dst_inplace =
+    dst.get_partitioner().get() == this->vector_partitioners[to_level].get();
+  if (dst_inplace == false)
+    {
+      if (this->ghosted_level_vector[to_level].get_partitioner().get() !=
+          this->vector_partitioners[to_level].get())
+        this->ghosted_level_vector[to_level].reinit(
+          this->vector_partitioners[to_level]);
+      AssertDimension(this->ghosted_level_vector[to_level].local_size(),
+                      dst.local_size());
+      this->ghosted_level_vector[to_level] = 0.;
+    }
+  else
+    dst = 0;
 
-  this->ghosted_level_vector[to_level - 1].copy_locally_owned_data_from(src);
-  this->ghosted_level_vector[to_level - 1].update_ghost_values();
-  this->ghosted_level_vector[to_level] = 0.;
+  const LinearAlgebra::distributed::Vector<Number> &src_vec =
+    src_inplace ? src : this->ghosted_level_vector[to_level - 1];
+  LinearAlgebra::distributed::Vector<Number> &dst_vec =
+    dst_inplace ? dst : this->ghosted_level_vector[to_level];
 
+  src_vec.update_ghost_values();
   // the implementation in do_prolongate_add is templated in the degree of the
   // element (for efficiency reasons), so we need to find the appropriate
   // kernel here...
   if (fe_degree == 0)
-    do_prolongate_add<0>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<0>(to_level, dst_vec, src_vec);
   else if (fe_degree == 1)
-    do_prolongate_add<1>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<1>(to_level, dst_vec, src_vec);
   else if (fe_degree == 2)
-    do_prolongate_add<2>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<2>(to_level, dst_vec, src_vec);
   else if (fe_degree == 3)
-    do_prolongate_add<3>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<3>(to_level, dst_vec, src_vec);
   else if (fe_degree == 4)
-    do_prolongate_add<4>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<4>(to_level, dst_vec, src_vec);
   else if (fe_degree == 5)
-    do_prolongate_add<5>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<5>(to_level, dst_vec, src_vec);
   else if (fe_degree == 6)
-    do_prolongate_add<6>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<6>(to_level, dst_vec, src_vec);
   else if (fe_degree == 7)
-    do_prolongate_add<7>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<7>(to_level, dst_vec, src_vec);
   else if (fe_degree == 8)
-    do_prolongate_add<8>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<8>(to_level, dst_vec, src_vec);
   else if (fe_degree == 9)
-    do_prolongate_add<9>(to_level,
-                         this->ghosted_level_vector[to_level],
-                         this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<9>(to_level, dst_vec, src_vec);
   else if (fe_degree == 10)
-    do_prolongate_add<10>(to_level,
-                          this->ghosted_level_vector[to_level],
-                          this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<10>(to_level, dst_vec, src_vec);
   else
-    do_prolongate_add<-1>(to_level,
-                          this->ghosted_level_vector[to_level],
-                          this->ghosted_level_vector[to_level - 1]);
+    do_prolongate_add<-1>(to_level, dst_vec, src_vec);
 
-  this->ghosted_level_vector[to_level].compress(VectorOperation::add);
-  dst.copy_locally_owned_data_from(this->ghosted_level_vector[to_level]);
+  dst_vec.compress(VectorOperation::add);
+  if (dst_inplace == false)
+    dst = dst_vec;
+
+  if (src_inplace == true)
+    src.zero_out_ghosts();
 }
 
 
@@ -240,67 +267,69 @@ MGTransferMatrixFree<dim, Number>::restrict_and_add(
   Assert((from_level >= 1) && (from_level <= level_dof_indices.size()),
          ExcIndexRange(from_level, 1, level_dof_indices.size() + 1));
 
-  AssertDimension(this->ghosted_level_vector[from_level].local_size(),
-                  src.local_size());
-  AssertDimension(this->ghosted_level_vector[from_level - 1].local_size(),
-                  dst.local_size());
+  const bool src_inplace =
+    src.get_partitioner().get() == this->vector_partitioners[from_level].get();
+  if (src_inplace == false)
+    {
+      if (this->ghosted_level_vector[from_level].get_partitioner().get() !=
+          this->vector_partitioners[from_level].get())
+        this->ghosted_level_vector[from_level].reinit(
+          this->vector_partitioners[from_level]);
+      this->ghosted_level_vector[from_level].copy_locally_owned_data_from(src);
+    }
+
+  const bool dst_inplace = dst.get_partitioner().get() ==
+                           this->vector_partitioners[from_level - 1].get();
+  if (dst_inplace == false)
+    {
+      if (this->ghosted_level_vector[from_level - 1].get_partitioner().get() !=
+          this->vector_partitioners[from_level - 1].get())
+        this->ghosted_level_vector[from_level - 1].reinit(
+          this->vector_partitioners[from_level - 1]);
+      AssertDimension(this->ghosted_level_vector[from_level - 1].local_size(),
+                      dst.local_size());
+      this->ghosted_level_vector[from_level - 1] = 0.;
+    }
 
-  this->ghosted_level_vector[from_level].copy_locally_owned_data_from(src);
-  this->ghosted_level_vector[from_level].update_ghost_values();
-  this->ghosted_level_vector[from_level - 1] = 0.;
+  const LinearAlgebra::distributed::Vector<Number> &src_vec =
+    src_inplace ? src : this->ghosted_level_vector[from_level];
+  LinearAlgebra::distributed::Vector<Number> &dst_vec =
+    dst_inplace ? dst : this->ghosted_level_vector[from_level - 1];
+
+  src_vec.update_ghost_values();
 
   if (fe_degree == 0)
-    do_restrict_add<0>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<0>(from_level, dst_vec, src_vec);
   else if (fe_degree == 1)
-    do_restrict_add<1>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<1>(from_level, dst_vec, src_vec);
   else if (fe_degree == 2)
-    do_restrict_add<2>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<2>(from_level, dst_vec, src_vec);
   else if (fe_degree == 3)
-    do_restrict_add<3>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<3>(from_level, dst_vec, src_vec);
   else if (fe_degree == 4)
-    do_restrict_add<4>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<4>(from_level, dst_vec, src_vec);
   else if (fe_degree == 5)
-    do_restrict_add<5>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<5>(from_level, dst_vec, src_vec);
   else if (fe_degree == 6)
-    do_restrict_add<6>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<6>(from_level, dst_vec, src_vec);
   else if (fe_degree == 7)
-    do_restrict_add<7>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<7>(from_level, dst_vec, src_vec);
   else if (fe_degree == 8)
-    do_restrict_add<8>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<8>(from_level, dst_vec, src_vec);
   else if (fe_degree == 9)
-    do_restrict_add<9>(from_level,
-                       this->ghosted_level_vector[from_level - 1],
-                       this->ghosted_level_vector[from_level]);
+    do_restrict_add<9>(from_level, dst_vec, src_vec);
   else if (fe_degree == 10)
-    do_restrict_add<10>(from_level,
-                        this->ghosted_level_vector[from_level - 1],
-                        this->ghosted_level_vector[from_level]);
+    do_restrict_add<10>(from_level, dst_vec, src_vec);
   else
     // go to the non-templated version of the evaluator
-    do_restrict_add<-1>(from_level,
-                        this->ghosted_level_vector[from_level - 1],
-                        this->ghosted_level_vector[from_level]);
+    do_restrict_add<-1>(from_level, dst_vec, src_vec);
+
+  dst_vec.compress(VectorOperation::add);
+  if (dst_inplace == false)
+    dst += dst_vec;
 
-  this->ghosted_level_vector[from_level - 1].compress(VectorOperation::add);
-  dst += this->ghosted_level_vector[from_level - 1];
+  if (src_inplace == true)
+    src.zero_out_ghosts();
 }
 
 

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.