ViennaCL - The Vienna Computing Library  1.5.2
gibbs_poole_stockmeyer.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_HPP
2 #define VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_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 <fstream>
30 #include <string>
31 #include <algorithm>
32 #include <map>
33 #include <vector>
34 #include <deque>
35 #include <cmath>
36 
37 #include "viennacl/forwards.h"
38 
40 
41 namespace viennacl
42 {
43 
44  namespace detail
45  {
46 
47  // calculates width of a node layering
48  inline int calc_layering_width(std::vector< std::vector<int> > const & l)
49  {
50  int w;
51 
52  w = 0;
53  for (vcl_size_t i = 0; i < l.size(); i++)
54  {
55  w = std::max(w, static_cast<int>(l[i].size()));
56  }
57 
58  return w;
59  }
60 
61  // function to decompose a list of nodes rg into connected components
62  // sorted by decreasing number of nodes per component
63  template <typename MatrixType>
64  std::vector< std::vector<int> > gps_rg_components(MatrixType const & matrix, int n,
65  std::vector<int> const & rg)
66  {
67  std::vector< std::vector<int> > rgc;
68  std::vector< std::vector<int> > rgc_sorted;
69  std::vector< std::vector<int> > sort_ind;
70  std::vector<int> ind(2);
71  std::vector<int> tmp;
72  int c;
73  std::vector<bool> inr(n, true);
74  std::deque<int> q;
75 
76  for (vcl_size_t i = 0; i < rg.size(); i++)
77  {
78  inr[rg[i]] = false;
79  }
80 
81  do
82  {
83  for (int i = 0; i < n; i++)
84  {
85  if (!inr[i])
86  {
87  q.push_front(i);
88  break;
89  }
90  }
91  if (q.size() == 0)
92  {
93  break;
94  }
95 
96  tmp.resize(0);
97  while (q.size() > 0)
98  {
99  c = q.front();
100  q.pop_front();
101 
102  if (!inr[c])
103  {
104  tmp.push_back(c);
105  inr[c] = true;
106 
107  for (typename MatrixType::value_type::const_iterator it = matrix[c].begin(); it != matrix[c].end(); it++)
108  {
109  if (it->first == c) continue;
110  if (inr[it->first]) continue;
111 
112  q.push_back(it->first);
113  }
114  }
115  }
116  rgc.push_back(tmp);
117  } while (true);
118 
119  for (vcl_size_t i = 0; i < rgc.size(); i++)
120  {
121  ind[0] = static_cast<int>(i);
122  ind[1] = static_cast<int>(rgc[i].size());
123  sort_ind.push_back(ind);
124  }
125  std::sort(sort_ind.begin(), sort_ind.end(), detail::cuthill_mckee_comp_func);
126  for (vcl_size_t i = 0; i < rgc.size(); i++)
127  {
128  rgc_sorted.push_back(rgc[sort_ind[rgc.size()-1-i][0]]);
129  }
130 
131  return rgc_sorted;
132  }
133 
134  } // namespace detail
135 
136 
139 
140 
154  template <typename MatrixType>
155  std::vector<int> reorder(MatrixType const & matrix,
157  {
158  vcl_size_t n = matrix.size();
159  std::vector<int> r(n);
160  std::vector< std::vector<int> > rl;
161  vcl_size_t l = 0;
162  int state;
163  bool state_end;
164  std::vector< std::vector<int> > nodes;
165  std::vector<bool> inr(n, false);
166  std::vector<bool> isn(n, false);
167  std::vector<int> tmp(2);
168  int g = 0;
169  int h = 0;
170  std::vector< std::vector<int> > lg;
171  std::vector< std::vector<int> > lh;
172  std::vector< std::vector<int> > ls;
173  std::map< int, std::vector<int> > lap;
174  std::vector<int> rg;
175  std::vector< std::vector<int> > rgc;
176  int m;
177  int m_min;
178  bool new_g = true;
179  int k1, k2, k3, k4;
180  std::vector<int> wvs;
181  std::vector<int> wvsg;
182  std::vector<int> wvsh;
183  int deg_min;
184  int deg;
185  int ind_min;
186 
187  nodes.reserve(n);
188 
189  int current_dof = 0;
190 
191  while (current_dof < static_cast<int>(n)) // for all components of the graph apply GPS algorithm
192  {
193  // determine node g with mimimal degree among all nodes which
194  // are not yet in result array r
195  deg_min = -1;
196  for (vcl_size_t i = 0; i < n; i++)
197  {
198  if (!inr[i])
199  {
200  deg = static_cast<int>(matrix[i].size() - 1); // node degree
201  if (deg_min < 0 || deg < deg_min)
202  {
203  g = static_cast<int>(i); // node number
204  deg_min = deg;
205  }
206  }
207  }
208 
209  // algorithm for determining nodes g, h as endpoints of a pseudo graph diameter
210  while (new_g)
211  {
212  lg.clear();
213  detail::generate_layering(matrix, lg, g);
214 
215  nodes.resize(0);
216  for (vcl_size_t i = 0; i < lg.back().size(); i++)
217  {
218  tmp[0] = lg.back()[i];
219  tmp[1] = static_cast<int>(matrix[lg.back()[i]].size() - 1);
220  nodes.push_back(tmp);
221  }
222  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
223  for (vcl_size_t i = 0; i < nodes.size(); i++)
224  {
225  lg.back()[i] = nodes[i][0];
226  }
227 
228  m_min = -1;
229  new_g = false;
230  for (vcl_size_t i = 0; i < lg.back().size(); i++)
231  {
232  lh.clear();
233  detail::generate_layering(matrix, lh, lg.back()[i]);
234  if (lh.size() > lg.size())
235  {
236  g = lg.back()[i];
237  new_g = true;
238  break;
239  }
241  if (m_min < 0 || m < m_min)
242  {
243  m_min = m;
244  h = lg.back()[i];
245  }
246  }
247  }
248 
249  lh.clear();
250  detail::generate_layering(matrix, lh, h);
251 
252  // calculate ls as layering intersection and rg as remaining
253  // graph
254  lap.clear();
255  for (vcl_size_t i = 0; i < lg.size(); i++)
256  {
257  for (vcl_size_t j = 0; j < lg[i].size(); j++)
258  {
259  lap[lg[i][j]].resize(2);
260  lap[lg[i][j]][0] = static_cast<int>(i);
261  }
262  }
263  for (vcl_size_t i = 0; i < lh.size(); i++)
264  {
265  for (vcl_size_t j = 0; j < lh[i].size(); j++)
266  {
267  lap[lh[i][j]][1] = static_cast<int>(lg.size() - 1 - i);
268  }
269  }
270  rg.clear();
271  ls.clear();
272  ls.resize(lg.size());
273  for (std::map< int, std::vector<int> >::iterator it = lap.begin();
274  it != lap.end(); it++)
275  {
276  if ((it->second)[0] == (it->second)[1])
277  {
278  ls[(it->second)[0]].push_back(it->first);
279  }
280  else
281  {
282  rg.push_back(it->first);
283  }
284  }
285  // partition remaining graph in connected components
286  rgc = detail::gps_rg_components(matrix, static_cast<int>(n), rg);
287 
288  // insert nodes of each component of rgc
291  wvs.resize(ls.size());
292  wvsg.resize(ls.size());
293  wvsh.resize(ls.size());
294  for (vcl_size_t i = 0; i < rgc.size(); i++)
295  {
296  for (vcl_size_t j = 0; j < ls.size(); j++)
297  {
298  wvs[j] = static_cast<int>(ls[j].size());
299  wvsg[j] = static_cast<int>(ls[j].size());
300  wvsh[j] = static_cast<int>(ls[j].size());
301  }
302  for (vcl_size_t j = 0; j < rgc[i].size(); j++)
303  {
304  (wvsg[lap[rgc[i][j]][0]])++;
305  (wvsh[lap[rgc[i][j]][1]])++;
306  }
307  k3 = 0;
308  k4 = 0;
309  for (vcl_size_t j = 0; j < ls.size(); j++)
310  {
311  if (wvsg[j] > wvs[j])
312  {
313  k3 = std::max(k3, wvsg[j]);
314  }
315  if (wvsh[j] > wvs[j])
316  {
317  k4 = std::max(k4, wvsh[j]);
318  }
319  }
320  if (k3 < k4 || (k3 == k4 && k1 <= k2) )
321  {
322  for (vcl_size_t j = 0; j < rgc[i].size(); j++)
323  {
324  ls[lap[rgc[i][j]][0]].push_back(rgc[i][j]);
325  }
326  }
327  else
328  {
329  for (vcl_size_t j = 0; j < rgc[i].size(); j++)
330  {
331  ls[lap[rgc[i][j]][1]].push_back(rgc[i][j]);
332  }
333  }
334  }
335 
336  // renumber nodes in ls
337  rl.clear();
338  rl.resize(ls.size());
339  state = 1;
340  state_end = false;
341  while (!state_end)
342  {
343  switch(state)
344  {
345  case 1:
346  l = 0;
347  state = 4;
348  break;
349 
350  case 2:
351  for (vcl_size_t i = 0; i < rl[l-1].size(); i++)
352  {
353  isn.assign(n, false);
354  for (std::map<int, double>::const_iterator it = matrix[rl[l-1][i]].begin();
355  it != matrix[rl[l-1][i]].end();
356  it++)
357  {
358  if (it->first == rl[l-1][i]) continue;
359  isn[it->first] = true;
360  }
361  nodes.resize(0);
362  for (vcl_size_t j = 0; j < ls[l].size(); j++)
363  {
364  if (inr[ls[l][j]]) continue;
365  if (!isn[ls[l][j]]) continue;
366  tmp[0] = ls[l][j];
367  tmp[1] = static_cast<int>(matrix[ls[l][j]].size() - 1);
368  nodes.push_back(tmp);
369  }
370  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
371  for (vcl_size_t j = 0; j < nodes.size(); j++)
372  {
373  rl[l].push_back(nodes[j][0]);
374  r[nodes[j][0]] = current_dof++;
375  inr[nodes[j][0]] = true;
376  }
377  }
378 
379  case 3:
380  for (vcl_size_t i = 0; i < rl[l].size(); i++)
381  {
382  isn.assign(n, false);
383  for (std::map<int, double>::const_iterator it = matrix[rl[l][i]].begin();
384  it != matrix[rl[l][i]].end();
385  it++)
386  {
387  if (it->first == rl[l][i]) continue;
388  isn[it->first] = true;
389  }
390  nodes.resize(0);
391  for (vcl_size_t j = 0; j < ls[l].size(); j++)
392  {
393  if (inr[ls[l][j]]) continue;
394  if (!isn[ls[l][j]]) continue;
395  tmp[0] = ls[l][j];
396  tmp[1] = static_cast<int>(matrix[ls[l][j]].size() - 1);
397  nodes.push_back(tmp);
398  }
399  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
400  for (vcl_size_t j = 0; j < nodes.size(); j++)
401  {
402  rl[l].push_back(nodes[j][0]);
403  r[nodes[j][0]] = current_dof++;
404  inr[nodes[j][0]] = true;
405  }
406  }
407 
408  case 4:
409  if (rl[l].size() < ls[l].size())
410  {
411  deg_min = -1;
412  for (vcl_size_t j = 0; j < ls[l].size(); j++)
413  {
414  if (inr[ls[l][j]]) continue;
415  deg = static_cast<int>(matrix[ls[l][j]].size() - 1);
416  if (deg_min < 0 || deg < deg_min)
417  {
418  ind_min = ls[l][j];
419  deg_min = deg;
420  }
421  }
422  rl[l].push_back(ind_min);
423  r[ind_min] = current_dof++;
424  inr[ind_min] = true;
425  state = 3;
426  break;
427  }
428 
429  case 5:
430  l++;
431  if (l < ls.size())
432  {
433  state = 2;
434  }
435  else
436  {
437  state_end = true;
438  }
439  break;
440 
441  default:
442  break;
443  }
444  }
445 
446  }
447 
448  return r;
449  }
450 
451 
452 } //namespace viennacl
453 
454 
455 #endif
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
This file provides the forward declarations for the main types used within ViennaCL.
A dense matrix class.
Definition: forwards.h:293
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
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
Tag class for identifying the Gibbs-Poole-Stockmeyer algorithm for reducing the bandwidth of a sparse...
Definition: gibbs_poole_stockmeyer.hpp:138
std::vector< std::vector< int > > gps_rg_components(MatrixType const &matrix, int n, std::vector< int > const &rg)
Definition: gibbs_poole_stockmeyer.hpp:64
Implementation of several flavors of the Cuthill-McKee algorithm. Experimental.
int calc_layering_width(std::vector< std::vector< int > > const &l)
Definition: gibbs_poole_stockmeyer.hpp:48
bool cuthill_mckee_comp_func(std::vector< int > const &a, std::vector< int > const &b)
Definition: cuthill_mckee.hpp:263