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