All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gibbs_poole_stockmeyer.hpp
Go to the documentation of this file.
#ifndef VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_HPP
#define VIENNACL_MISC_GIBBS_POOLE_STOCKMEYER_HPP
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <map>
#include <vector>
#include <deque>
#include <cmath>
#include "viennacl/forwards.h"
#include "viennacl/misc/cuthill_mckee.hpp"
namespace viennacl
{
namespace detail
{
// calculates width of a node layering
inline int calc_layering_width(std::vector< std::vector<int> > const & l)
{
int w;
w = 0;
for (vcl_size_t i = 0; i < l.size(); i++)
{
w = std::max(w, static_cast<int>(l[i].size()));
}
return w;
}
// function to decompose a list of nodes rg into connected components
// sorted by decreasing number of nodes per component
template<typename MatrixType>
std::vector< std::vector<int> > gps_rg_components(MatrixType const & matrix, int n,
std::vector<int> const & rg)
{
std::vector< std::vector<int> > rgc;
std::vector< std::vector<int> > rgc_sorted;
std::vector< std::vector<int> > sort_ind;
std::vector<int> ind(2);
std::vector<int> tmp;
int c;
std::vector<bool> inr(static_cast<vcl_size_t>(n), true);
std::deque<int> q;
for (vcl_size_t i = 0; i < rg.size(); i++)
{
inr[static_cast<vcl_size_t>(rg[i])] = false;
}
do
{
for (int i = 0; i < n; i++)
{
if (!inr[static_cast<vcl_size_t>(i)])
{
q.push_front(i);
break;
}
}
if (q.size() == 0)
break;
tmp.resize(0);
while (q.size() > 0)
{
c = q.front();
q.pop_front();
if (!inr[static_cast<vcl_size_t>(c)])
{
tmp.push_back(c);
inr[static_cast<vcl_size_t>(c)] = true;
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++)
{
if (it->first == c) continue;
if (inr[static_cast<vcl_size_t>(it->first)]) continue;
q.push_back(it->first);
}
}
}
rgc.push_back(tmp);
} while (true);
for (vcl_size_t i = 0; i < rgc.size(); i++)
{
ind[0] = static_cast<int>(i);
ind[1] = static_cast<int>(rgc[i].size());
sort_ind.push_back(ind);
}
std::sort(sort_ind.begin(), sort_ind.end(), detail::cuthill_mckee_comp_func);
for (vcl_size_t i = 0; i < rgc.size(); i++)
rgc_sorted.push_back(rgc[static_cast<vcl_size_t>(sort_ind[rgc.size()-1-i][0])]);
return rgc_sorted;
}
} // namespace detail
struct gibbs_poole_stockmeyer_tag {};
template<typename MatrixType>
std::vector<int> reorder(MatrixType const & matrix,
gibbs_poole_stockmeyer_tag)
{
vcl_size_t n = matrix.size();
std::vector<int> r(n);
std::vector< std::vector<int> > rl;
vcl_size_t l = 0;
int state;
bool state_end;
std::vector< std::vector<int> > nodes;
std::vector<bool> inr(n, false);
std::vector<bool> isn(n, false);
std::vector<int> tmp(2);
int g = 0;
int h = 0;
std::vector< std::vector<int> > lg;
std::vector< std::vector<int> > lh;
std::vector< std::vector<int> > ls;
std::map< int, std::vector<int> > lap;
std::vector<int> rg;
std::vector< std::vector<int> > rgc;
int m;
int m_min;
bool new_g = true;
int k1, k2, k3, k4;
std::vector<int> wvs;
std::vector<int> wvsg;
std::vector<int> wvsh;
int deg_min;
int deg;
int ind_min = 0;
nodes.reserve(n);
int current_dof = 0;
while (current_dof < static_cast<int>(n)) // for all components of the graph apply GPS algorithm
{
// determine node g with mimimal degree among all nodes which
// are not yet in result array r
deg_min = -1;
for (vcl_size_t i = 0; i < n; i++)
{
if (!inr[i])
{
deg = static_cast<int>(matrix[i].size() - 1); // node degree
if (deg_min < 0 || deg < deg_min)
{
g = static_cast<int>(i); // node number
deg_min = deg;
}
}
}
// algorithm for determining nodes g, h as endpoints of a pseudo graph diameter
while (new_g)
{
lg.clear();
detail::generate_layering(matrix, lg, g);
nodes.resize(0);
for (vcl_size_t i = 0; i < lg.back().size(); i++)
{
tmp[0] = lg.back()[i];
tmp[1] = static_cast<int>(matrix[static_cast<vcl_size_t>(lg.back()[i])].size() - 1);
nodes.push_back(tmp);
}
std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
for (vcl_size_t i = 0; i < nodes.size(); i++)
{
lg.back()[i] = nodes[i][0];
}
m_min = -1;
new_g = false;
for (vcl_size_t i = 0; i < lg.back().size(); i++)
{
lh.clear();
detail::generate_layering(matrix, lh, lg.back()[i]);
if (lh.size() > lg.size())
{
g = lg.back()[i];
new_g = true;
break;
}
m = detail::calc_layering_width(lh);
if (m_min < 0 || m < m_min)
{
m_min = m;
h = lg.back()[i];
}
}
}
lh.clear();
detail::generate_layering(matrix, lh, h);
// calculate ls as layering intersection and rg as remaining
// graph
lap.clear();
for (vcl_size_t i = 0; i < lg.size(); i++)
{
for (vcl_size_t j = 0; j < lg[i].size(); j++)
{
lap[lg[i][j]].resize(2);
lap[lg[i][j]][0] = static_cast<int>(i);
}
}
for (vcl_size_t i = 0; i < lh.size(); i++)
for (vcl_size_t j = 0; j < lh[i].size(); j++)
lap[lh[i][j]][1] = static_cast<int>(lg.size() - 1 - i);
rg.clear();
ls.clear();
ls.resize(lg.size());
for (std::map< int, std::vector<int> >::iterator it = lap.begin();
it != lap.end(); it++)
{
if ((it->second)[0] == (it->second)[1])
ls[static_cast<vcl_size_t>((it->second)[0])].push_back(it->first);
else
rg.push_back(it->first);
}
// partition remaining graph in connected components
rgc = detail::gps_rg_components(matrix, static_cast<int>(n), rg);
// insert nodes of each component of rgc
k1 = detail::calc_layering_width(lg);
k2 = detail::calc_layering_width(lh);
wvs.resize(ls.size());
wvsg.resize(ls.size());
wvsh.resize(ls.size());
for (vcl_size_t i = 0; i < rgc.size(); i++)
{
for (vcl_size_t j = 0; j < ls.size(); j++)
{
wvs[j] = static_cast<int>(ls[j].size());
wvsg[j] = static_cast<int>(ls[j].size());
wvsh[j] = static_cast<int>(ls[j].size());
}
for (vcl_size_t j = 0; j < rgc[i].size(); j++)
{
(wvsg[static_cast<vcl_size_t>(lap[rgc[i][j]][0])])++;
(wvsh[static_cast<vcl_size_t>(lap[rgc[i][j]][1])])++;
}
k3 = 0;
k4 = 0;
for (vcl_size_t j = 0; j < ls.size(); j++)
{
if (wvsg[j] > wvs[j])
k3 = std::max(k3, wvsg[j]);
if (wvsh[j] > wvs[j])
k4 = std::max(k4, wvsh[j]);
}
if (k3 < k4 || (k3 == k4 && k1 <= k2) )
for (vcl_size_t j = 0; j < rgc[i].size(); j++)
ls[static_cast<vcl_size_t>(lap[rgc[i][j]][0])].push_back(rgc[i][j]);
else
for (vcl_size_t j = 0; j < rgc[i].size(); j++)
ls[static_cast<vcl_size_t>(lap[rgc[i][j]][1])].push_back(rgc[i][j]);
}
// renumber nodes in ls
rl.clear();
rl.resize(ls.size());
state = 1;
state_end = false;
while (!state_end)
{
switch (state)
{
case 1:
l = 0;
state = 4;
break;
case 2:
for (vcl_size_t i = 0; i < rl[l-1].size(); i++)
{
isn.assign(n, false);
for (std::map<int, double>::const_iterator it = matrix[static_cast<vcl_size_t>(rl[l-1][i])].begin();
it != matrix[static_cast<vcl_size_t>(rl[l-1][i])].end();
it++)
{
if (it->first == rl[l-1][i]) continue;
isn[static_cast<vcl_size_t>(it->first)] = true;
}
nodes.resize(0);
for (vcl_size_t j = 0; j < ls[l].size(); j++)
{
if (inr[static_cast<vcl_size_t>(ls[l][j])]) continue;
if (!isn[static_cast<vcl_size_t>(ls[l][j])]) continue;
tmp[0] = ls[l][j];
tmp[1] = static_cast<int>(matrix[static_cast<vcl_size_t>(ls[l][j])].size() - 1);
nodes.push_back(tmp);
}
std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
for (vcl_size_t j = 0; j < nodes.size(); j++)
{
rl[l].push_back(nodes[j][0]);
r[static_cast<vcl_size_t>(nodes[j][0])] = current_dof++;
inr[static_cast<vcl_size_t>(nodes[j][0])] = true;
}
}
case 3:
for (vcl_size_t i = 0; i < rl[l].size(); i++)
{
isn.assign(n, false);
for (std::map<int, double>::const_iterator it = matrix[static_cast<vcl_size_t>(rl[l][i])].begin();
it != matrix[static_cast<vcl_size_t>(rl[l][i])].end();
it++)
{
if (it->first == rl[l][i]) continue;
isn[static_cast<vcl_size_t>(it->first)] = true;
}
nodes.resize(0);
for (vcl_size_t j = 0; j < ls[l].size(); j++)
{
if (inr[static_cast<vcl_size_t>(ls[l][j])]) continue;
if (!isn[static_cast<vcl_size_t>(ls[l][j])]) continue;
tmp[0] = ls[l][j];
tmp[1] = static_cast<int>(matrix[static_cast<vcl_size_t>(ls[l][j])].size() - 1);
nodes.push_back(tmp);
}
std::sort(nodes.begin(), nodes.end(), detail::cuthill_mckee_comp_func);
for (vcl_size_t j = 0; j < nodes.size(); j++)
{
rl[l].push_back(nodes[j][0]);
r[static_cast<vcl_size_t>(nodes[j][0])] = current_dof++;
inr[static_cast<vcl_size_t>(nodes[j][0])] = true;
}
}
case 4:
if (rl[l].size() < ls[l].size())
{
deg_min = -1;
for (vcl_size_t j = 0; j < ls[l].size(); j++)
{
if (inr[static_cast<vcl_size_t>(ls[l][j])]) continue;
deg = static_cast<int>(matrix[static_cast<vcl_size_t>(ls[l][j])].size() - 1);
if (deg_min < 0 || deg < deg_min)
{
ind_min = ls[l][j];
deg_min = deg;
}
}
rl[l].push_back(ind_min);
r[static_cast<vcl_size_t>(ind_min)] = current_dof++;
inr[static_cast<vcl_size_t>(ind_min)] = true;
state = 3;
break;
}
case 5:
l++;
if (l < ls.size())
state = 2;
else
state_end = true;
break;
default:
break;
}
}
}
return r;
}
} //namespace viennacl
#endif
viennacl::reorder
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:367
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::detail::generate_layering
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
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::gibbs_poole_stockmeyer_tag
Tag class for identifying the Gibbs-Poole-Stockmeyer algorithm for reducing the bandwidth of a sparse...
Definition: gibbs_poole_stockmeyer.hpp:133
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::detail::gps_rg_components
std::vector< std::vector< int > > gps_rg_components(MatrixType const &matrix, int n, std::vector< int > const &rg)
Definition: gibbs_poole_stockmeyer.hpp:63
cuthill_mckee.hpp
Implementation of several flavors of the Cuthill-McKee algorithm. Experimental.
viennacl::detail::calc_layering_width
int calc_layering_width(std::vector< std::vector< int > > const &l)
Definition: gibbs_poole_stockmeyer.hpp:47
viennacl::detail::cuthill_mckee_comp_func
bool cuthill_mckee_comp_func(std::vector< int > const &a, std::vector< int > const &b)
Definition: cuthill_mckee.hpp:263