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