ViennaCL - The Vienna Computing Library  1.2.0
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-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
20 
27 #include <iostream>
28 #include <fstream>
29 #include <string>
30 #include <algorithm>
31 #include <map>
32 #include <vector>
33 #include <deque>
34 #include <cmath>
35 
36 
37 namespace viennacl
38 {
39 
40  namespace detail
41  {
42 
43  // function to calculate the increment of combination comb.
44  // parameters:
45  // comb: pointer to vector<int> of size m, m <= n
46  // 1 <= comb[i] <= n for 0 <= i < m
47  // comb[i] < comb[i+1] for 0 <= i < m - 1
48  // comb represents an unordered selection of m values out of n
49  // n: int
50  // total number of values out of which comb is taken as selection
51  inline bool comb_inc(std::vector<int> & comb, int n)
52  {
53  int m;
54  int k;
55 
56  m = comb.size();
57  // calculate k as highest possible index such that (*comb)[k-1] can be incremented
58  k = m;
59  while ( (k > 0) && ( ((k == m) && (comb[k-1] == n)) ||
60  ((k < m) && (comb[k-1] == comb[k] - 1) )) )
61  {
62  k--;
63  }
64  if (k == 0) // no further increment of comb possible -> return false
65  {
66  return false;
67  }
68  else
69  {
70  (comb[k-1])++; // increment (*comb)[k-1],
71  for (int i = k; i < m; i++) // and all higher index positions of comb are
72  // calculated just as directly following integer values (lowest possible values)
73  {
74  comb[i] = comb[k-1] + (i - k + 1);
75  }
76  return true;
77  }
78  }
79 
80 
81  // function to generate a node layering as a tree structure rooted at
82  // node s
83  template <typename MatrixType>
84  void generate_layering(MatrixType const & matrix,
85  std::vector< std::vector<int> > & l,
86  int s)
87  {
88  std::size_t n = matrix.size();
89  //std::vector< std::vector<int> > l;
90  std::vector<bool> inr(n, false);
91  std::vector<int> nlist;
92 
93  nlist.push_back(s);
94  inr[s] = true;
95  l.push_back(nlist);
96 
97  for (;;)
98  {
99  nlist.clear();
100  for (std::vector<int>::iterator it = l.back().begin();
101  it != l.back().end();
102  it++)
103  {
104  for (typename MatrixType::value_type::const_iterator it2 = matrix[*it].begin();
105  it2 != matrix[*it].end();
106  it2++)
107  {
108  if (it2->first == *it) continue;
109  if (inr[it2->first]) continue;
110 
111  nlist.push_back(it2->first);
112  inr[it2->first] = true;
113  }
114  }
115 
116  if (nlist.size() == 0)
117  break;
118 
119  l.push_back(nlist);
120  }
121 
122  }
123 
124 
125  // comparison function for comparing two vector<int> values by their
126  // [1]-element
127  inline bool cuthill_mckee_comp_func(std::vector<int> const & a,
128  std::vector<int> const & b)
129  {
130  return (a[1] < b[1]);
131  }
132 
133  }
134 
135  //
136  // Part 1: The original Cuthill-McKee algorithm
137  //
138 
139 
140  struct cuthill_mckee_tag {};
141 
156  template <typename MatrixType>
157  std::vector<int> reorder(MatrixType const & matrix, cuthill_mckee_tag)
158  {
159  std::size_t n = matrix.size();
160  std::vector<int> r;
161  std::vector<bool> inr(n, false); // status array which remembers which nodes have been added to r
162  std::deque<int> q;
163  std::vector< std::vector<int> > nodes;
164  std::vector<int> tmp(2);
165  int p = 0;
166  int c;
167 
168  int deg;
169  int deg_min;
170 
171  r.reserve(n);
172  nodes.reserve(n);
173 
174  do
175  {
176  // under all nodes not yet in r determine one with minimal degree
177  deg_min = -1;
178  for (std::size_t i = 0; i < n; i++)
179  {
180  if (!inr[i])
181  {
182  deg = matrix[i].size() - 1; // node degree
183  if (deg_min < 0 || deg < deg_min)
184  {
185  p = i; // node number
186  deg_min = deg;
187  }
188  }
189  }
190  q.push_back(p); // push that node p with minimal degree on q
191 
192  do
193  {
194  c = q.front();
195  q.pop_front();
196  if (!inr[c])
197  {
198  r.push_back(c);
199  inr[c] = true;
200 
201  // add all neighbouring nodes of c which are not yet in r to nodes
202  nodes.resize(0);
203  for (typename MatrixType::value_type::const_iterator it = matrix[c].begin(); it != matrix[c].end(); it++)
204  {
205  if (it->first == c) continue;
206  if (inr[it->first]) continue;
207 
208  tmp[0] = it->first;
209  tmp[1] = matrix[it->first].size() - 1;
210  nodes.push_back(tmp);
211  }
212 
213  // sort nodes by node degree
214  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
215  for (std::vector< std::vector<int> >::iterator it = nodes.begin(); it != nodes.end(); it++)
216  {
217  q.push_back((*it)[0]);
218  }
219  }
220  } while (q.size() != 0);
221  } while (r.size() < n);
222 
223  return r;
224  }
225 
226 
227  //
228  // Part 2: Advanced Cuthill McKee
229  //
230 
233  {
234  public:
253  advanced_cuthill_mckee_tag(double a = 0.0, std::size_t gmax = 1) : starting_node_param_(a), max_root_nodes_(gmax) {}
254 
255  double starting_node_param() const { return starting_node_param_;}
256  void starting_node_param(double a) { if (a >= 0) starting_node_param_ = a; }
257 
258  std::size_t max_root_nodes() const { return max_root_nodes_; }
259  void max_root_nodes(std::size_t gmax) { max_root_nodes_ = gmax; }
260 
261  private:
262  double starting_node_param_;
263  std::size_t max_root_nodes_;
264  };
265 
266 
267 
276  template <typename MatrixType>
277  std::vector<int> reorder(MatrixType const & matrix,
278  advanced_cuthill_mckee_tag const & tag)
279  {
280  std::size_t n = matrix.size();
281  double a = tag.starting_node_param();
282  std::size_t gmax = tag.max_root_nodes();
283  std::vector<int> r;
284  std::vector<int> r_tmp;
285  std::vector<int> r_best;
286  std::vector<int> r2(n);
287  std::vector<bool> inr(n, false);
288  std::vector<bool> inr_tmp(n);
289  std::vector<bool> inr_best(n);
290  std::deque<int> q;
291  std::vector< std::vector<int> > nodes;
292  std::vector<int> nodes_p;
293  std::vector<int> tmp(2);
294  std::vector< std::vector<int> > l;
295  int deg_min;
296  int deg_max;
297  int deg_a;
298  int deg;
299  int bw;
300  int bw_best;
301  std::vector<int> comb;
302  std::size_t g;
303  int c;
304 
305  r.reserve(n);
306  r_tmp.reserve(n);
307  r_best.reserve(n);
308  nodes.reserve(n);
309  nodes_p.reserve(n);
310  comb.reserve(n);
311 
312  do
313  {
314  // add to nodes_p all nodes not yet in r which are candidates for the root node layer
315  // search unnumbered node and generate layering
316  for (std::size_t i = 0; i < n; i++)
317  {
318  if (!inr[i])
319  {
320  detail::generate_layering(matrix, l, i);
321  break;
322  }
323  }
324  nodes.resize(0);
325  for (std::vector< std::vector<int> >::iterator it = l.begin();
326  it != l.end(); it++)
327  {
328  for (std::vector<int>::iterator it2 = it->begin();
329  it2 != it->end(); it2++)
330  {
331  tmp[0] = *it2;
332  tmp[1] = matrix[*it2].size() - 1;
333  nodes.push_back(tmp);
334  }
335  }
336  // determine minimum and maximum node degree
337  deg_min = -1;
338  deg_max = -1;
339  for (std::vector< std::vector<int> >::iterator it = nodes.begin();
340  it != nodes.end(); it++)
341  {
342  deg = (*it)[1];
343  if (deg_min < 0 || deg < deg_min)
344  {
345  deg_min = deg;
346  }
347  if (deg_max < 0 || deg > deg_max)
348  {
349  deg_max = deg;
350  }
351  }
352  deg_a = deg_min + (int) (a * (deg_max - deg_min));
353  nodes_p.resize(0);
354  for (std::vector< std::vector<int> >::iterator it = nodes.begin();
355  it != nodes.end(); it++)
356  {
357  if ((*it)[1] <= deg_a)
358  {
359  nodes_p.push_back((*it)[0]);
360  }
361  }
362 
363  inr_tmp = inr;
364  g = 1;
365  comb.resize(1);
366  comb[0] = 1;
367  bw_best = -1;
368 
369  for (;;) // for all combinations of g <= gmax root nodes repeat
370  {
371  inr = inr_tmp;
372  r_tmp.resize(0);
373 
374  // add the selected root nodes according to actual combination comb to q
375  for (std::vector<int>::iterator it = comb.begin();
376  it != comb.end(); it++)
377  {
378  q.push_back(nodes_p[(*it)-1]);
379  }
380 
381  do // perform normal CutHill-McKee algorithm for given root nodes with
382  // resulting numbering stored in r_tmp
383  {
384  c = q.front();
385  q.pop_front();
386  if (!inr[c])
387  {
388  r_tmp.push_back(c);
389  inr[c] = true;
390 
391  nodes.resize(0);
392  for (typename MatrixType::value_type::const_iterator it = matrix[c].begin(); it != matrix[c].end(); it++)
393  {
394  if (it->first == c) continue;
395  if (inr[it->first]) continue;
396 
397  tmp[0] = it->first;
398  tmp[1] = matrix[it->first].size() - 1;
399  nodes.push_back(tmp);
400  }
401  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
402  for (std::vector< std::vector<int> >::iterator it =
403  nodes.begin(); it != nodes.end(); it++)
404  {
405  q.push_back((*it)[0]);
406  }
407  }
408  } while (q.size() != 0);
409 
410  // calculate resulting bandwith for root node combination
411  // comb for current numbered component of the node graph
412  for (std::size_t i = 0; i < r_tmp.size(); i++)
413  {
414  r2[r_tmp[i]] = r.size() + i;
415  }
416  bw = 0;
417  for (std::size_t i = 0; i < r_tmp.size(); i++)
418  {
419  for (typename MatrixType::value_type::const_iterator it = matrix[r_tmp[i]].begin();
420  it != matrix[r_tmp[i]].end();
421  it++)
422  {
423  bw = std::max(bw, std::abs(static_cast<int>(r.size() + i) - r2[it->first]));
424  }
425  }
426 
427  // remember ordering r_tmp in r_best for smallest bandwith
428  if (bw_best < 0 || bw < bw_best)
429  {
430  r_best = r_tmp;
431  bw_best = bw;
432  inr_best = inr;
433  }
434 
435  // calculate next combination comb, if not existing
436  // increment g if g stays <= gmax, or else terminate loop
437  if (!detail::comb_inc(comb, nodes_p.size()))
438  {
439  g++;
440  if ( (gmax > 0 && g > gmax) || g > nodes_p.size())
441  {
442  break;
443  }
444  comb.resize(g);
445  for (std::size_t i = 0; i < g; i++)
446  {
447  comb[i] = i + 1;
448  }
449  }
450  }
451 
452  // store best order r_best in result array r
453  for (std::vector<int>::iterator it = r_best.begin();
454  it != r_best.end(); it++)
455  {
456  r.push_back((*it));
457  }
458  inr = inr_best;
459 
460  } while (r.size() < n);
461 
462  return r;
463  }
464 
465 } //namespace viennacl
466 
467 
468 #endif