]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement hierarchical reordering of a SparsityPattern
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 10 Apr 2015 08:27:06 +0000 (10:27 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 11 Apr 2015 06:22:14 +0000 (08:22 +0200)
include/deal.II/lac/sparsity_tools.h
source/lac/sparsity_tools.cc
tests/lac/sparsity_tools_02.cc [new file with mode: 0644]
tests/lac/sparsity_tools_02.output [new file with mode: 0644]

index 1c5fb872485591ea9091f571caff4cd6c5a7c3aa..d82e5ec8e771e91c2b37747b8f8e49fac5a964cc 100644 (file)
@@ -144,6 +144,30 @@ namespace SparsityTools
                          std::vector<SparsityPattern::size_type> &new_indices,
                          const std::vector<SparsityPattern::size_type> &starting_indices = std::vector<SparsityPattern::size_type>()) DEAL_II_DEPRECATED;
 
+  /**
+   * For a given sparsity pattern, compute a re-enumeration of row/column
+   * indices in a hierarchical way, similar to what
+   * DoFRenumbering::hierarchical does for degrees of freedom on
+   * hierarchically refined meshes.
+   *
+   * This algorithm first selects a node with the minimum number of neighbors
+   * and puts that node and its direct neighbors into one chunk. Next, it
+   * selects one of the neighbors of the already selected nodes, adds the node
+   * and its direct neighbors that are not part of one of the previous chunks,
+   * into the next. After this sweep, neighboring nodes are grouped
+   * together. To ensure a similar grouping on a more global level, this
+   * grouping is called recursively on the groups so formed. The recursion
+   * stops when no further grouping is possible. Eventually, the ordering
+   * obtained by this method passes through the indices represented in the
+   * sparsity pattern in a z-like way.
+   *
+   * If the graph has two or more unconnected components, the algorithm will
+   * number each component consecutively, starting with the components with
+   * the lowest number of nodes.
+   */
+  void
+  reorder_hierarchical (const DynamicSparsityPattern                   &sparsity,
+                        std::vector<DynamicSparsityPattern::size_type> &new_indices);
 
 #ifdef DEAL_II_WITH_MPI
   /**
index 1fba415d92fffe1ba465a3599c4632bc65e58a11..905edda0c50f9ec5eaa985471d122266c0d60263 100644 (file)
@@ -21,6 +21,7 @@
 
 #include <algorithm>
 #include <functional>
+#include <set>
 
 #ifdef DEAL_II_WITH_MPI
 #include <deal.II/base/utilities.h>
@@ -341,7 +342,7 @@ namespace SparsityTools
     for (unsigned int row=0; row<sparsity.n_rows(); ++row)
       {
         for (SparsityPattern::iterator it=sparsity.begin(row); it!=sparsity.end(row)
-               && it->is_valid_entry() ; ++it)
+             && it->is_valid_entry() ; ++it)
           dsp.add(row, it->column());
       }
     reorder_Cuthill_McKee(dsp, new_indices, starting_indices);
@@ -349,6 +350,159 @@ namespace SparsityTools
 
 
 
+  namespace internal
+  {
+    void
+    reorder_hierarchical (const DynamicSparsityPattern                   &connectivity,
+                          std::vector<DynamicSparsityPattern::size_type> &renumbering)
+    {
+      AssertDimension (connectivity.n_rows(), connectivity.n_cols());
+      AssertDimension (connectivity.n_rows(), renumbering.size());
+
+      std::vector<types::global_dof_index> touched_nodes(connectivity.n_rows(),
+                                                         numbers::invalid_dof_index);
+      std::set<types::global_dof_index> current_neighbors;
+      std::vector<std::vector<types::global_dof_index> > groups;
+
+      // This outer loop is typically traversed only once, unless the global
+      // graph is not connected
+      while (true)
+        {
+          // Find cell with the minimal number of neighbors (typically a corner
+          // node when based on FEM meshes). If no cell is left, we are done.
+          std::pair<types::global_dof_index,types::global_dof_index> min_neighbors
+          (numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);
+          for (types::global_dof_index i=0; i<touched_nodes.size(); ++i)
+            if (touched_nodes[i] == numbers::invalid_unsigned_int)
+              if (connectivity.row_length(i) < min_neighbors.second)
+                min_neighbors = std::make_pair(i, connectivity.row_length(i));
+          if (min_neighbors.first == numbers::invalid_unsigned_int)
+            break;
+
+          current_neighbors.clear();
+          current_neighbors.insert(min_neighbors.first);
+          while (!current_neighbors.empty())
+            {
+              // Find cell with minimum number of untouched neighbors among the
+              // next set of possible neighbors
+              min_neighbors = std::make_pair (numbers::invalid_unsigned_int,
+                                              numbers::invalid_unsigned_int);
+              for (std::set<types::global_dof_index>::iterator it=current_neighbors.begin();
+                   it != current_neighbors.end(); ++it)
+                {
+                  AssertThrow(touched_nodes[*it] == numbers::invalid_unsigned_int,
+                              ExcInternalError());
+                  types::global_dof_index active_row_length = 0;
+                  for (CompressedSimpleSparsityPattern::row_iterator rowit
+                       = connectivity.row_begin(*it);
+                       rowit != connectivity.row_end(*it); ++rowit)
+                    if (touched_nodes[*rowit] == numbers::invalid_unsigned_int)
+                      ++active_row_length;
+                  if (active_row_length < min_neighbors.second)
+                    min_neighbors = std::make_pair(*it, active_row_length);
+                }
+              // Among the set of cells with the minimal number of neighbors,
+              // choose the one with the largest number of touched neighbors,
+              // i.e., the one with the largest row length
+              const types::global_dof_index best_row_length = min_neighbors.second;
+              for (std::set<types::global_dof_index>::iterator it=current_neighbors.begin();
+                   it != current_neighbors.end(); ++it)
+                {
+                  types::global_dof_index active_row_length = 0;
+                  for (CompressedSimpleSparsityPattern::row_iterator rowit
+                       = connectivity.row_begin(*it);
+                       rowit != connectivity.row_end(*it); ++rowit)
+                    if (touched_nodes[*rowit] == numbers::invalid_unsigned_int)
+                      ++active_row_length;
+                  if (active_row_length == best_row_length)
+                    if (connectivity.row_length(*it) > min_neighbors.second)
+                      min_neighbors = std::make_pair(*it, connectivity.row_length(*it));
+                }
+
+              // Add the pivot cell to the current list
+              groups.push_back(std::vector<types::global_dof_index>());
+              std::vector<types::global_dof_index> &new_entries = groups.back();
+              new_entries.push_back(min_neighbors.first);
+              touched_nodes[min_neighbors.first] = groups.size()-1;
+
+              // Add all direct neighbors of the pivot cell not yet touched to
+              // the current list
+              for (CompressedSimpleSparsityPattern::row_iterator it
+                   = connectivity.row_begin(min_neighbors.first);
+                   it != connectivity.row_end(min_neighbors.first); ++it)
+                {
+                  if (touched_nodes[*it] == numbers::invalid_unsigned_int)
+                    {
+                      new_entries.push_back(*it);
+                      touched_nodes[*it] = groups.size()-1;
+                    }
+                }
+
+              // Add all neighbors of the current list not yet touched to the
+              // set of possible next pivots. Delete the entries of the current
+              // list from the set of possible next pivots.
+              for (types::global_dof_index i=0; i<new_entries.size(); ++i)
+                {
+                  for (CompressedSimpleSparsityPattern::row_iterator it
+                       = connectivity.row_begin(new_entries[i]);
+                       it != connectivity.row_end(new_entries[i]); ++it)
+                    if (touched_nodes[*it] == numbers::invalid_unsigned_int)
+                      current_neighbors.insert(*it);
+                  current_neighbors.erase(new_entries[i]);
+                }
+            }
+        }
+      // If the number of groups is smaller than the number of nodes, we can
+      // continue by recursively calling this method
+      if (groups.size() < connectivity.n_rows())
+        {
+          // Form the connectivity of the groups
+          DynamicSparsityPattern connectivity_next(groups.size(),
+                                                   groups.size());
+          for (types::global_dof_index i=0; i<groups.size(); ++i)
+            for (types::global_dof_index col=0; col<groups[i].size(); ++col)
+              for (CompressedSimpleSparsityPattern::row_iterator it
+                   = connectivity.row_begin(groups[i][col]);
+                   it != connectivity.row_end(groups[i][col]); ++it)
+                {
+                  if (touched_nodes[*it] != i)
+                    connectivity_next.add(i, touched_nodes[*it]);
+                }
+
+          // Recursively call the reordering
+          std::vector<types::global_dof_index> renumbering_next(groups.size());
+          reorder_hierarchical(connectivity_next, renumbering_next);
+
+          // Renumber the indices group by group according to the incoming
+          // ordering for the groups
+          for (types::global_dof_index i=0,count=0; i<groups.size(); ++i)
+            for (types::global_dof_index col=0; col<groups[renumbering_next[i]].size(); ++col, ++count)
+              renumbering[count] = groups[renumbering_next[i]][col];
+        }
+      else
+        {
+          // All groups should have size one and no more recursion is possible,
+          // so use the numbering of the groups
+          for (types::global_dof_index i=0,count=0; i<groups.size(); ++i)
+            for (types::global_dof_index col=0; col<groups[i].size(); ++col, ++count)
+              renumbering[count] = groups[i][col];
+        }
+    }
+  }
+
+  void
+  reorder_hierarchical (const DynamicSparsityPattern                   &connectivity,
+                        std::vector<DynamicSparsityPattern::size_type> &renumbering)
+  {
+    // the internal renumbering keeps the numbering the wrong way around (but
+    // we cannot invert the numbering inside that method because it is used
+    // recursively), so invert it here
+    internal::reorder_hierarchical(connectivity, renumbering);
+    renumbering = Utilities::invert_permutation(renumbering);
+  }
+
+
+
 #ifdef DEAL_II_WITH_MPI
   template <class DSP_t>
   void distribute_sparsity_pattern(DSP_t &dsp,
diff --git a/tests/lac/sparsity_tools_02.cc b/tests/lac/sparsity_tools_02.cc
new file mode 100644 (file)
index 0000000..8987406
--- /dev/null
@@ -0,0 +1,67 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// apply SparsityTools::reorder_hierarchical to a graph that consists
+// of two non-connected parts.
+
+#include "../tests.h"
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include <iomanip>
+#include <fstream>
+
+int main ()
+{
+  std::ofstream logfile("output");
+  logfile.setf(std::ios::fixed);
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  DynamicSparsityPattern dsp (8,8);
+  for (unsigned int i=0; i<8; ++i)
+    dsp.add (i,i);
+
+  // create a graph with components 0,2 and 1,3 that are disconnected
+  dsp.add (0,2);
+  dsp.add (2,0);
+
+  dsp.add (1,3);
+  dsp.add (3,1);
+
+  // couple all indices between 3 and 7 except 4
+  for (unsigned int i=3; i<7; ++i)
+    for (unsigned int j=3; j<7; ++j)
+      if (i != 4 && j != 4)
+        dsp.add(i,j);
+
+  // indices 4,7 couple
+  dsp.add (4,7);
+  dsp.add (7,4);
+
+  // now find permutation
+  std::vector<types::global_dof_index> permutation(8);
+  SparsityTools::reorder_hierarchical (dsp, permutation);
+
+  for (unsigned int i=0; i<permutation.size(); ++i)
+    deallog << permutation[i] << std::endl;
+}
+
+
+
diff --git a/tests/lac/sparsity_tools_02.output b/tests/lac/sparsity_tools_02.output
new file mode 100644 (file)
index 0000000..6b70711
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL::0
+DEAL::4
+DEAL::1
+DEAL::5
+DEAL::2
+DEAL::6
+DEAL::7
+DEAL::3

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.