Reference documentation for deal.II version GIT 58febcd5cf 2023-09-30 20:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
graph_coloring.h
Go to the documentation of this file.
1 
2 // ---------------------------------------------------------------------
3 //
4 // Copyright (C) 2013 - 2023 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE.md at
13 // the top level directory of deal.II.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef dealii_graph_coloring_h
18 # define dealii_graph_coloring_h
19 
20 
21 # include <deal.II/base/config.h>
22 
24 
25 # include <algorithm>
26 # include <functional>
27 # include <set>
28 # include <unordered_map>
29 # include <unordered_set>
30 # include <vector>
31 
32 
34 
35 class SparsityPattern;
36 
40 namespace GraphColoring
41 {
42  namespace internal
43  {
53  inline bool
55  const std::vector<types::global_dof_index> &indices1,
56  const std::vector<types::global_dof_index> &indices2)
57  {
58  // we assume that both arrays are sorted, so we can walk
59  // them in lockstep and see if we encounter an index that's
60  // in both arrays. once we reach the end of either array,
61  // we know that there is no intersection
62  std::vector<types::global_dof_index>::const_iterator p = indices1.begin(),
63  q = indices2.begin();
64  while ((p != indices1.end()) && (q != indices2.end()))
65  {
66  if (*p < *q)
67  ++p;
68  else if (*p > *q)
69  ++q;
70  else
71  // conflict found!
72  return true;
73  }
74 
75  // no conflict found!
76  return false;
77  }
78 
79 
108  template <typename Iterator>
109  std::vector<std::vector<Iterator>>
111  const Iterator &begin,
113  const std::function<std::vector<types::global_dof_index>(
114  const Iterator &)> &get_conflict_indices)
115  {
116  // Number of iterators.
117  unsigned int n_iterators = 0;
118 
119  // Create a map from conflict indices to iterators
120  std::unordered_map<types::global_dof_index, std::vector<Iterator>>
121  indices_to_iterators;
122  for (Iterator it = begin; it != end; ++it)
123  {
124  const std::vector<types::global_dof_index> conflict_indices =
125  get_conflict_indices(it);
126  const unsigned int n_conflict_indices = conflict_indices.size();
127  for (unsigned int i = 0; i < n_conflict_indices; ++i)
128  indices_to_iterators[conflict_indices[i]].push_back(it);
129  ++n_iterators;
130  }
131 
132  // create the very first zone which contains only the first
133  // iterator. then create the other zones. keep track of all the
134  // iterators that have already been assigned to a zone
135  std::vector<std::vector<Iterator>> zones(1,
136  std::vector<Iterator>(1, begin));
137  std::set<Iterator> used_it;
138  used_it.insert(begin);
139  while (used_it.size() != n_iterators)
140  {
141  // loop over the elements of the previous zone. for each element of
142  // the previous zone, get the conflict indices and from there get
143  // those iterators that are conflicting with the current element
144  typename std::vector<Iterator>::iterator previous_zone_it(
145  zones.back().begin());
146  typename std::vector<Iterator>::iterator previous_zone_end(
147  zones.back().end());
148  std::vector<Iterator> new_zone;
149  for (; previous_zone_it != previous_zone_end; ++previous_zone_it)
150  {
151  const std::vector<types::global_dof_index> conflict_indices =
152  get_conflict_indices(*previous_zone_it);
153 
154  const unsigned int n_conflict_indices(conflict_indices.size());
155  for (unsigned int i = 0; i < n_conflict_indices; ++i)
156  {
157  const std::vector<Iterator> &conflicting_elements =
158  indices_to_iterators[conflict_indices[i]];
159  for (unsigned int j = 0; j < conflicting_elements.size(); ++j)
160  {
161  // check that the iterator conflicting with the current
162  // one is not associated to a zone yet and if so, assign
163  // it to the current zone. mark it as used
164  //
165  // we can shortcut this test if the conflicting iterator
166  // is the current iterator
167  if ((conflicting_elements[j] != *previous_zone_it) &&
168  (used_it.count(conflicting_elements[j]) == 0))
169  {
170  new_zone.push_back(conflicting_elements[j]);
171  used_it.insert(conflicting_elements[j]);
172  }
173  }
174  }
175  }
176 
177  // If there are iterators in the new zone, then the zone is added to
178  // the partition. Otherwise, the graph is disconnected and we need to
179  // find an iterator on the other part of the graph. start the whole
180  // process again with the first iterator that hasn't been assigned to
181  // a zone yet
182  if (new_zone.size() != 0)
183  zones.push_back(new_zone);
184  else
185  for (Iterator it = begin; it != end; ++it)
186  if (used_it.count(it) == 0)
187  {
188  zones.push_back(std::vector<Iterator>(1, it));
189  used_it.insert(it);
190  break;
191  }
192  }
193 
194  return zones;
195  }
196 
197 
198 
221  template <typename Iterator>
222  void
224  std::vector<Iterator> &partition,
225  const std::function<std::vector<types::global_dof_index>(
226  const Iterator &)> &get_conflict_indices,
227  std::vector<std::vector<Iterator>> &partition_coloring)
228  {
229  partition_coloring.clear();
230 
231  // Number of zones composing the partitioning.
232  const unsigned int partition_size(partition.size());
233  std::vector<unsigned int> sorted_vertices(partition_size);
234  std::vector<int> degrees(partition_size);
235  std::vector<std::vector<types::global_dof_index>> conflict_indices(
236  partition_size);
237  std::vector<std::vector<unsigned int>> graph(partition_size);
238 
239  // Get the conflict indices associated to each iterator. The
240  // conflict_indices have to be sorted so we can more easily find conflicts
241  // later on
242  for (unsigned int i = 0; i < partition_size; ++i)
243  {
244  conflict_indices[i] = get_conflict_indices(partition[i]);
245  std::sort(conflict_indices[i].begin(), conflict_indices[i].end());
246  }
247 
248  // Compute the degree of each vertex of the graph using the
249  // intersection of the conflict indices.
250  for (unsigned int i = 0; i < partition_size; ++i)
251  for (unsigned int j = i + 1; j < partition_size; ++j)
252  // If the two iterators share indices then we increase the degree of
253  // the vertices and create an ''edge'' in the graph.
254  if (have_nonempty_intersection(conflict_indices[i],
255  conflict_indices[j]))
256  {
257  ++degrees[i];
258  ++degrees[j];
259  graph[i].push_back(j);
260  graph[j].push_back(i);
261  }
262 
263  // Sort the vertices by decreasing degree.
264  std::vector<int>::iterator degrees_it;
265  for (unsigned int i = 0; i < partition_size; ++i)
266  {
267  // Find the largest element.
268  degrees_it = std::max_element(degrees.begin(), degrees.end());
269  sorted_vertices[i] = degrees_it - degrees.begin();
270  // Put the largest element to -1 so it cannot be chosen again.
271  *degrees_it = -1;
272  }
273 
274  // Color the graph.
275  std::vector<std::unordered_set<unsigned int>> colors_used;
276  for (unsigned int i = 0; i < partition_size; ++i)
277  {
278  const unsigned int current_vertex(sorted_vertices[i]);
279  bool new_color(true);
280  // Try to use an existing color, i.e., try to find a color which is
281  // not associated to one of the vertices linked to current_vertex.
282  // Loop over the color.
283  for (unsigned int j = 0; j < partition_coloring.size(); ++j)
284  {
285  // Loop on the vertices linked to current_vertex. If one vertex
286  // linked to current_vertex is already using the color j, this
287  // color cannot be used anymore.
288  bool unused_color(true);
289  for (const auto adjacent_vertex : graph[current_vertex])
290  if (colors_used[j].count(adjacent_vertex) == 1)
291  {
292  unused_color = false;
293  break;
294  }
295  if (unused_color)
296  {
297  partition_coloring[j].push_back(partition[current_vertex]);
298  colors_used[j].insert(current_vertex);
299  new_color = false;
300  break;
301  }
302  }
303  // Add a new color.
304  if (new_color)
305  {
306  partition_coloring.push_back(
307  std::vector<Iterator>(1, partition[current_vertex]));
308  std::unordered_set<unsigned int> tmp;
309  tmp.insert(current_vertex);
310  colors_used.push_back(tmp);
311  }
312  }
313  }
314 
315 
316 
326  template <typename Iterator>
327  std::vector<std::vector<Iterator>>
329  const std::vector<std::vector<std::vector<Iterator>>> &partition_coloring)
330  {
331  std::vector<std::vector<Iterator>> coloring;
332 
333  // Count the number of iterators in each color.
334  const unsigned int partition_size(partition_coloring.size());
335  std::vector<std::vector<unsigned int>> colors_counter(partition_size);
336  for (unsigned int i = 0; i < partition_size; ++i)
337  {
338  const unsigned int n_colors(partition_coloring[i].size());
339  colors_counter[i].resize(n_colors);
340  for (unsigned int j = 0; j < n_colors; ++j)
341  colors_counter[i][j] = partition_coloring[i][j].size();
342  }
343 
344  // Find the partition with the largest number of colors for the even
345  // partition.
346  unsigned int i_color(0);
347  unsigned int max_even_n_colors(0);
348  const unsigned int colors_size(colors_counter.size());
349  for (unsigned int i = 0; i < colors_size; i += 2)
350  {
351  if (max_even_n_colors < colors_counter[i].size())
352  {
353  max_even_n_colors = colors_counter[i].size();
354  i_color = i;
355  }
356  }
357  coloring.resize(max_even_n_colors);
358  for (unsigned int j = 0; j < colors_counter[i_color].size(); ++j)
359  coloring[j] = partition_coloring[i_color][j];
360 
361  for (unsigned int i = 0; i < partition_size; i += 2)
362  {
363  if (i != i_color)
364  {
365  std::unordered_set<unsigned int> used_k;
366  for (unsigned int j = 0; j < colors_counter[i].size(); ++j)
367  {
368  // Find the color in the current partition with the largest
369  // number of iterators.
370  std::vector<unsigned int>::iterator it;
371  it = std::max_element(colors_counter[i].begin(),
372  colors_counter[i].end());
373  unsigned int min_iterators(static_cast<unsigned int>(-1));
374  unsigned int pos(0);
375  // Find the color of coloring with the least number of colors
376  // among the colors that have not been used yet.
377  for (unsigned int k = 0; k < max_even_n_colors; ++k)
378  if (used_k.count(k) == 0)
379  if (colors_counter[i_color][k] < min_iterators)
380  {
381  min_iterators = colors_counter[i_color][k];
382  pos = k;
383  }
384  colors_counter[i_color][pos] += *it;
385  // Concatenate the current color with the existing coloring.
386  coloring[pos].insert(
387  coloring[pos].end(),
388  partition_coloring[i][it - colors_counter[i].begin()]
389  .begin(),
390  partition_coloring[i][it - colors_counter[i].begin()]
391  .end());
392  used_k.insert(pos);
393  // Put the number of iterators to the current color to zero.
394  *it = 0;
395  }
396  }
397  }
398 
399  // If there is more than one partition, do the same thing that we did for
400  // the even partitions to the odd partitions
401  if (partition_size > 1)
402  {
403  unsigned int max_odd_n_colors(0);
404  for (unsigned int i = 1; i < partition_size; i += 2)
405  {
406  if (max_odd_n_colors < colors_counter[i].size())
407  {
408  max_odd_n_colors = colors_counter[i].size();
409  i_color = i;
410  }
411  }
412  coloring.resize(max_even_n_colors + max_odd_n_colors);
413  for (unsigned int j = 0; j < colors_counter[i_color].size(); ++j)
414  coloring[max_even_n_colors + j] = partition_coloring[i_color][j];
415 
416  for (unsigned int i = 1; i < partition_size; i += 2)
417  {
418  if (i != i_color)
419  {
420  std::unordered_set<unsigned int> used_k;
421  for (unsigned int j = 0; j < colors_counter[i].size(); ++j)
422  {
423  // Find the color in the current partition with the
424  // largest number of iterators.
425  std::vector<unsigned int>::iterator it;
426  it = std::max_element(colors_counter[i].begin(),
427  colors_counter[i].end());
428  unsigned int min_iterators(static_cast<unsigned int>(-1));
429  unsigned int pos(0);
430  // Find the color of coloring with the least number of
431  // colors among the colors that have not been used yet.
432  for (unsigned int k = 0; k < max_odd_n_colors; ++k)
433  if (used_k.count(k) == 0)
434  if (colors_counter[i_color][k] < min_iterators)
435  {
436  min_iterators = colors_counter[i_color][k];
437  pos = k;
438  }
439  colors_counter[i_color][pos] += *it;
440  // Concatenate the current color with the existing
441  // coloring.
442  coloring[max_even_n_colors + pos].insert(
443  coloring[max_even_n_colors + pos].end(),
444  partition_coloring[i][it - colors_counter[i].begin()]
445  .begin(),
446  partition_coloring[i][it - colors_counter[i].begin()]
447  .end());
448  used_k.insert(pos);
449  // Put the number of iterators to the current color to
450  // zero.
451  *it = 0;
452  }
453  }
454  }
455  }
456 
457  return coloring;
458  }
459  } // namespace internal
460 
461 
539  template <typename Iterator>
540  std::vector<std::vector<Iterator>>
542  const Iterator &begin,
544  const std::function<std::vector<types::global_dof_index>(
545  const std_cxx20::type_identity_t<Iterator> &)> &get_conflict_indices)
546  {
547  Assert(begin != end,
548  ExcMessage(
549  "GraphColoring is not prepared to deal with empty ranges!"));
550 
551  // Create the partitioning.
552  std::vector<std::vector<Iterator>> partitioning =
553  internal::create_partitioning(begin, end, get_conflict_indices);
554 
555  // Color the iterators within each partition.
556  // Run the coloring algorithm on each zone in parallel
557  const unsigned int partitioning_size(partitioning.size());
558  std::vector<std::vector<std::vector<Iterator>>> partition_coloring(
559  partitioning_size);
560 
561  Threads::TaskGroup<> tasks;
562  for (unsigned int i = 0; i < partitioning_size; ++i)
563  tasks += Threads::new_task(&internal::make_dsatur_coloring<Iterator>,
564  partitioning[i],
565  get_conflict_indices,
566  partition_coloring[i]);
567  tasks.join_all();
568 
569  // Gather the colors together.
570  return internal::gather_colors(partition_coloring);
571  }
572 
579  unsigned int
580  color_sparsity_pattern(const SparsityPattern &sparsity_pattern,
581  std::vector<unsigned int> &color_indices);
582 
583 } // namespace GraphColoring
584 
586 
587 
588 //---------------------------- graph_coloring.h ---------------------------
589 // end of #ifndef dealii_graph_coloring_h
590 #endif
591 //---------------------------- graph_coloring.h ---------------------------
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcMessage(std::string arg1)
Task< RT > new_task(const std::function< RT()> &function)
std::vector< std::vector< Iterator > > create_partitioning(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const std::function< std::vector< types::global_dof_index >(const Iterator &)> &get_conflict_indices)
bool have_nonempty_intersection(const std::vector< types::global_dof_index > &indices1, const std::vector< types::global_dof_index > &indices2)
void make_dsatur_coloring(std::vector< Iterator > &partition, const std::function< std::vector< types::global_dof_index >(const Iterator &)> &get_conflict_indices, std::vector< std::vector< Iterator >> &partition_coloring)
std::vector< std::vector< Iterator > > gather_colors(const std::vector< std::vector< std::vector< Iterator >>> &partition_coloring)
unsigned int color_sparsity_pattern(const SparsityPattern &sparsity_pattern, std::vector< unsigned int > &color_indices)
std::vector< std::vector< Iterator > > make_graph_coloring(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const std::function< std::vector< types::global_dof_index >(const std_cxx20::type_identity_t< Iterator > &)> &get_conflict_indices)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
typename type_identity< T >::type type_identity_t
Definition: type_traits.h:96