From 7ad517777070ba4716c30a5437f4784b8aa1cc7e Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 6 Mar 2008 21:38:12 +0000
Subject: [PATCH] Interface to the minimum degree algorithm in boost to reorder
 matrices.

git-svn-id: https://svn.dealii.org/trunk@15865 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/dofs/dof_renumbering.h    |  36 +++++
 .../deal.II/source/dofs/dof_renumbering.cc    | 130 ++++++++++++++++++
 2 files changed, 166 insertions(+)

diff --git a/deal.II/deal.II/include/dofs/dof_renumbering.h b/deal.II/deal.II/include/dofs/dof_renumbering.h
index 8ba22cb130..b72a0344b8 100644
--- a/deal.II/deal.II/include/dofs/dof_renumbering.h
+++ b/deal.II/deal.II/include/dofs/dof_renumbering.h
@@ -298,6 +298,42 @@ class DoFRenumbering
 			       const DH&,
 			       const bool reversed_numbering = false,
 			       const bool use_constraints    = false);
+
+					 /**
+					  * Renumber the degrees of
+					  * freedom based on the BOOST
+					  * implementation of the
+					  * minimum degree
+					  * algorithm. This often
+					  * results in slightly larger
+					  * (by a few percent)
+					  * bandwidths than the
+					  * Cuthill-McKee algorithm,
+					  * but sparse ILUs are often
+					  * slightly (also by a few
+					  * percent) better
+					  * preconditioners.
+					  */
+	template <class DH>
+	static void 
+	minimum_degree (DH&                              dof_handler,
+			const bool                       reversed_numbering = false,
+			const bool                       use_constraints    = false);
+
+					 /**
+					  * Compute the renumbering
+					  * for the minimum degree
+					  * algorithm but do not
+					  * actually renumber the
+					  * degrees of freedom in the
+					  * DoF handler argument.
+					  */
+	template <class DH>
+	static void
+	compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
+				const DH&,
+				const bool reversed_numbering = false,
+				const bool use_constraints    = false);
     };
     
 				     /**
diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc
index 1d93fd6cc5..b0e1eebd95 100644
--- a/deal.II/deal.II/source/dofs/dof_renumbering.cc
+++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc
@@ -39,6 +39,7 @@
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/cuthill_mckee_ordering.hpp>
 #include <boost/graph/king_ordering.hpp>
+#include <boost/graph/minimum_degree_ordering.hpp>
 #include <boost/graph/properties.hpp>
 #include <boost/graph/bandwidth.hpp>
 
@@ -314,6 +315,126 @@ compute_king_ordering (std::vector<unsigned int>& new_dof_indices,
 }
 
 
+
+template <class DH>
+void
+DoFRenumbering::boost::
+minimum_degree (DH&              dof_handler,
+		const bool       reversed_numbering,
+		const bool       use_constraints)
+{
+  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+					DH::invalid_dof_index);
+  compute_minimum_degree(renumbering, dof_handler, reversed_numbering,
+			 use_constraints);
+
+				   // actually perform renumbering;
+				   // this is dimension specific and
+				   // thus needs an own function
+  dof_handler.renumber_dofs (renumbering);
+}
+
+
+template <class DH>
+void
+DoFRenumbering::boost::
+compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
+			const DH        &dof_handler,
+			const bool       reversed_numbering,
+			const bool       use_constraints)
+{
+  Assert (use_constraints == false, ExcNotImplemented());
+
+				   // the following code is pretty
+				   // much a verbatim copy of the
+				   // sample code for the
+				   // minimum_degree_ordering manual
+				   // page from the BOOST Graph
+				   // Library
+  using namespace boost;
+
+  int delta = 0;
+
+  typedef double Type;
+
+				   // must be BGL directed graph now
+  typedef adjacency_list<vecS, vecS, directedS>  Graph;
+  typedef graph_traits<Graph>::vertex_descriptor Vertex;
+
+  int n = dof_handler.n_dofs();
+
+  Graph G(n);
+
+  std::vector<unsigned int> dofs_on_this_cell;
+
+  typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+				    endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+
+      const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+      dofs_on_this_cell.resize (dofs_per_cell);
+
+      cell->get_dof_indices (dofs_on_this_cell);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	for (unsigned int j=0; j<dofs_per_cell; ++j)
+	  if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
+	    {
+	      add_edge (dofs_on_this_cell[i], dofs_on_this_cell[j], G);
+	      add_edge (dofs_on_this_cell[j], dofs_on_this_cell[i], G);
+	    }
+    }
+  
+  
+  typedef std::vector<int> Vector;
+
+
+  Vector inverse_perm(n, 0);
+
+  Vector perm(n, 0);
+
+
+  Vector supernode_sizes(n, 1);
+				   // init has to be 1
+
+  ::boost::property_map<Graph, vertex_index_t>::type
+    id = get(vertex_index, G);
+
+
+  Vector degree(n, 0);
+
+
+    minimum_degree_ordering
+      (G,
+       make_iterator_property_map(&degree[0], id, degree[0]),
+       &inverse_perm[0],
+       &perm[0],
+       make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]),
+       delta, id);
+
+
+    for (int i=0; i<n; ++i)
+      {
+	Assert (std::find (perm.begin(), perm.end(), i)
+		!= inverse_perm.end(),
+		ExcInternalError());
+	Assert (std::find (inverse_perm.begin(), inverse_perm.end(), i)
+		!= inverse_perm.end(),
+		ExcInternalError());
+	Assert (inverse_perm[perm[i]] == i, ExcInternalError());
+      }
+
+    if (reversed_numbering == true)
+      std::copy (perm.begin(), perm.end(),
+		 new_dof_indices.begin());
+    else
+      std::copy (inverse_perm.begin(), inverse_perm.end(),
+		 new_dof_indices.begin());
+}
+
+
   
 
 
@@ -1714,6 +1835,15 @@ void
 DoFRenumbering::boost::compute_king_ordering (std::vector<unsigned int> &,
 					      const DoFHandler<deal_II_dimension> &, bool, bool);
 
+template
+void
+DoFRenumbering::boost::minimum_degree (DoFHandler<deal_II_dimension> &, bool, bool);
+
+template
+void
+DoFRenumbering::boost::compute_minimum_degree (std::vector<unsigned int> &,
+					       const DoFHandler<deal_II_dimension> &, bool, bool);
+
 
 // non-boost functions:
 
-- 
2.39.5