ViennaCL - The Vienna Computing Library  1.5.2
cuthill_mckee.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_MISC_CUTHILL_MCKEE_HPP
2 #define VIENNACL_MISC_CUTHILL_MCKEE_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
21 
28 #include <iostream>
29 #include <iterator>
30 #include <fstream>
31 #include <string>
32 #include <algorithm>
33 #include <map>
34 #include <vector>
35 #include <deque>
36 #include <cmath>
37 
38 #include "viennacl/forwards.h"
39 
40 namespace viennacl
41 {
42 
43  namespace detail
44  {
45 
46  // Calculate the bandwidth of a reordered matrix
47  template <typename IndexT, typename ValueT>
48  IndexT calc_reordered_bw(std::vector< std::map<IndexT, ValueT> > const & matrix,
49  std::vector<bool> & dof_assigned_to_node,
50  std::vector<IndexT> const & permutation)
51  {
52  IndexT bw = 0;
53 
54  for (vcl_size_t i = 0; i < permutation.size(); i++)
55  {
56  if (!dof_assigned_to_node[i])
57  continue;
58 
59  IndexT min_index = static_cast<IndexT>(matrix.size());
60  IndexT max_index = 0;
61  for (typename std::map<IndexT, ValueT>::const_iterator it = matrix[i].begin(); it != matrix[i].end(); it++)
62  {
63  if (!dof_assigned_to_node[it->first])
64  continue;
65 
66  if (permutation[it->first] > max_index)
67  max_index = permutation[it->first];
68  if (permutation[it->first] < min_index)
69  min_index = permutation[it->first];
70  }
71  if (max_index > min_index)
72  bw = std::max(bw, max_index - min_index);
73  }
74 
75  return bw;
76  }
77 
78 
79 
80  // function to calculate the increment of combination comb.
81  // parameters:
82  // comb: pointer to vector<int> of size m, m <= n
83  // 1 <= comb[i] <= n for 0 <= i < m
84  // comb[i] < comb[i+1] for 0 <= i < m - 1
85  // comb represents an ordered selection of m values out of n
86  // n: int
87  // total number of values out of which comb is taken as selection
88  template <typename IndexT>
89  bool comb_inc(std::vector<IndexT> & comb, vcl_size_t n)
90  {
91  IndexT m;
92  IndexT k;
93 
94  m = static_cast<IndexT>(comb.size());
95  // calculate k as highest possible index such that (*comb)[k-1] can be incremented
96  k = m;
97  while ( (k > 0) && ( ((k == m) && (comb[k-1] == static_cast<IndexT>(n)-1)) ||
98  ((k < m) && (comb[k-1] == comb[k] - 1) )) )
99  {
100  k--;
101  }
102 
103  if (k == 0) // no further increment of comb possible -> return false
104  return false;
105 
106  comb[k-1] += 1;
107 
108  // and all higher index positions of comb are calculated just as directly following integer values
109  // Example (1, 4, 7) -> (1, 5, 6) -> (1, 5, 7) -> (1, 6, 7) -> done for n=7
110  for (IndexT i = k; i < m; i++)
111  comb[i] = comb[k-1] + (i - k);
112  return true;
113  }
114 
115 
120  // node s
121  template <typename MatrixT, typename IndexT>
122  void generate_layering(MatrixT const & matrix,
123  std::vector< std::vector<IndexT> > & layer_list)
124  {
125  std::vector<bool> node_visited_already(matrix.size(), false);
126 
127  //
128  // Step 1: Set root nodes to visited
129  //
130  for (vcl_size_t i=0; i<layer_list.size(); ++i)
131  {
132  for (typename std::vector<IndexT>::iterator it = layer_list[i].begin();
133  it != layer_list[i].end();
134  it++)
135  node_visited_already[*it] = true;
136  }
137 
138  //
139  // Step 2: Fill next layers
140  //
141  while (layer_list.back().size() > 0)
142  {
143  vcl_size_t layer_index = layer_list.size(); //parent nodes are at layer 0
144  layer_list.push_back(std::vector<IndexT>());
145 
146  for (typename std::vector<IndexT>::iterator it = layer_list[layer_index].begin();
147  it != layer_list[layer_index].end();
148  it++)
149  {
150  for (typename MatrixT::value_type::const_iterator it2 = matrix[*it].begin();
151  it2 != matrix[*it].end();
152  it2++)
153  {
154  if (it2->first == *it) continue;
155  if (node_visited_already[it2->first]) continue;
156 
157  layer_list.back().push_back(it2->first);
158  node_visited_already[it2->first] = true;
159  }
160  }
161  }
162 
163  // remove last (empty) nodelist:
164  layer_list.resize(layer_list.size()-1);
165  }
166 
167 
168  // function to generate a node layering as a tree structure rooted at node s
169  template <typename MatrixType>
170  void generate_layering(MatrixType const & matrix,
171  std::vector< std::vector<int> > & l,
172  int s)
173  {
174  vcl_size_t n = matrix.size();
175  //std::vector< std::vector<int> > l;
176  std::vector<bool> inr(n, false);
177  std::vector<int> nlist;
178 
179  nlist.push_back(s);
180  inr[s] = true;
181  l.push_back(nlist);
182 
183  for (;;)
184  {
185  nlist.clear();
186  for (std::vector<int>::iterator it = l.back().begin();
187  it != l.back().end();
188  it++)
189  {
190  for (typename MatrixType::value_type::const_iterator it2 = matrix[*it].begin();
191  it2 != matrix[*it].end();
192  it2++)
193  {
194  if (it2->first == *it) continue;
195  if (inr[it2->first]) continue;
196 
197  nlist.push_back(it2->first);
198  inr[it2->first] = true;
199  }
200  }
201 
202  if (nlist.size() == 0)
203  break;
204 
205  l.push_back(nlist);
206  }
207 
208  }
209 
214  template <typename MatrixT, typename IndexT>
216  std::vector<IndexT> & node_list)
217  {
218  std::vector<bool> node_visited_already(matrix.size(), false);
219  std::deque<IndexT> node_queue;
220 
221  //
222  // Step 1: Push root nodes to queue:
223  //
224  for (typename std::vector<IndexT>::iterator it = node_list.begin();
225  it != node_list.end();
226  it++)
227  {
228  node_queue.push_back(*it);
229  }
230  node_list.resize(0);
231 
232  //
233  // Step 2: Fill with remaining nodes of strongly connected compontent
234  //
235  while (!node_queue.empty())
236  {
237  IndexT node_id = node_queue.front();
238  node_queue.pop_front();
239 
240  if (!node_visited_already[node_id])
241  {
242  node_list.push_back(node_id);
243  node_visited_already[node_id] = true;
244 
245  for (typename MatrixT::value_type::const_iterator it = matrix[node_id].begin();
246  it != matrix[node_id].end();
247  it++)
248  {
249  IndexT neighbor_node_id = it->first;
250  if (neighbor_node_id == node_id) continue;
251  if (node_visited_already[neighbor_node_id]) continue;
252 
253  node_queue.push_back(neighbor_node_id);
254  }
255  }
256  }
257 
258  }
259 
260 
261  // comparison function for comparing two vector<int> values by their
262  // [1]-element
263  inline bool cuthill_mckee_comp_func(std::vector<int> const & a,
264  std::vector<int> const & b)
265  {
266  return (a[1] < b[1]);
267  }
268 
269  template <typename IndexT>
270  bool cuthill_mckee_comp_func_pair(std::pair<IndexT, IndexT> const & a,
271  std::pair<IndexT, IndexT> const & b)
272  {
273  return (a.second < b.second);
274  }
275 
286  template <typename IndexT, typename ValueT>
287  vcl_size_t cuthill_mckee_on_strongly_connected_component(std::vector< std::map<IndexT, ValueT> > const & matrix,
288  std::deque<IndexT> & node_assignment_queue,
289  std::vector<bool> & dof_assigned_to_node,
290  std::vector<IndexT> & permutation,
291  vcl_size_t current_dof)
292  {
293  typedef std::pair<IndexT, IndexT> NodeIdDegreePair; //first member is the node ID, second member is the node degree
294 
295  std::vector< NodeIdDegreePair > local_neighbor_nodes(matrix.size());
296 
297  while (!node_assignment_queue.empty())
298  {
299  // Grab first node from queue
300  vcl_size_t node_id = node_assignment_queue.front();
301  node_assignment_queue.pop_front();
302 
303  // Assign dof if a new dof hasn't been assigned yet
304  if (!dof_assigned_to_node[node_id])
305  {
306  permutation[node_id] = static_cast<IndexT>(current_dof); //TODO: Invert this!
307  ++current_dof;
308  dof_assigned_to_node[node_id] = true;
309 
310  //
311  // Get all neighbors of that node:
312  //
313  vcl_size_t num_neighbors = 0;
314  for (typename std::map<IndexT, ValueT>::const_iterator neighbor_it = matrix[node_id].begin();
315  neighbor_it != matrix[node_id].end();
316  ++neighbor_it)
317  {
318  if (!dof_assigned_to_node[neighbor_it->first])
319  {
320  local_neighbor_nodes[num_neighbors] = NodeIdDegreePair(neighbor_it->first, static_cast<IndexT>(matrix[neighbor_it->first].size()));
321  ++num_neighbors;
322  }
323  }
324 
325  // Sort neighbors by increasing node degree
326  std::sort(local_neighbor_nodes.begin(), local_neighbor_nodes.begin() + num_neighbors, detail::cuthill_mckee_comp_func_pair<IndexT>);
327 
328  // Push neighbors to queue
329  for (vcl_size_t i=0; i<num_neighbors; ++i)
330  node_assignment_queue.push_back(local_neighbor_nodes[i].first);
331 
332  } // if node doesn't have a new dof yet
333 
334  } // while nodes in queue
335 
336  return current_dof;
337 
338  }
339 
340  } //namespace detail
341 
342  //
343  // Part 1: The original Cuthill-McKee algorithm
344  //
345 
347  struct cuthill_mckee_tag {};
348 
363  template <typename IndexT, typename ValueT>
364  std::vector<IndexT> reorder(std::vector< std::map<IndexT, ValueT> > const & matrix, cuthill_mckee_tag)
365  {
366  std::vector<IndexT> permutation(matrix.size());
367  std::vector<bool> dof_assigned_to_node(matrix.size(), false); //flag vector indicating whether node i has received a new dof
368  std::deque<IndexT> node_assignment_queue;
369 
370  vcl_size_t current_dof = 0; //the dof to be assigned
371 
372  while (current_dof < matrix.size()) //outer loop for each strongly connected component (there may be more than one)
373  {
374  //
375  // preprocessing: Determine node degrees for nodes which have not been assigned
376  //
377  vcl_size_t current_min_degree = matrix.size();
378  vcl_size_t node_with_minimum_degree = 0;
379  bool found_unassigned_node = false;
380  for (vcl_size_t i=0; i<matrix.size(); ++i)
381  {
382  if (!dof_assigned_to_node[i])
383  {
384  if (matrix[i].size() == 1) //This is an isolated node, so assign DOF right away
385  {
386  permutation[i] = static_cast<IndexT>(current_dof);
387  dof_assigned_to_node[i] = true;
388  ++current_dof;
389  continue;
390  }
391 
392  if (!found_unassigned_node) //initialize minimum degree on first node without new dof
393  {
394  current_min_degree = matrix[i].size();
395  node_with_minimum_degree = i;
396  found_unassigned_node = true;
397  }
398 
399  if (matrix[i].size() < current_min_degree) //found a node with smaller degree
400  {
401  current_min_degree = matrix[i].size();
402  node_with_minimum_degree = i;
403  }
404  }
405  }
406 
407  //
408  // Stage 2: Distribute dofs on this closely connected (sub-)graph in a breath-first manner using one root node
409  //
410  if (found_unassigned_node) // there's work to be done
411  {
412  node_assignment_queue.push_back(static_cast<IndexT>(node_with_minimum_degree));
413  current_dof = detail::cuthill_mckee_on_strongly_connected_component(matrix, node_assignment_queue, dof_assigned_to_node, permutation, current_dof);
414  }
415  }
416 
417  return permutation;
418  }
419 
420 
421  //
422  // Part 2: Advanced Cuthill McKee
423  //
424 
427  {
428  public:
447  advanced_cuthill_mckee_tag(double a = 0.0, vcl_size_t gmax = 1) : starting_node_param_(a), max_root_nodes_(gmax) {}
448 
449  double starting_node_param() const { return starting_node_param_;}
450  void starting_node_param(double a) { if (a >= 0) starting_node_param_ = a; }
451 
452  vcl_size_t max_root_nodes() const { return max_root_nodes_; }
453  void max_root_nodes(vcl_size_t gmax) { max_root_nodes_ = gmax; }
454 
455  private:
456  double starting_node_param_;
457  vcl_size_t max_root_nodes_;
458  };
459 
460 
461 
470  template <typename IndexT, typename ValueT>
471  std::vector<IndexT> reorder(std::vector< std::map<IndexT, ValueT> > const & matrix,
472  advanced_cuthill_mckee_tag const & tag)
473  {
474  vcl_size_t n = matrix.size();
475  double a = tag.starting_node_param();
476  vcl_size_t gmax = tag.max_root_nodes();
477  std::vector<IndexT> permutation(n);
478  std::vector<bool> dof_assigned_to_node(n, false);
479  std::vector<IndexT> nodes_in_strongly_connected_component;
480  std::vector<IndexT> parent_nodes;
481  vcl_size_t deg_min;
482  vcl_size_t deg_max;
483  vcl_size_t deg_a;
484  vcl_size_t deg;
485  std::vector<IndexT> comb;
486 
487  nodes_in_strongly_connected_component.reserve(n);
488  parent_nodes.reserve(n);
489  comb.reserve(n);
490 
491  vcl_size_t current_dof = 0;
492 
493  while (current_dof < matrix.size()) // for all strongly connected components
494  {
495  // get all nodes of the strongly connected component:
496  nodes_in_strongly_connected_component.resize(0);
497  for (vcl_size_t i = 0; i < n; i++)
498  {
499  if (!dof_assigned_to_node[i])
500  {
501  nodes_in_strongly_connected_component.push_back(static_cast<IndexT>(i));
502  detail::nodes_of_strongly_connected_component(matrix, nodes_in_strongly_connected_component);
503  break;
504  }
505  }
506 
507  // determine minimum and maximum node degree
508  deg_min = 0;
509  deg_max = 0;
510  for (typename std::vector<IndexT>::iterator it = nodes_in_strongly_connected_component.begin();
511  it != nodes_in_strongly_connected_component.end();
512  it++)
513  {
514  deg = matrix[*it].size();
515  if (deg_min == 0 || deg < deg_min)
516  deg_min = deg;
517  if (deg_max == 0 || deg > deg_max)
518  deg_max = deg;
519  }
520  deg_a = deg_min + static_cast<vcl_size_t>(a * (deg_max - deg_min));
521 
522  // fill array of parent nodes:
523  parent_nodes.resize(0);
524  for (typename std::vector<IndexT>::iterator it = nodes_in_strongly_connected_component.begin();
525  it != nodes_in_strongly_connected_component.end();
526  it++)
527  {
528  if (matrix[*it].size() <= deg_a)
529  parent_nodes.push_back(*it);
530  }
531 
532  //
533  // backup current state in order to restore for every new combination of parent nodes below
534  //
535  std::vector<bool> dof_assigned_to_node_backup = dof_assigned_to_node;
536  std::vector<bool> dof_assigned_to_node_best;
537 
538  std::vector<IndexT> permutation_backup = permutation;
539  std::vector<IndexT> permutation_best = permutation;
540 
541  vcl_size_t current_dof_backup = current_dof;
542 
543  vcl_size_t g = 1;
544  comb.resize(1);
545  comb[0] = 0;
546 
547  IndexT bw_best = 0;
548 
549  //
550  // Loop over all combinations of g <= gmax root nodes
551  //
552 
553  for (;;)
554  {
555  dof_assigned_to_node = dof_assigned_to_node_backup;
556  permutation = permutation_backup;
557  current_dof = current_dof_backup;
558 
559  std::deque<IndexT> node_queue;
560 
561  // add the selected root nodes according to actual combination comb to q
562  for (typename std::vector<IndexT>::iterator it = comb.begin(); it != comb.end(); it++)
563  node_queue.push_back(parent_nodes[*it]);
564 
565  current_dof = detail::cuthill_mckee_on_strongly_connected_component(matrix, node_queue, dof_assigned_to_node, permutation, current_dof);
566 
567  // calculate resulting bandwith for root node combination
568  // comb for current numbered component of the node graph
569  IndexT bw = detail::calc_reordered_bw(matrix, dof_assigned_to_node, permutation);
570 
571  // remember best ordering:
572  if (bw_best == 0 || bw < bw_best)
573  {
574  permutation_best = permutation;
575  bw_best = bw;
576  dof_assigned_to_node_best = dof_assigned_to_node;
577  }
578 
579  // calculate next combination comb, if not existing
580  // increment g if g stays <= gmax, or else terminate loop
581  if (!detail::comb_inc(comb, parent_nodes.size()))
582  {
583  ++g;
584  if ( (gmax > 0 && g > gmax) || g > parent_nodes.size())
585  break;
586 
587  comb.resize(g);
588  for (vcl_size_t i = 0; i < g; i++)
589  comb[i] = static_cast<IndexT>(i);
590  }
591  }
592 
593  //
594  // restore best permutation
595  //
596  permutation = permutation_best;
597  dof_assigned_to_node = dof_assigned_to_node_best;
598 
599  }
600 
601  return permutation;
602  }
603 
604 
605 } //namespace viennacl
606 
607 
608 #endif
advanced_cuthill_mckee_tag(double a=0.0, vcl_size_t gmax=1)
CTOR which may take the additional parameters for the advanced algorithm.
Definition: cuthill_mckee.hpp:447
std::vector< IndexT > reorder(std::vector< std::map< IndexT, ValueT > > const &matrix, cuthill_mckee_tag)
Function for the calculation of a node number permutation to reduce the bandwidth of an incidence mat...
Definition: cuthill_mckee.hpp:364
std::size_t vcl_size_t
Definition: forwards.h:58
void nodes_of_strongly_connected_component(MatrixT const &matrix, std::vector< IndexT > &node_list)
Fills the provided nodelist with all nodes of the same strongly connected component as the nodes in t...
Definition: cuthill_mckee.hpp:215
A dense matrix class.
Definition: forwards.h:293
This file provides the forward declarations for the main types used within ViennaCL.
void starting_node_param(double a)
Definition: cuthill_mckee.hpp:450
void generate_layering(MatrixT const &matrix, std::vector< std::vector< IndexT > > &layer_list)
Function to generate a node layering as a tree structure.
Definition: cuthill_mckee.hpp:122
void max_root_nodes(vcl_size_t gmax)
Definition: cuthill_mckee.hpp:453
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
vcl_size_t max_root_nodes() const
Definition: cuthill_mckee.hpp:452
double starting_node_param() const
Definition: cuthill_mckee.hpp:449
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
vcl_size_t cuthill_mckee_on_strongly_connected_component(std::vector< std::map< IndexT, ValueT > > const &matrix, std::deque< IndexT > &node_assignment_queue, std::vector< bool > &dof_assigned_to_node, std::vector< IndexT > &permutation, vcl_size_t current_dof)
Runs the Cuthill-McKee algorithm on a strongly connected component of a graph.
Definition: cuthill_mckee.hpp:287
IndexT calc_reordered_bw(std::vector< std::map< IndexT, ValueT > > const &matrix, std::vector< bool > &dof_assigned_to_node, std::vector< IndexT > const &permutation)
Definition: cuthill_mckee.hpp:48
bool comb_inc(std::vector< IndexT > &comb, vcl_size_t n)
Definition: cuthill_mckee.hpp:89
bool cuthill_mckee_comp_func_pair(std::pair< IndexT, IndexT > const &a, std::pair< IndexT, IndexT > const &b)
Definition: cuthill_mckee.hpp:270
Tag for the advanced Cuthill-McKee algorithm (i.e. running the 'standard' Cuthill-McKee algorithm for...
Definition: cuthill_mckee.hpp:426
A tag class for selecting the Cuthill-McKee algorithm for reducing the bandwidth of a sparse matrix...
Definition: cuthill_mckee.hpp:347
bool cuthill_mckee_comp_func(std::vector< int > const &a, std::vector< int > const &b)
Definition: cuthill_mckee.hpp:263