1 #ifndef VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_HPP
2 #define VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_HPP
55 w = std::max(w, static_cast<int>(l[i].
size()));
63 template <
typename MatrixType>
65 std::vector<int>
const & rg)
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);
73 std::vector<bool> inr(n,
true);
83 for (
int i = 0; i < n; i++)
107 for (
typename MatrixType::value_type::const_iterator it = matrix[c].begin(); it != matrix[c].end(); it++)
109 if (it->first == c)
continue;
110 if (inr[it->first])
continue;
112 q.push_back(it->first);
121 ind[0] =
static_cast<int>(i);
122 ind[1] =
static_cast<int>(rgc[i].size());
123 sort_ind.push_back(ind);
128 rgc_sorted.push_back(rgc[sort_ind[rgc.size()-1-i][0]]);
154 template <
typename MatrixType>
159 std::vector<int> r(n);
160 std::vector< std::vector<int> > rl;
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);
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;
175 std::vector< std::vector<int> > rgc;
180 std::vector<int> wvs;
181 std::vector<int> wvsg;
182 std::vector<int> wvsh;
191 while (current_dof < static_cast<int>(n))
200 deg =
static_cast<int>(matrix[i].size() - 1);
201 if (deg_min < 0 || deg < deg_min)
203 g =
static_cast<int>(i);
216 for (
vcl_size_t i = 0; i < lg.back().size(); i++)
218 tmp[0] = lg.back()[i];
219 tmp[1] =
static_cast<int>(matrix[lg.back()[i]].size() - 1);
220 nodes.push_back(tmp);
225 lg.back()[i] = nodes[i][0];
230 for (
vcl_size_t i = 0; i < lg.back().size(); i++)
234 if (lh.size() > lg.size())
241 if (m_min < 0 || m < m_min)
259 lap[lg[i][j]].resize(2);
260 lap[lg[i][j]][0] =
static_cast<int>(i);
267 lap[lh[i][j]][1] =
static_cast<int>(lg.size() - 1 - i);
272 ls.resize(lg.size());
273 for (std::map<
int, std::vector<int> >::iterator it = lap.begin();
274 it != lap.end(); it++)
276 if ((it->second)[0] == (it->second)[1])
278 ls[(it->second)[0]].push_back(it->first);
282 rg.push_back(it->first);
291 wvs.resize(ls.size());
292 wvsg.resize(ls.size());
293 wvsh.resize(ls.size());
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());
302 for (
vcl_size_t j = 0; j < rgc[i].size(); j++)
304 (wvsg[lap[rgc[i][j]][0]])++;
305 (wvsh[lap[rgc[i][j]][1]])++;
311 if (wvsg[j] > wvs[j])
313 k3 = std::max(k3, wvsg[j]);
315 if (wvsh[j] > wvs[j])
317 k4 = std::max(k4, wvsh[j]);
320 if (k3 < k4 || (k3 == k4 && k1 <= k2) )
322 for (
vcl_size_t j = 0; j < rgc[i].size(); j++)
324 ls[lap[rgc[i][j]][0]].push_back(rgc[i][j]);
329 for (
vcl_size_t j = 0; j < rgc[i].size(); j++)
331 ls[lap[rgc[i][j]][1]].push_back(rgc[i][j]);
338 rl.resize(ls.size());
351 for (
vcl_size_t i = 0; i < rl[l-1].size(); i++)
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();
358 if (it->first == rl[l-1][i])
continue;
359 isn[it->first] =
true;
364 if (inr[ls[l][j]])
continue;
365 if (!isn[ls[l][j]])
continue;
367 tmp[1] =
static_cast<int>(matrix[ls[l][j]].size() - 1);
368 nodes.push_back(tmp);
373 rl[l].push_back(nodes[j][0]);
374 r[nodes[j][0]] = current_dof++;
375 inr[nodes[j][0]] =
true;
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();
387 if (it->first == rl[l][i])
continue;
388 isn[it->first] =
true;
393 if (inr[ls[l][j]])
continue;
394 if (!isn[ls[l][j]])
continue;
396 tmp[1] =
static_cast<int>(matrix[ls[l][j]].size() - 1);
397 nodes.push_back(tmp);
402 rl[l].push_back(nodes[j][0]);
403 r[nodes[j][0]] = current_dof++;
404 inr[nodes[j][0]] =
true;
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)
422 rl[l].push_back(ind_min);
423 r[ind_min] = current_dof++;
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