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