1 #ifndef VIENNACL_MISC_CUTHILL_MCKEE_HPP
2 #define VIENNACL_MISC_CUTHILL_MCKEE_HPP
47 template <
typename IndexT,
typename ValueT>
49 std::vector<bool> & dof_assigned_to_node,
50 std::vector<IndexT>
const & permutation)
54 for (
vcl_size_t i = 0; i < permutation.size(); i++)
56 if (!dof_assigned_to_node[i])
59 IndexT min_index =
static_cast<IndexT
>(
matrix.size());
61 for (
typename std::map<IndexT, ValueT>::const_iterator it =
matrix[i].begin(); it !=
matrix[i].end(); it++)
63 if (!dof_assigned_to_node[it->first])
66 if (permutation[it->first] > max_index)
67 max_index = permutation[it->first];
68 if (permutation[it->first] < min_index)
69 min_index = permutation[it->first];
71 if (max_index > min_index)
72 bw = std::max(bw, max_index - min_index);
88 template <
typename IndexT>
94 m =
static_cast<IndexT
>(comb.size());
97 while ( (k > 0) && ( ((k == m) && (comb[k-1] == static_cast<IndexT>(n)-1)) ||
98 ((k < m) && (comb[k-1] == comb[k] - 1) )) )
110 for (IndexT i = k; i < m; i++)
111 comb[i] = comb[k-1] + (i - k);
121 template <
typename MatrixT,
typename IndexT>
123 std::vector< std::vector<IndexT> > & layer_list)
125 std::vector<bool> node_visited_already(matrix.size(),
false);
130 for (
vcl_size_t i=0; i<layer_list.size(); ++i)
132 for (
typename std::vector<IndexT>::iterator it = layer_list[i].begin();
133 it != layer_list[i].end();
135 node_visited_already[*it] =
true;
141 while (layer_list.back().size() > 0)
144 layer_list.push_back(std::vector<IndexT>());
146 for (
typename std::vector<IndexT>::iterator it = layer_list[layer_index].begin();
147 it != layer_list[layer_index].end();
150 for (
typename MatrixT::value_type::const_iterator it2 = matrix[*it].begin();
151 it2 != matrix[*it].end();
154 if (it2->first == *it)
continue;
155 if (node_visited_already[it2->first])
continue;
157 layer_list.back().push_back(it2->first);
158 node_visited_already[it2->first] =
true;
164 layer_list.resize(layer_list.size()-1);
169 template <
typename MatrixType>
171 std::vector< std::vector<int> > & l,
176 std::vector<bool> inr(n,
false);
177 std::vector<int> nlist;
186 for (std::vector<int>::iterator it = l.back().begin();
187 it != l.back().end();
190 for (
typename MatrixType::value_type::const_iterator it2 = matrix[*it].begin();
191 it2 != matrix[*it].end();
194 if (it2->first == *it)
continue;
195 if (inr[it2->first])
continue;
197 nlist.push_back(it2->first);
198 inr[it2->first] =
true;
202 if (nlist.size() == 0)
214 template <
typename MatrixT,
typename IndexT>
216 std::vector<IndexT> & node_list)
218 std::vector<bool> node_visited_already(matrix.size(),
false);
219 std::deque<IndexT> node_queue;
224 for (
typename std::vector<IndexT>::iterator it = node_list.begin();
225 it != node_list.end();
228 node_queue.push_back(*it);
235 while (!node_queue.empty())
237 IndexT node_id = node_queue.front();
238 node_queue.pop_front();
240 if (!node_visited_already[node_id])
242 node_list.push_back(node_id);
243 node_visited_already[node_id] =
true;
245 for (
typename MatrixT::value_type::const_iterator it = matrix[node_id].begin();
246 it != matrix[node_id].end();
249 IndexT neighbor_node_id = it->first;
250 if (neighbor_node_id == node_id)
continue;
251 if (node_visited_already[neighbor_node_id])
continue;
253 node_queue.push_back(neighbor_node_id);
264 std::vector<int>
const & b)
266 return (a[1] < b[1]);
269 template <
typename IndexT>
271 std::pair<IndexT, IndexT>
const & b)
273 return (a.second < b.second);
286 template <
typename IndexT,
typename ValueT>
288 std::deque<IndexT> & node_assignment_queue,
289 std::vector<bool> & dof_assigned_to_node,
290 std::vector<IndexT> & permutation,
293 typedef std::pair<IndexT, IndexT> NodeIdDegreePair;
295 std::vector< NodeIdDegreePair > local_neighbor_nodes(
matrix.size());
297 while (!node_assignment_queue.empty())
300 vcl_size_t node_id = node_assignment_queue.front();
301 node_assignment_queue.pop_front();
304 if (!dof_assigned_to_node[node_id])
306 permutation[node_id] =
static_cast<IndexT
>(current_dof);
308 dof_assigned_to_node[node_id] =
true;
314 for (
typename std::map<IndexT, ValueT>::const_iterator neighbor_it =
matrix[node_id].begin();
315 neighbor_it !=
matrix[node_id].end();
318 if (!dof_assigned_to_node[neighbor_it->first])
320 local_neighbor_nodes[num_neighbors] = NodeIdDegreePair(neighbor_it->first, static_cast<IndexT>(
matrix[neighbor_it->first].size()));
326 std::sort(local_neighbor_nodes.begin(), local_neighbor_nodes.begin() + num_neighbors, detail::cuthill_mckee_comp_func_pair<IndexT>);
330 node_assignment_queue.push_back(local_neighbor_nodes[i].first);
363 template <
typename IndexT,
typename ValueT>
366 std::vector<IndexT> permutation(
matrix.size());
367 std::vector<bool> dof_assigned_to_node(
matrix.size(),
false);
368 std::deque<IndexT> node_assignment_queue;
372 while (current_dof <
matrix.size())
379 bool found_unassigned_node =
false;
382 if (!dof_assigned_to_node[i])
386 permutation[i] =
static_cast<IndexT
>(current_dof);
387 dof_assigned_to_node[i] =
true;
392 if (!found_unassigned_node)
394 current_min_degree =
matrix[i].size();
395 node_with_minimum_degree = i;
396 found_unassigned_node =
true;
401 current_min_degree =
matrix[i].size();
402 node_with_minimum_degree = i;
410 if (found_unassigned_node)
412 node_assignment_queue.push_back(static_cast<IndexT>(node_with_minimum_degree));
456 double starting_node_param_;
470 template <
typename IndexT,
typename ValueT>
471 std::vector<IndexT>
reorder(std::vector< std::map<IndexT, ValueT> >
const &
matrix,
477 std::vector<IndexT> permutation(n);
478 std::vector<bool> dof_assigned_to_node(n,
false);
479 std::vector<IndexT> nodes_in_strongly_connected_component;
480 std::vector<IndexT> parent_nodes;
485 std::vector<IndexT> comb;
487 nodes_in_strongly_connected_component.reserve(n);
488 parent_nodes.reserve(n);
493 while (current_dof <
matrix.size())
496 nodes_in_strongly_connected_component.resize(0);
499 if (!dof_assigned_to_node[i])
501 nodes_in_strongly_connected_component.push_back(static_cast<IndexT>(i));
510 for (
typename std::vector<IndexT>::iterator it = nodes_in_strongly_connected_component.begin();
511 it != nodes_in_strongly_connected_component.end();
515 if (deg_min == 0 || deg < deg_min)
517 if (deg_max == 0 || deg > deg_max)
520 deg_a = deg_min +
static_cast<vcl_size_t>(a * (deg_max - deg_min));
523 parent_nodes.resize(0);
524 for (
typename std::vector<IndexT>::iterator it = nodes_in_strongly_connected_component.begin();
525 it != nodes_in_strongly_connected_component.end();
529 parent_nodes.push_back(*it);
535 std::vector<bool> dof_assigned_to_node_backup = dof_assigned_to_node;
536 std::vector<bool> dof_assigned_to_node_best;
538 std::vector<IndexT> permutation_backup = permutation;
539 std::vector<IndexT> permutation_best = permutation;
555 dof_assigned_to_node = dof_assigned_to_node_backup;
556 permutation = permutation_backup;
557 current_dof = current_dof_backup;
559 std::deque<IndexT> node_queue;
562 for (
typename std::vector<IndexT>::iterator it = comb.begin(); it != comb.end(); it++)
563 node_queue.push_back(parent_nodes[*it]);
572 if (bw_best == 0 || bw < bw_best)
574 permutation_best = permutation;
576 dof_assigned_to_node_best = dof_assigned_to_node;
584 if ( (gmax > 0 && g > gmax) || g > parent_nodes.size())
589 comb[i] = static_cast<IndexT>(i);
596 permutation = permutation_best;
597 dof_assigned_to_node = dof_assigned_to_node_best;
advanced_cuthill_mckee_tag(double a=0.0, vcl_size_t gmax=1)
CTOR which may take the additional parameters for the advanced algorithm.
Definition: cuthill_mckee.hpp:447
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
void nodes_of_strongly_connected_component(MatrixT const &matrix, std::vector< IndexT > &node_list)
Fills the provided nodelist with all nodes of the same strongly connected component as the nodes in t...
Definition: cuthill_mckee.hpp:215
This file provides the forward declarations for the main types used within ViennaCL.
A dense matrix class.
Definition: forwards.h:293
void starting_node_param(double a)
Definition: cuthill_mckee.hpp:450
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
void max_root_nodes(vcl_size_t gmax)
Definition: cuthill_mckee.hpp:453
vcl_size_t max_root_nodes() const
Definition: cuthill_mckee.hpp:452
double starting_node_param() const
Definition: cuthill_mckee.hpp:449
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
vcl_size_t cuthill_mckee_on_strongly_connected_component(std::vector< std::map< IndexT, ValueT > > const &matrix, std::deque< IndexT > &node_assignment_queue, std::vector< bool > &dof_assigned_to_node, std::vector< IndexT > &permutation, vcl_size_t current_dof)
Runs the Cuthill-McKee algorithm on a strongly connected component of a graph.
Definition: cuthill_mckee.hpp:287
IndexT calc_reordered_bw(std::vector< std::map< IndexT, ValueT > > const &matrix, std::vector< bool > &dof_assigned_to_node, std::vector< IndexT > const &permutation)
Definition: cuthill_mckee.hpp:48
bool comb_inc(std::vector< IndexT > &comb, vcl_size_t n)
Definition: cuthill_mckee.hpp:89
bool cuthill_mckee_comp_func_pair(std::pair< IndexT, IndexT > const &a, std::pair< IndexT, IndexT > const &b)
Definition: cuthill_mckee.hpp:270
Tag for the advanced Cuthill-McKee algorithm (i.e. running the 'standard' Cuthill-McKee algorithm for...
Definition: cuthill_mckee.hpp:426
A tag class for selecting the Cuthill-McKee algorithm for reducing the bandwidth of a sparse matrix...
Definition: cuthill_mckee.hpp:347
bool cuthill_mckee_comp_func(std::vector< int > const &a, std::vector< int > const &b)
Definition: cuthill_mckee.hpp:263