From: Martin Kronbichler Date: Fri, 10 Apr 2015 08:27:06 +0000 (+0200) Subject: Implement hierarchical reordering of a SparsityPattern X-Git-Tag: v8.3.0-rc1~299^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=14b8b43ec5633a46419aa1922b2eb2b98391e17e;p=dealii.git Implement hierarchical reordering of a SparsityPattern --- diff --git a/include/deal.II/lac/sparsity_tools.h b/include/deal.II/lac/sparsity_tools.h index 1c5fb87248..d82e5ec8e7 100644 --- a/include/deal.II/lac/sparsity_tools.h +++ b/include/deal.II/lac/sparsity_tools.h @@ -144,6 +144,30 @@ namespace SparsityTools std::vector &new_indices, const std::vector &starting_indices = std::vector()) 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 &new_indices); #ifdef DEAL_II_WITH_MPI /** diff --git a/source/lac/sparsity_tools.cc b/source/lac/sparsity_tools.cc index 1fba415d92..905edda0c5 100644 --- a/source/lac/sparsity_tools.cc +++ b/source/lac/sparsity_tools.cc @@ -21,6 +21,7 @@ #include #include +#include #ifdef DEAL_II_WITH_MPI #include @@ -341,7 +342,7 @@ namespace SparsityTools for (unsigned int row=0; rowis_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 &renumbering) + { + AssertDimension (connectivity.n_rows(), connectivity.n_cols()); + AssertDimension (connectivity.n_rows(), renumbering.size()); + + std::vector touched_nodes(connectivity.n_rows(), + numbers::invalid_dof_index); + std::set current_neighbors; + std::vector > 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 min_neighbors + (numbers::invalid_unsigned_int, numbers::invalid_unsigned_int); + for (types::global_dof_index i=0; i::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::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()); + std::vector &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 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 &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 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 index 0000000000..89874064c7 --- /dev/null +++ b/tests/lac/sparsity_tools_02.cc @@ -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 +#include + +#include +#include + +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 permutation(8); + SparsityTools::reorder_hierarchical (dsp, permutation); + + for (unsigned int i=0; i