From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 26 Oct 2000 15:48:17 +0000 (+0000)
Subject: Speed up the compute_integrid_constraints function by some parallelization.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b52a1c71af96dfb8da5b38881d3e427cacf9b94c;p=dealii-svn.git

Speed up the compute_integrid_constraints function by some parallelization.


git-svn-id: https://svn.dealii.org/trunk@3461 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h
index 1786321975..5b227216b9 100644
--- a/deal.II/deal.II/include/dofs/dof_tools.h
+++ b/deal.II/deal.II/include/dofs/dof_tools.h
@@ -870,6 +870,44 @@ class DoFTools
 				      * Exception
 				      */
     DeclException0 (ExcInvalidBoundaryIndicator);
+
+  private:
+				     /**
+				      * This is a helper function that
+				      * is used in the computation of
+				      * integrid constraints. See the
+				      * function for a thorough
+				      * description of how it works.
+				      */
+    template <int dim>
+    static void
+    compute_intergrid_weights (const DoFHandler<dim>              &coarse_grid,
+			       const unsigned int                  coarse_component,
+			       const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+			       const vector<Vector<double> >      &parameter_dofs,
+			       const vector<int>                  &weight_mapping,
+			       FullMatrix<double>                 &weights);
+
+				     /**
+				      * This is a function that is
+				      * called by the previous one and
+				      * that operates on a range of
+				      * cells only. It is used to
+				      * split up the whole range of
+				      * cells into chunks which are
+				      * then worked on in parallel, if
+				      * multithreading is available.
+				      */
+    template <int dim>
+    static void
+    compute_intergrid_weights_1 (const DoFHandler<dim>              &coarse_grid,
+				 const unsigned int                  coarse_component,
+				 const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+				 const vector<Vector<double> >      &parameter_dofs,
+				 const vector<int>                  &weight_mapping,
+				 FullMatrix<double>                 &weights,
+				 const typename DoFHandler<dim>::active_cell_iterator &begin,
+				 const typename DoFHandler<dim>::active_cell_iterator &end);
 };
 
 
diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index bc5f965742..9c6aa0b71d 100644
--- a/deal.II/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/deal.II/source/dofs/dof_tools.cc
@@ -12,6 +12,8 @@
 //----------------------------  dof_tools.cc  ---------------------------
 
 
+#include <base/multithread_info.h>
+#include <base/thread_management.h>
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <grid/intergrid_map.h>
@@ -1155,11 +1157,6 @@ DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coa
 					    true);
     };
 
-				   // vector to hold the representation of
-				   // a single degree of freedom on the
-				   // coarse grid (for the selected fe)
-				   // on the fine grid
-  Vector<double> global_parameter_representation (n_fine_dofs);
 
 				   // store the weights with which a dof
 				   // on the parameter grid contributes
@@ -1219,6 +1216,201 @@ DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coa
     };
   
   
+				   // for each cell on the parameter grid:
+				   // find out which degrees of freedom on the
+				   // fine grid correspond in which way to
+				   // the degrees of freedom on the parameter
+				   // grid
+				   //
+				   // do this in a separate function
+				   // to allow for multithreading
+				   // there. see this function also if
+				   // you want to read more
+				   // information on the algorithm
+				   // used.
+  compute_intergrid_weights (coarse_grid, coarse_component,
+			     coarse_to_fine_grid_map, parameter_dofs,
+			     weight_mapping, weights);
+
+  
+				   // ok, now we have all weights for each
+				   // dof on the fine grid. if in debug
+				   // mode lets see if everything went smooth,
+				   // i.e. each dof has sum of weights one
+				   //
+				   // in other words this means that
+				   // if the sum of all shape
+				   // functions on the parameter grid
+				   // is one (which is always the
+				   // case), then the representation
+				   // on the state grid should be as
+				   // well (division of unity)
+				   //
+				   // if the parameter grid has more
+				   // than one component, then the
+				   // respective dofs of the other
+				   // components have sum of weights
+				   // zero, of course. we do not
+				   // explicitely ask which component
+				   // a dof belongs to, but this at
+				   // least tests some errors
+#ifdef DEBUG
+  for (unsigned int col=0; col<weights.n(); ++col)
+    {
+      double sum=0;
+      for (unsigned int row=0; row<weights.m(); ++row)
+	sum += weights(row,col);
+      Assert ((sum==1) ||
+	      ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
+    };
+#endif
+
+
+				   // now we know that the weights in
+				   // each row constitute a
+				   // constraint. enter this into the
+				   // constraints object
+				   //
+				   // first task: for each parameter
+				   // dof on the parameter grid, find
+				   // a representant on the fine,
+				   // global grid. this is possible
+				   // since we use conforming finite
+				   // element. we take this
+				   // representant to be the first
+				   // element in this row with weight
+				   // identical to one. the
+				   // representant will become an
+				   // unconstrained degree of freedom,
+				   // while all others will be
+				   // constrained to this dof (and
+				   // possibly others)
+  vector<int> representants(weights.m(), -1);
+  for (unsigned int parameter_dof=0; parameter_dof<weights.m(); ++parameter_dof)
+    {
+      unsigned int column=0;
+      for (; column<weights.n(); ++column)
+	if (weights(parameter_dof,column) == 1)
+	  break;
+      Assert (column < weights.n(), ExcInternalError());
+
+				       // now we know in which column of
+				       // weights the representant is, but
+				       // we don't know its global index. get
+				       // it using the inverse operation of
+				       // the weight_mapping
+      unsigned int global_dof=0;
+      for (; global_dof<weight_mapping.size(); ++global_dof)
+	if (weight_mapping[global_dof] == static_cast<int>(column))
+	  break;
+      Assert (global_dof < weight_mapping.size(), ExcInternalError());
+
+				       // now enter the representants global
+				       // index into our list
+      representants[parameter_dof] = global_dof;
+    };
+
+
+				   // note for people that want to
+				   // optimize this function: the
+				   // largest part of the computing
+				   // time is spent in the following,
+				   // rather innocent block of
+				   // code. basically, it must be the
+				   // ConstraintMatrix::add_entry call
+				   // which takes the bulk of the
+				   // time, but it is not known to the
+				   // author how to make it faster...
+  vector<pair<unsigned int,double> > constraint_line;
+  for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
+    if (weight_mapping[global_dof] != -1)
+				       // this global dof is a parameter
+				       // dof, so it may carry a constraint
+				       // note that for each global dof,
+				       // the sum of weights shall be one,
+				       // so we can find out whether this
+				       // dof is constrained in the following
+				       // way: if the only weight in this row
+				       // is a one, and the representant for
+				       // the parameter dof of the line in
+				       // which this one is is the present
+				       // dof, then we consider this dof
+				       // to be unconstrained. otherwise,
+				       // all other dofs are constrained
+      {
+	unsigned int first_used_row=0;
+	for (; first_used_row<weights.m(); ++first_used_row)
+	  if (weights(first_used_row,weight_mapping[global_dof]) != 0)
+	    break;
+
+	if ((weights(first_used_row,weight_mapping[global_dof]) == 1) &&
+	    (representants[first_used_row] == static_cast<int>(global_dof)))
+					   // dof unconstrained
+	  continue;
+
+					 // otherwise enter all constraints
+	constraints.add_line (global_dof);
+
+	constraint_line.clear ();
+	for (unsigned int row=first_used_row; row<weights.m(); ++row)
+	  if (weights(row,weight_mapping[global_dof]) != 0)
+	    constraint_line.push_back (make_pair(representants[row],
+						 weights(row,weight_mapping[global_dof])));
+
+	constraints.add_entries (global_dof, constraint_line);
+      };
+};
+
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_weights (const DoFHandler<dim>              &coarse_grid,
+				     const unsigned int                  coarse_component,
+				     const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+				     const vector<Vector<double> >      &parameter_dofs,
+				     const vector<int>                  &weight_mapping,
+				     FullMatrix<double>                 &weights)
+{
+				   // simply distribute the range of
+				   // cells to different threads
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+  vector<pair<active_cell_iterator,active_cell_iterator> >
+    cell_intervals = Threads::split_range<active_cell_iterator> (coarse_grid.begin_active(),
+								 coarse_grid.end(),
+								 multithread_info.n_default_threads);
+  
+  Threads::ThreadManager thread_manager;
+  for (unsigned int i=0; i<multithread_info.n_default_threads; ++i)
+    Threads::spawn (thread_manager,
+		    Threads::encapsulate (&DoFTools::template compute_intergrid_weights_1<dim>)
+		    .collect_args (coarse_grid, coarse_component,
+				   coarse_to_fine_grid_map, parameter_dofs,
+				   weight_mapping, weights,
+				   cell_intervals[i].first,
+				   cell_intervals[i].second));
+
+				   // wait for the threads to finish
+  thread_manager.wait ();
+};
+
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_weights_1 (const DoFHandler<dim>              &coarse_grid,
+				       const unsigned int                  coarse_component,
+				       const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+				       const vector<Vector<double> >      &parameter_dofs,
+				       const vector<int>                  &weight_mapping,
+				       FullMatrix<double>                 &weights,
+				       const typename DoFHandler<dim>::active_cell_iterator &begin,
+				       const typename DoFHandler<dim>::active_cell_iterator &end)
+{
+				   // aliases to the finite elements
+				   // used by the dof handlers:
+  const FiniteElement<dim> &coarse_fe = coarse_grid.get_fe();    
+  
 				   // for each cell on the parameter grid:
 				   // find out which degrees of freedom on the
 				   // fine grid correspond in which way to
@@ -1283,10 +1475,18 @@ DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coa
 				   // both cells 1 and 2, but the
 				   // correct weight is nevertheless
 				   // only 1.
+
+				   // vector to hold the representation of
+				   // a single degree of freedom on the
+				   // coarse grid (for the selected fe)
+				   // on the fine grid
+  const unsigned int n_fine_dofs = weight_mapping.size();
+  Vector<double> global_parameter_representation (n_fine_dofs);
+  
   typename DoFHandler<dim>::active_cell_iterator cell;
   vector<unsigned int> parameter_dof_indices (coarse_fe.dofs_per_cell);
   
-  for (cell=coarse_grid.begin_active(); cell!=coarse_grid.end(); ++cell)
+  for (cell=begin; cell!=end; ++cell)
     {
 				       // get the global indices of the
 				       // parameter dofs on this parameter
@@ -1358,118 +1558,6 @@ DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coa
 			ExcInternalError());
 	  };
     };
-  
-				   // ok, now we have all weights for each
-				   // dof on the fine grid. if in debug
-				   // mode lets see if everything went smooth,
-				   // i.e. each dof has sum of weights one
-				   //
-				   // in other words this means that
-				   // if the sum of all shape
-				   // functions on the parameter grid
-				   // is one (which is always the
-				   // case), then the representation
-				   // on the state grid should be as
-				   // well (division of unity)
-				   //
-				   // if the parameter grid has more
-				   // than one component, then the
-				   // respective dofs of the other
-				   // components have sum of weights
-				   // zero, of course. we do not
-				   // explicitely ask which component
-				   // a dof belongs to, but this at
-				   // least tests some errors
-#ifdef DEBUG
-  for (unsigned int col=0; col<weights.n(); ++col)
-    {
-      double sum=0;
-      for (unsigned int row=0; row<weights.m(); ++row)
-	sum += weights(row,col);
-      Assert ((sum==1) ||
-	      ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
-    };
-#endif
-
-
-				   // now we know that the weights in
-				   // each row constitute a
-				   // constraint. enter this into the
-				   // constraints object
-				   //
-				   // first task: for each parameter
-				   // dof on the parameter grid, find
-				   // a representant on the fine,
-				   // global grid. this is possible
-				   // since we use conforming finite
-				   // element. we take this
-				   // representant to be the first
-				   // element in this row with weight
-				   // identical to one. the
-				   // representant will become an
-				   // unconstrained degree of freedom,
-				   // while all others will be
-				   // constrained to this dof (and
-				   // possibly others)
-  vector<int> representants(weights.m(), -1);
-  for (unsigned int parameter_dof=0; parameter_dof<weights.m(); ++parameter_dof)
-    {
-      unsigned int column=0;
-      for (; column<weights.n(); ++column)
-	if (weights(parameter_dof,column) == 1)
-	  break;
-      Assert (column < weights.n(), ExcInternalError());
-
-				       // now we know in which column of
-				       // weights the representant is, but
-				       // we don't know its global index. get
-				       // it using the inverse operation of
-				       // the weight_mapping
-      unsigned int global_dof=0;
-      for (; global_dof<weight_mapping.size(); ++global_dof)
-	if (weight_mapping[global_dof] == static_cast<int>(column))
-	  break;
-      Assert (global_dof < weight_mapping.size(), ExcInternalError());
-
-				       // now enter the representants global
-				       // index into our list
-      representants[parameter_dof] = global_dof;
-    };
-  
-  for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
-    if (weight_mapping[global_dof] != -1)
-				       // this global dof is a parameter
-				       // dof, so it may carry a constraint
-				       // note that for each global dof,
-				       // the sum of weights shall be one,
-				       // so we can find out whether this
-				       // dof is constrained in the following
-				       // way: if the only weight in this row
-				       // is a one, and the representant for
-				       // the parameter dof of the line in
-				       // which this one is is the present
-				       // dof, then we consider this dof
-				       // to be unconstrained. otherwise,
-				       // all other dofs are constrained
-      {
-	unsigned int first_used_row=0;
-	for (; first_used_row<weights.m(); ++first_used_row)
-	  if (weights(first_used_row,weight_mapping[global_dof]) != 0)
-	    break;
-
-	if ((weights(first_used_row,weight_mapping[global_dof]) == 1) &&
-	    (representants[first_used_row] == static_cast<int>(global_dof)))
-					   // dof unconstrained
-	  continue;
-
-					 // otherwise enter all constraints
-	constraints.add_line (global_dof);
-	for (unsigned int row=first_used_row; row<weights.m(); ++row)
-	  if (weights(row,weight_mapping[global_dof]) != 0)
-	    constraints.add_entry (global_dof,
-				   representants[row],
-				   weights(row,weight_mapping[global_dof]));
-      };
 };
 
 
@@ -1628,13 +1716,13 @@ template
 void
 DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
 					 const Vector<float>  &cell_data,
-					 Vector<double>       &dof_data) const;
+					 Vector<double>       &dof_data);
 
 template
 void
 DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
 					 const Vector<double> &cell_data,
-					 Vector<double>       &dof_data) const;
+					 Vector<double>       &dof_data);
 
 
 template void DoFTools::extract_dofs(const DoFHandler<deal_II_dimension>& dof,