ViennaCL - The Vienna Computing Library  1.2.0
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-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 
37 
38 namespace viennacl
39 {
40 
41  namespace detail
42  {
43 
44  // calculates width of a node layering
45  inline int calc_layering_width(std::vector< std::vector<int> > const & l)
46  {
47  int w;
48 
49  w = 0;
50  for (std::size_t i = 0; i < l.size(); i++)
51  {
52  w = std::max(w, static_cast<int>(l[i].size()));
53  }
54 
55  return w;
56  }
57 
58  // function to decompose a list of nodes rg into connected components
59  // sorted by decreasing number of nodes per component
60  template <typename MatrixType>
61  std::vector< std::vector<int> > gps_rg_components(MatrixType const & matrix, int n,
62  std::vector<int> const & rg)
63  {
64  std::vector< std::vector<int> > rgc;
65  std::vector< std::vector<int> > rgc_sorted;
66  std::vector< std::vector<int> > sort_ind;
67  std::vector<int> ind(2);
68  std::vector<int> tmp;
69  int c;
70  std::vector<bool> inr(n, true);
71  std::deque<int> q;
72 
73  for (std::size_t i = 0; i < rg.size(); i++)
74  {
75  inr[rg[i]] = false;
76  }
77 
78  do
79  {
80  for (int i = 0; i < n; i++)
81  {
82  if (!inr[i])
83  {
84  q.push_front(i);
85  break;
86  }
87  }
88  if (q.size() == 0)
89  {
90  break;
91  }
92 
93  tmp.resize(0);
94  while (q.size() > 0)
95  {
96  c = q.front();
97  q.pop_front();
98 
99  if (!inr[c])
100  {
101  tmp.push_back(c);
102  inr[c] = true;
103 
104  for (typename MatrixType::value_type::const_iterator it = matrix[c].begin(); it != matrix[c].end(); it++)
105  {
106  if (it->first == c) continue;
107  if (inr[it->first]) continue;
108 
109  q.push_back(it->first);
110  }
111  }
112  }
113  rgc.push_back(tmp);
114  } while (true);
115 
116  for (std::size_t i = 0; i < rgc.size(); i++)
117  {
118  ind[0] = i;
119  ind[1] = 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 (std::size_t i = 0; i < rgc.size(); i++)
124  {
125  rgc_sorted.push_back(rgc[sort_ind[rgc.size()-1-i][0]]);
126  }
127 
128  return rgc_sorted;
129  }
130 
131  } // namespace detail
132 
133 
135 
136 
150  template <typename MatrixType>
151  std::vector<int> reorder(MatrixType const & matrix,
153  {
154  std::size_t n = matrix.size();
155  std::vector<int> r;
156  std::vector< std::vector<int> > rl;
157  std::size_t l = 0;
158  int state;
159  bool state_end;
160  std::vector< std::vector<int> > nodes;
161  std::vector<bool> inr(n, false);
162  std::vector<bool> isn(n, false);
163  std::vector<int> tmp(2);
164  int g = 0;
165  int h = 0;
166  std::vector< std::vector<int> > lg;
167  std::vector< std::vector<int> > lh;
168  std::vector< std::vector<int> > ls;
169  std::map< int, std::vector<int> > lap;
170  std::vector<int> rg;
171  std::vector< std::vector<int> > rgc;
172  int m;
173  int m_min;
174  bool new_g = true;
175  int k1, k2, k3, k4;
176  std::vector<int> wvs;
177  std::vector<int> wvsg;
178  std::vector<int> wvsh;
179  int deg_min;
180  int deg;
181  int ind_min;
182 
183  r.reserve(n);
184  nodes.reserve(n);
185 
186  while (r.size() < 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 (std::size_t i = 0; i < n; i++)
192  {
193  if (!inr[i])
194  {
195  deg = matrix[i].size() - 1; // node degree
196  if (deg_min < 0 || deg < deg_min)
197  {
198  g = 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 (std::size_t i = 0; i < lg.back().size(); i++)
212  {
213  tmp[0] = lg.back()[i];
214  tmp[1] = matrix[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 (std::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 (std::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 (std::size_t i = 0; i < lg.size(); i++)
251  {
252  for (std::size_t j = 0; j < lg[i].size(); j++)
253  {
254  lap[lg[i][j]].resize(2);
255  lap[lg[i][j]][0] = i;
256  }
257  }
258  for (std::size_t i = 0; i < lh.size(); i++)
259  {
260  for (std::size_t j = 0; j < lh[i].size(); j++)
261  {
262  lap[lh[i][j]][1] = lg.size() - 1 - i;
263  }
264  }
265  rg.clear();
266  ls.clear();
267  ls.resize(lg.size());
268  for (std::map< int, std::vector<int> >::iterator it = lap.begin();
269  it != lap.end(); it++)
270  {
271  if ((it->second)[0] == (it->second)[1])
272  {
273  ls[(it->second)[0]].push_back(it->first);
274  }
275  else
276  {
277  rg.push_back(it->first);
278  }
279  }
280  // partition remaining graph in connected components
281  rgc = detail::gps_rg_components(matrix, n, rg);
282 
283  // insert nodes of each component of rgc
286  wvs.resize(ls.size());
287  wvsg.resize(ls.size());
288  wvsh.resize(ls.size());
289  for (std::size_t i = 0; i < rgc.size(); i++)
290  {
291  for (std::size_t j = 0; j < ls.size(); j++)
292  {
293  wvs[j] = ls[j].size();
294  wvsg[j] = ls[j].size();
295  wvsh[j] = ls[j].size();
296  }
297  for (std::size_t j = 0; j < rgc[i].size(); j++)
298  {
299  (wvsg[lap[rgc[i][j]][0]])++;
300  (wvsh[lap[rgc[i][j]][1]])++;
301  }
302  k3 = 0;
303  k4 = 0;
304  for (std::size_t j = 0; j < ls.size(); j++)
305  {
306  if (wvsg[j] > wvs[j])
307  {
308  k3 = std::max(k3, wvsg[j]);
309  }
310  if (wvsh[j] > wvs[j])
311  {
312  k4 = std::max(k4, wvsh[j]);
313  }
314  }
315  if (k3 < k4 || (k3 == k4 && k1 <= k2) )
316  {
317  for (std::size_t j = 0; j < rgc[i].size(); j++)
318  {
319  ls[lap[rgc[i][j]][0]].push_back(rgc[i][j]);
320  }
321  }
322  else
323  {
324  for (std::size_t j = 0; j < rgc[i].size(); j++)
325  {
326  ls[lap[rgc[i][j]][1]].push_back(rgc[i][j]);
327  }
328  }
329  }
330 
331  // renumber nodes in ls
332  rl.clear();
333  rl.resize(ls.size());
334  state = 1;
335  state_end = false;
336  while (!state_end)
337  {
338  switch(state)
339  {
340  case 1:
341  l = 0;
342  state = 4;
343  break;
344 
345  case 2:
346  for (std::size_t i = 0; i < rl[l-1].size(); i++)
347  {
348  isn.assign(n, false);
349  for (std::map<int, double>::const_iterator it = matrix[rl[l-1][i]].begin();
350  it != matrix[rl[l-1][i]].end();
351  it++)
352  {
353  if (it->first == rl[l-1][i]) continue;
354  isn[it->first] = true;
355  }
356  nodes.resize(0);
357  for (std::size_t j = 0; j < ls[l].size(); j++)
358  {
359  if (inr[ls[l][j]]) continue;
360  if (!isn[ls[l][j]]) continue;
361  tmp[0] = ls[l][j];
362  tmp[1] = matrix[ls[l][j]].size() - 1;
363  nodes.push_back(tmp);
364  }
365  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
366  for (std::size_t j = 0; j < nodes.size(); j++)
367  {
368  rl[l].push_back(nodes[j][0]);
369  r.push_back(nodes[j][0]);
370  inr[nodes[j][0]] = true;
371  }
372  }
373 
374  case 3:
375  for (std::size_t i = 0; i < rl[l].size(); i++)
376  {
377  isn.assign(n, false);
378  for (std::map<int, double>::const_iterator it = matrix[rl[l][i]].begin();
379  it != matrix[rl[l][i]].end();
380  it++)
381  {
382  if (it->first == rl[l][i]) continue;
383  isn[it->first] = true;
384  }
385  nodes.resize(0);
386  for (std::size_t j = 0; j < ls[l].size(); j++)
387  {
388  if (inr[ls[l][j]]) continue;
389  if (!isn[ls[l][j]]) continue;
390  tmp[0] = ls[l][j];
391  tmp[1] = matrix[ls[l][j]].size() - 1;
392  nodes.push_back(tmp);
393  }
394  std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
395  for (std::size_t j = 0; j < nodes.size(); j++)
396  {
397  rl[l].push_back(nodes[j][0]);
398  r.push_back(nodes[j][0]);
399  inr[nodes[j][0]] = true;
400  }
401  }
402 
403  case 4:
404  if (rl[l].size() < ls[l].size())
405  {
406  deg_min = -1;
407  for (std::size_t j = 0; j < ls[l].size(); j++)
408  {
409  if (inr[ls[l][j]]) continue;
410  deg = matrix[ls[l][j]].size() - 1;
411  if (deg_min < 0 || deg < deg_min)
412  {
413  ind_min = ls[l][j];
414  deg_min = deg;
415  }
416  }
417  rl[l].push_back(ind_min);
418  r.push_back(ind_min);
419  inr[ind_min] = true;
420  state = 3;
421  break;
422  }
423 
424  case 5:
425  l++;
426  if (l < ls.size())
427  {
428  state = 2;
429  }
430  else
431  {
432  state_end = true;
433  }
434  break;
435 
436  default:
437  break;
438  }
439  }
440 
441  }
442 
443  return r;
444  }
445 
446 
447 } //namespace viennacl
448 
449 
450 #endif