deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20: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
sparsity_tools.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
21
22#include <algorithm>
23#include <functional>
24#include <memory>
25#include <set>
26
27#ifdef DEAL_II_WITH_MPI
28# include <deal.II/base/mpi.h>
30
33#endif
34
35#ifdef DEAL_II_WITH_METIS
36extern "C"
37{
38# include <metis.h>
39}
40#endif
41
42#ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
43# include <zoltan_cpp.h>
44#endif
45
46#include <string>
47
48
50
51namespace SparsityTools
52{
53 namespace
54 {
55 void
56 partition_metis(const SparsityPattern &sparsity_pattern,
57 const std::vector<unsigned int> &cell_weights,
58 const unsigned int n_partitions,
59 std::vector<unsigned int> &partition_indices)
60 {
61 // Make sure that METIS is actually
62 // installed and detected
63#ifndef DEAL_II_WITH_METIS
64 (void)sparsity_pattern;
65 (void)cell_weights;
66 (void)n_partitions;
67 (void)partition_indices;
69#else
70
71 // Generate the data structures for METIS. Note that this is particularly
72 // simple, since METIS wants exactly our compressed row storage format.
73 // We only have to set up a few auxiliary arrays and convert from our
74 // unsigned cell weights to signed ones.
75 idx_t n = static_cast<signed int>(sparsity_pattern.n_rows());
76
77 idx_t ncon = 1; // number of balancing constraints (should be >0)
78
79 // We can not partition n items into more than n parts. METIS will
80 // generate non-sensical output (everything is owned by a single process)
81 // and complain with a message (but won't return an error code!):
82 // ***Cannot bisect a graph with 0 vertices!
83 // ***You are trying to partition a graph into too many parts!
84 idx_t nparts =
85 std::min(n,
86 static_cast<idx_t>(
87 n_partitions)); // number of subdomains to create
88
89 // use default options for METIS
90 idx_t options[METIS_NOPTIONS];
91 METIS_SetDefaultOptions(options);
92
93 // one more nuisance: we have to copy our own data to arrays that store
94 // signed integers :-(
95 std::vector<idx_t> int_rowstart(1);
96 int_rowstart.reserve(sparsity_pattern.n_rows() + 1);
97 std::vector<idx_t> int_colnums;
98 int_colnums.reserve(sparsity_pattern.n_nonzero_elements());
99 for (SparsityPattern::size_type row = 0; row < sparsity_pattern.n_rows();
100 ++row)
101 {
102 for (SparsityPattern::iterator col = sparsity_pattern.begin(row);
103 col < sparsity_pattern.end(row);
104 ++col)
105 int_colnums.push_back(col->column());
106 int_rowstart.push_back(int_colnums.size());
107 }
108
109 std::vector<idx_t> int_partition_indices(sparsity_pattern.n_rows());
110
111 // Set up cell weighting option
112 std::vector<idx_t> int_cell_weights;
113 if (cell_weights.size() > 0)
114 {
115 Assert(cell_weights.size() == sparsity_pattern.n_rows(),
116 ExcDimensionMismatch(cell_weights.size(),
117 sparsity_pattern.n_rows()));
118 int_cell_weights.resize(cell_weights.size());
119 std::copy(cell_weights.begin(),
120 cell_weights.end(),
121 int_cell_weights.begin());
122 }
123 // Set a pointer to the optional cell weighting information.
124 // METIS expects a null pointer if there are no weights to be considered.
125 idx_t *const p_int_cell_weights =
126 (cell_weights.size() > 0 ? int_cell_weights.data() : nullptr);
127
128
129 // Make use of METIS' error code.
130 int ierr;
131
132 // Select which type of partitioning to create
133
134 // Use recursive if the number of partitions is less than or equal to 8
135 idx_t dummy; // output: # of edges cut by the resulting partition
136 if (nparts <= 8)
137 ierr = METIS_PartGraphRecursive(&n,
138 &ncon,
139 int_rowstart.data(),
140 int_colnums.data(),
141 p_int_cell_weights,
142 nullptr,
143 nullptr,
144 &nparts,
145 nullptr,
146 nullptr,
147 options,
148 &dummy,
149 int_partition_indices.data());
150
151 // Otherwise use kway
152 else
153 ierr = METIS_PartGraphKway(&n,
154 &ncon,
155 int_rowstart.data(),
156 int_colnums.data(),
157 p_int_cell_weights,
158 nullptr,
159 nullptr,
160 &nparts,
161 nullptr,
162 nullptr,
163 options,
164 &dummy,
165 int_partition_indices.data());
166
167 // If metis returns normally, an error code METIS_OK=1 is returned from
168 // the above functions (see metish.h)
169 AssertThrow(ierr == 1, ExcMETISError(ierr));
170
171 // now copy back generated indices into the output array
172 std::copy(int_partition_indices.begin(),
173 int_partition_indices.end(),
174 partition_indices.begin());
175#endif
176 }
177
178
179// Query functions unused if zoltan is not installed
180#ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
181 // Query functions for partition_zoltan
182 int
183 get_number_of_objects(void *data, int *ierr)
184 {
185 SparsityPattern *graph = reinterpret_cast<SparsityPattern *>(data);
186
187 *ierr = ZOLTAN_OK;
188
189 return graph->n_rows();
190 }
191
192
193 void
194 get_object_list(void *data,
195 int /*sizeGID*/,
196 int /*sizeLID*/,
197 ZOLTAN_ID_PTR globalID,
198 ZOLTAN_ID_PTR localID,
199 int /*wgt_dim*/,
200 float * /*obj_wgts*/,
201 int *ierr)
202 {
203 SparsityPattern *graph = reinterpret_cast<SparsityPattern *>(data);
204 *ierr = ZOLTAN_OK;
205
206 Assert(globalID != nullptr, ExcInternalError());
207 Assert(localID != nullptr, ExcInternalError());
208
209 // set global degrees of freedom
210 auto n_dofs = graph->n_rows();
211
212 for (unsigned int i = 0; i < n_dofs; ++i)
213 {
214 globalID[i] = i;
215 localID[i] = i; // Same as global ids.
216 }
217 }
218
219
220 void
221 get_num_edges_list(void *data,
222 int /*sizeGID*/,
223 int /*sizeLID*/,
224 int num_obj,
225 ZOLTAN_ID_PTR globalID,
226 ZOLTAN_ID_PTR /*localID*/,
227 int *numEdges,
228 int *ierr)
229 {
230 SparsityPattern *graph = reinterpret_cast<SparsityPattern *>(data);
231
232 *ierr = ZOLTAN_OK;
233
234 Assert(numEdges != nullptr, ExcInternalError());
235
236 for (int i = 0; i < num_obj; ++i)
237 {
238 if (graph->exists(i, i)) // Check if diagonal element is present
239 numEdges[i] = graph->row_length(globalID[i]) - 1;
240 else
241 numEdges[i] = graph->row_length(globalID[i]);
242 }
243 }
244
245
246
247 void
248 get_edge_list(void *data,
249 int /*sizeGID*/,
250 int /*sizeLID*/,
251 int num_obj,
252 ZOLTAN_ID_PTR /*globalID*/,
253 ZOLTAN_ID_PTR /*localID*/,
254 int * /*num_edges*/,
255 ZOLTAN_ID_PTR nborGID,
256 int *nborProc,
257 int /*wgt_dim*/,
258 float * /*ewgts*/,
259 int *ierr)
260 {
261 SparsityPattern *graph = reinterpret_cast<SparsityPattern *>(data);
262 *ierr = ZOLTAN_OK;
263
264 ZOLTAN_ID_PTR nextNborGID = nborGID;
265 int *nextNborProc = nborProc;
266
267 // Loop through rows corresponding to indices in globalID implicitly
269 i < static_cast<SparsityPattern::size_type>(num_obj);
270 ++i)
271 {
272 // Loop through each column to find neighbours
273 for (SparsityPattern::iterator col = graph->begin(i);
274 col < graph->end(i);
275 ++col)
276 // Ignore diagonal entries. Not needed for partitioning.
277 if (i != col->column())
278 {
279 Assert(nextNborGID != nullptr, ExcInternalError());
280 Assert(nextNborProc != nullptr, ExcInternalError());
281
282 *nextNborGID++ = col->column();
283 *nextNborProc++ = 0; // All the vertices on processor 0
284 }
285 }
286 }
287#endif
288
289
290 void
291 partition_zoltan(const SparsityPattern &sparsity_pattern,
292 const std::vector<unsigned int> &cell_weights,
293 const unsigned int n_partitions,
294 std::vector<unsigned int> &partition_indices)
295 {
296 // Make sure that ZOLTAN is actually
297 // installed and detected
298#ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
299 (void)sparsity_pattern;
300 (void)cell_weights;
301 (void)n_partitions;
302 (void)partition_indices;
304#else
305
306 Assert(
307 cell_weights.empty(),
309 "The cell weighting functionality for Zoltan has not yet been implemented."));
310 (void)cell_weights;
311
312 // MPI environment must have been initialized by this point.
313 std::unique_ptr<Zoltan> zz = std::make_unique<Zoltan>(MPI_COMM_SELF);
314
315 // General parameters
316 // DEBUG_LEVEL call must precede the call to LB_METHOD
317 zz->Set_Param("DEBUG_LEVEL", "0"); // set level of debug info
318 zz->Set_Param(
319 "LB_METHOD",
320 "GRAPH"); // graph based partition method (LB-load balancing)
321 zz->Set_Param("NUM_LOCAL_PARTS",
322 std::to_string(n_partitions)); // set number of partitions
323
324 // The PHG partitioner is a hypergraph partitioner that Zoltan could use
325 // for graph partitioning.
326 // If number of vertices in hyperedge divided by total vertices in
327 // hypergraph exceeds PHG_EDGE_SIZE_THRESHOLD,
328 // then the hyperedge will be omitted as such (dense) edges will likely
329 // incur high communication costs regardless of the partition.
330 // PHG_EDGE_SIZE_THRESHOLD value is raised to 0.5 from the default
331 // value of 0.25 so that the PHG partitioner doesn't throw warning saying
332 // "PHG_EDGE_SIZE_THRESHOLD is low ..." after removing all dense edges.
333 // For instance, in two dimensions if the triangulation being partitioned
334 // is two quadrilaterals sharing an edge and if PHG_EDGE_SIZE_THRESHOLD
335 // value is set to 0.25, PHG will remove all the edges throwing the
336 // above warning.
337 zz->Set_Param("PHG_EDGE_SIZE_THRESHOLD", "0.5");
338
339 // Need a non-const object equal to sparsity_pattern
340 SparsityPattern graph;
341 graph.copy_from(sparsity_pattern);
342
343 // Set query functions
344 zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
345 zz->Set_Obj_List_Fn(get_object_list, &graph);
346 zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
347 zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
348
349 // Variables needed by partition function
350 int changes = 0;
351 int num_gid_entries = 1;
352 int num_lid_entries = 1;
353 int num_import = 0;
354 ZOLTAN_ID_PTR import_global_ids = nullptr;
355 ZOLTAN_ID_PTR import_local_ids = nullptr;
356 int *import_procs = nullptr;
357 int *import_to_part = nullptr;
358 int num_export = 0;
359 ZOLTAN_ID_PTR export_global_ids = nullptr;
360 ZOLTAN_ID_PTR export_local_ids = nullptr;
361 int *export_procs = nullptr;
362 int *export_to_part = nullptr;
363
364 // call partitioner
365 const int rc = zz->LB_Partition(changes,
366 num_gid_entries,
367 num_lid_entries,
368 num_import,
369 import_global_ids,
370 import_local_ids,
371 import_procs,
372 import_to_part,
373 num_export,
374 export_global_ids,
375 export_local_ids,
376 export_procs,
377 export_to_part);
378 (void)rc;
379
380 // check for error code in partitioner
381 Assert(rc == ZOLTAN_OK, ExcInternalError());
382
383 // By default, all indices belong to part 0. After zoltan partition
384 // some are migrated to different part ID, which is stored in
385 // export_to_part array.
386 std::fill(partition_indices.begin(), partition_indices.end(), 0);
387
388 // copy from export_to_part to partition_indices, whose part_ids != 0.
389 Assert(export_to_part != nullptr, ExcInternalError());
390 for (int i = 0; i < num_export; ++i)
391 partition_indices[export_local_ids[i]] = export_to_part[i];
392#endif
393 }
394 } // namespace
395
396
397 void
398 partition(const SparsityPattern &sparsity_pattern,
399 const unsigned int n_partitions,
400 std::vector<unsigned int> &partition_indices,
401 const Partitioner partitioner)
402 {
403 std::vector<unsigned int> cell_weights;
404
405 // Call the other more general function
406 partition(sparsity_pattern,
407 cell_weights,
408 n_partitions,
409 partition_indices,
410 partitioner);
411 }
412
413
414 void
415 partition(const SparsityPattern &sparsity_pattern,
416 const std::vector<unsigned int> &cell_weights,
417 const unsigned int n_partitions,
418 std::vector<unsigned int> &partition_indices,
419 const Partitioner partitioner)
420 {
421 Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
423 Assert(sparsity_pattern.is_compressed(),
425
426 Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
427 Assert(partition_indices.size() == sparsity_pattern.n_rows(),
428 ExcInvalidArraySize(partition_indices.size(),
429 sparsity_pattern.n_rows()));
430
431 // check for an easy return
432 if (n_partitions == 1 || (sparsity_pattern.n_rows() == 1))
433 {
434 std::fill_n(partition_indices.begin(), partition_indices.size(), 0U);
435 return;
436 }
437
438 if (partitioner == Partitioner::metis)
439 partition_metis(sparsity_pattern,
440 cell_weights,
441 n_partitions,
442 partition_indices);
443 else if (partitioner == Partitioner::zoltan)
444 partition_zoltan(sparsity_pattern,
445 cell_weights,
446 n_partitions,
447 partition_indices);
448 else
450 }
451
452
453 unsigned int
455 std::vector<unsigned int> &color_indices)
456 {
457 // Make sure that ZOLTAN is actually
458 // installed and detected
459#ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
460 (void)sparsity_pattern;
461 (void)color_indices;
463 return 0;
464#else
465 // coloring algorithm is run in serial by each processor.
466 std::unique_ptr<Zoltan> zz = std::make_unique<Zoltan>(MPI_COMM_SELF);
467
468 // Coloring parameters
469 // DEBUG_LEVEL must precede all other calls
470 zz->Set_Param("DEBUG_LEVEL", "0"); // level of debug info
471 zz->Set_Param("COLORING_PROBLEM", "DISTANCE-1"); // Standard coloring
472 zz->Set_Param("NUM_GID_ENTRIES", "1"); // 1 entry represents global ID
473 zz->Set_Param("NUM_LID_ENTRIES", "1"); // 1 entry represents local ID
474 zz->Set_Param("OBJ_WEIGHT_DIM", "0"); // object weights not used
475 zz->Set_Param("RECOLORING_NUM_OF_ITERATIONS", "0");
476
477 // Zoltan::Color function requires a non-const SparsityPattern object
478 SparsityPattern graph;
479 graph.copy_from(sparsity_pattern);
480
481 // Set query functions required by coloring function
482 zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
483 zz->Set_Obj_List_Fn(get_object_list, &graph);
484 zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
485 zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
486
487 // Variables needed by coloring function
488 int num_gid_entries = 1;
489 const int num_objects = graph.n_rows();
490
491 // Preallocate input variables. Element type fixed by ZOLTAN.
492 std::vector<ZOLTAN_ID_TYPE> global_ids(num_objects);
493 std::vector<int> color_exp(num_objects);
494
495 // Set ids for which coloring needs to be done
496 for (int i = 0; i < num_objects; ++i)
497 global_ids[i] = i;
498
499 // Call ZOLTAN coloring algorithm
500 int rc = zz->Color(num_gid_entries,
501 num_objects,
502 global_ids.data(),
503 color_exp.data());
504
505 (void)rc;
506 // Check for error code
507 Assert(rc == ZOLTAN_OK, ExcInternalError());
508
509 // Allocate and assign color indices
510 color_indices.resize(num_objects);
511 Assert(color_exp.size() == color_indices.size(),
512 ExcDimensionMismatch(color_exp.size(), color_indices.size()));
513
514 std::copy(color_exp.begin(), color_exp.end(), color_indices.begin());
515
516 unsigned int n_colors =
517 *(std::max_element(color_indices.begin(), color_indices.end()));
518 return n_colors;
519#endif
520 }
521
522
523 namespace internal
524 {
532 const DynamicSparsityPattern &sparsity,
533 const std::vector<DynamicSparsityPattern::size_type> &new_indices)
534 {
535 DynamicSparsityPattern::size_type starting_point =
537 DynamicSparsityPattern::size_type min_coordination = sparsity.n_rows();
538 for (DynamicSparsityPattern::size_type row = 0; row < sparsity.n_rows();
539 ++row)
540 // look over all as-yet unnumbered indices
541 if (new_indices[row] == numbers::invalid_size_type)
542 {
543 if (sparsity.row_length(row) < min_coordination)
544 {
545 min_coordination = sparsity.row_length(row);
546 starting_point = row;
547 }
548 }
549
550 // now we still have to care for the case that no unnumbered dof has a
551 // coordination number less than sparsity.n_rows(). this rather exotic
552 // case only happens if we only have one cell, as far as I can see,
553 // but there may be others as well.
554 //
555 // if that should be the case, we can chose an arbitrary dof as
556 // starting point, e.g. the first unnumbered one
557 if (starting_point == numbers::invalid_size_type)
558 {
559 for (DynamicSparsityPattern::size_type i = 0; i < new_indices.size();
560 ++i)
561 if (new_indices[i] == numbers::invalid_size_type)
562 {
563 starting_point = i;
564 break;
565 }
566
567 Assert(starting_point != numbers::invalid_size_type,
569 }
570
571 return starting_point;
572 }
573 } // namespace internal
574
575
576
577 void
579 const DynamicSparsityPattern &sparsity,
580 std::vector<DynamicSparsityPattern::size_type> &new_indices,
581 const std::vector<DynamicSparsityPattern::size_type> &starting_indices)
582 {
583 Assert(sparsity.n_rows() == sparsity.n_cols(),
584 ExcDimensionMismatch(sparsity.n_rows(), sparsity.n_cols()));
585 Assert(sparsity.n_rows() == new_indices.size(),
586 ExcDimensionMismatch(sparsity.n_rows(), new_indices.size()));
587 Assert(starting_indices.size() <= sparsity.n_rows(),
589 "You can't specify more starting indices than there are rows"));
590 Assert(sparsity.row_index_set().size() == 0 ||
591 sparsity.row_index_set().size() == sparsity.n_rows(),
593 "Only valid for sparsity patterns which store all rows."));
594 for (const auto starting_index : starting_indices)
595 {
596 (void)starting_index;
597 Assert(starting_index < sparsity.n_rows(),
598 ExcMessage("Invalid starting index: All starting indices need "
599 "to be between zero and the number of rows in the "
600 "sparsity pattern."));
601 }
602
603 // store the indices of the dofs renumbered in the last round. Default to
604 // starting points
605 std::vector<DynamicSparsityPattern::size_type> last_round_dofs(
606 starting_indices);
607
608 // initialize the new_indices array with invalid values
609 std::fill(new_indices.begin(),
610 new_indices.end(),
612
613 // if no starting indices were given: find dof with lowest coordination
614 // number
615 if (last_round_dofs.empty())
616 last_round_dofs.push_back(
617 internal::find_unnumbered_starting_index(sparsity, new_indices));
618
619 // store next free dof index
620 DynamicSparsityPattern::size_type next_free_number = 0;
621
622 // enumerate the first round dofs
623 for (const auto &last_round_dof : last_round_dofs)
624 new_indices[last_round_dof] = next_free_number++;
625
626 // store the indices of the dofs to be renumbered in the next round
627 std::vector<DynamicSparsityPattern::size_type> next_round_dofs;
628
629 // store for each coordination number the dofs with these coordination
630 // number
631 std::vector<std::pair<unsigned int, DynamicSparsityPattern::size_type>>
632 dofs_by_coordination;
633
634 // now do as many steps as needed to renumber all dofs
635 while (true)
636 {
637 next_round_dofs.clear();
638
639 // find all neighbors of the dofs numbered in the last round
640 for (const auto dof : last_round_dofs)
641 {
642 const unsigned int row_length = sparsity.row_length(dof);
643 for (unsigned int i = 0; i < row_length; ++i)
644 {
645 // skip dofs which are already numbered
646 const auto column = sparsity.column_number(dof, i);
647 if (new_indices[column] == numbers::invalid_size_type)
648 {
649 next_round_dofs.push_back(column);
650
651 // assign a dummy value to 'new_indices' to avoid adding
652 // the same index again; those will get the right number
653 // at the end of the outer 'while' loop
654 new_indices[column] = 0;
655 }
656 }
657 }
658
659 // check whether there are any new dofs in the list. if there are
660 // none, then we have completely numbered the current component of the
661 // graph. check if there are as yet unnumbered components of the graph
662 // that we would then have to do next
663 if (next_round_dofs.empty())
664 {
665 if (std::find(new_indices.begin(),
666 new_indices.end(),
667 numbers::invalid_size_type) == new_indices.end())
668 // no unnumbered indices, so we can leave now
669 break;
670
671 // otherwise find a valid starting point for the next component of
672 // the graph and continue with numbering that one. we only do so
673 // if no starting indices were provided by the user (see the
674 // documentation of this function) so produce an error if we got
675 // here and starting indices were given
676 Assert(starting_indices.empty(),
677 ExcMessage("The input graph appears to have more than one "
678 "component, but as stated in the documentation "
679 "we only want to reorder such graphs if no "
680 "starting indices are given. The function was "
681 "called with starting indices, however."));
682
683 next_round_dofs.push_back(
684 internal::find_unnumbered_starting_index(sparsity, new_indices));
685 }
686
687
688 // find coordination number for each of these dofs
689 dofs_by_coordination.clear();
690 for (const types::global_dof_index next_round_dof : next_round_dofs)
691 dofs_by_coordination.emplace_back(sparsity.row_length(next_round_dof),
692 next_round_dof);
693 std::sort(dofs_by_coordination.begin(), dofs_by_coordination.end());
694
695 // assign new DoF numbers to the elements of the present front:
696 for (const auto &i : dofs_by_coordination)
697 new_indices[i.second] = next_free_number++;
698
699 // after that: use this round's dofs for the next round
700 last_round_dofs.swap(next_round_dofs);
701 }
702
703 // test for all indices numbered. this mostly tests whether the
704 // front-marching-algorithm (which Cuthill-McKee actually is) has reached
705 // all points.
706 Assert((std::find(new_indices.begin(),
707 new_indices.end(),
708 numbers::invalid_size_type) == new_indices.end()) &&
709 (next_free_number == sparsity.n_rows()),
711 }
712
713
714
715 namespace internal
716 {
717 void
719 const DynamicSparsityPattern &connectivity,
720 std::vector<DynamicSparsityPattern::size_type> &renumbering)
721 {
722 AssertDimension(connectivity.n_rows(), connectivity.n_cols());
723 AssertDimension(connectivity.n_rows(), renumbering.size());
724 Assert(connectivity.row_index_set().size() == 0 ||
725 connectivity.row_index_set().size() == connectivity.n_rows(),
727 "Only valid for sparsity patterns which store all rows."));
728
729 // The algorithm below works by partitioning the rows in the
730 // connectivity graph, called nodes, into groups. The groups are defined
731 // as those nodes in immediate neighborhood of some pivot node, which we
732 // choose by minimal adjacency below.
733
734 // We define two types of node categories for nodes not yet classified,
735 // one consisting of all nodes we've not seen at all, and one for nodes
736 // identified as neighbors (variable current_neighbors below) but not
737 // yet grouped. We use this classification in combination with an
738 // unsorted vector, which is much faster than keeping a sorted data
739 // structure (e.g. std::set)
740 constexpr types::global_dof_index unseen_node =
742 constexpr types::global_dof_index available_node = unseen_node - 1;
743 const types::global_dof_index n_nodes = connectivity.n_rows();
744 std::vector<types::global_dof_index> touched_nodes(n_nodes, unseen_node);
745
746 std::vector<unsigned int> row_lengths(n_nodes);
747 std::vector<types::global_dof_index> current_neighbors;
748 std::vector<types::global_dof_index> group_starts(1);
749 std::vector<types::global_dof_index> group_indices;
750 group_indices.reserve(n_nodes);
751
752 // First collect the number of neighbors for each node. We use this
753 // field to find the next node with the minimum number of non-touched
754 // neighbors in the field n_remaining_neighbors, so we will count down
755 // on this field. We also cache the row lengths because we need this
756 // data frequently and getting it from the sparsity pattern is more
757 // expensive.
758 for (types::global_dof_index row = 0; row < n_nodes; ++row)
759 {
760 row_lengths[row] = connectivity.row_length(row);
761 Assert(row_lengths[row] > 0, ExcInternalError());
762 }
763 std::vector<unsigned int> n_remaining_neighbors(row_lengths);
764
765 // This outer loop is typically traversed only once, unless the global
766 // graph is not connected
767 while (true)
768 {
769 // Find node with the minimal number of neighbors (typically a
770 // corner node when based on FEM meshes). If no node is left, we are
771 // done. Together with the outer while loop, this loop can possibly
772 // be of quadratic complexity in the number of disconnected
773 // partitions, i.e. up to n_nodes in the worst case,
774 // but that is not the usual use case of this loop and thus not
775 // optimized for.
776 {
777 unsigned int candidate_valence = numbers::invalid_unsigned_int;
778 types::global_dof_index candidate_index =
780 for (types::global_dof_index i = 0; i < n_nodes; ++i)
781 if (touched_nodes[i] == unseen_node)
782 if (row_lengths[i] < candidate_valence)
783 {
784 candidate_index = i;
785 candidate_valence = n_remaining_neighbors[i];
786 if (candidate_valence <= 1)
787 break;
788 }
789 if (candidate_index == numbers::invalid_dof_index)
790 break;
791
792 Assert(candidate_valence > 0, ExcInternalError());
793
794 current_neighbors = {candidate_index};
795 touched_nodes[candidate_index] = available_node;
796 }
797
798 while (true)
799 {
800 // Find node with minimum number of untouched neighbors among
801 // the next set of possible neighbors (= valence), and among the
802 // set of nodes with the minimal number of neighbors, choose the
803 // one with the largest number of touched neighbors (i.e., the
804 // largest row length).
805 //
806 // This loop is typically the most expensive part for large
807 // graphs and thus only run once. We also do some cleanup, i.e.,
808 // the indices added to a group in the previous round need to be
809 // removed at this point.
810 unsigned int candidate_valence = numbers::invalid_unsigned_int;
811 types::global_dof_index candidate_index =
813 unsigned int candidate_row_length = 0;
814 const unsigned int loop_length = current_neighbors.size();
815 unsigned int write_index = 0;
816 for (unsigned int i = 0; i < loop_length; ++i)
817 {
818 const types::global_dof_index node = current_neighbors[i];
819 Assert(touched_nodes[node] != unseen_node,
821 if (touched_nodes[node] == available_node)
822 {
823 current_neighbors[write_index] = node;
824 ++write_index;
825 if (n_remaining_neighbors[node] < candidate_valence ||
826 (n_remaining_neighbors[node] == candidate_valence &&
827 (row_lengths[node] > candidate_row_length ||
828 (row_lengths[node] == candidate_row_length &&
829 node < candidate_index))))
830 {
831 candidate_index = node;
832 candidate_valence = n_remaining_neighbors[node];
833 candidate_row_length = row_lengths[node];
834 }
835 }
836 }
837 current_neighbors.resize(write_index);
838
839#ifdef DEBUG
840 for (const types::global_dof_index node : current_neighbors)
841 Assert(touched_nodes[node] == available_node,
843#endif
844
845 // No more neighbors left -> terminate loop
846 if (current_neighbors.empty())
847 break;
848
849 // Add the pivot and all direct neighbors of the pivot node not
850 // yet touched to the list of new entries.
851 group_indices.push_back(candidate_index);
852 touched_nodes[candidate_index] = group_starts.size() - 1;
853 const auto end_it = connectivity.end(candidate_index);
854 for (auto it = connectivity.begin(candidate_index); it != end_it;
855 ++it)
856 if (touched_nodes[it->column()] >= available_node)
857 {
858 group_indices.push_back(it->column());
859 touched_nodes[it->column()] = group_starts.size() - 1;
860 }
861 group_starts.push_back(group_indices.size());
862
863 // Add all neighbors of the current list not yet seen to the set
864 // of possible next nodes. The added node is grouped and thus no
865 // longer a valid neighbor (here we assume symmetry of the
866 // connectivity). It will be removed from the list of neighbors
867 // by the code further up in the next iteration of the
868 // surrounding loop.
869 for (types::global_dof_index index =
870 group_starts[group_starts.size() - 2];
871 index < group_starts.back();
872 ++index)
873 {
874 auto it = connectivity.begin(group_indices[index]);
875 const auto end_row = connectivity.end(group_indices[index]);
876 for (; it != end_row; ++it)
877 {
878 if (touched_nodes[it->column()] == unseen_node)
879 {
880 current_neighbors.push_back(it->column());
881 touched_nodes[it->column()] = available_node;
882 }
883 n_remaining_neighbors[it->column()]--;
884 }
885 }
886 }
887 }
888
889 // Sanity check: for all nodes, there should not be any neighbors left
890 for (types::global_dof_index row = 0; row < n_nodes; ++row)
891 Assert(n_remaining_neighbors[row] == 0, ExcInternalError());
892
893 // If the number of groups is smaller than the number of nodes, we
894 // continue by recursively calling this method
895 const unsigned int n_groups = group_starts.size() - 1;
896 if (n_groups < n_nodes)
897 {
898 // Form the connectivity of the groups
899 DynamicSparsityPattern connectivity_next(n_groups, n_groups);
900 for (types::global_dof_index row = 0; row < n_groups; ++row)
901 for (types::global_dof_index index = group_starts[row];
902 index < group_starts[row + 1];
903 ++index)
904 {
905 auto it = connectivity.begin(group_indices[index]);
906 const auto end_it = connectivity.end(group_indices[index]);
907 for (; it != end_it; ++it)
908 connectivity_next.add(row, touched_nodes[it->column()]);
909 }
910
911 // Recursively call the reordering
912 std::vector<types::global_dof_index> renumbering_next(n_groups);
913 reorder_hierarchical(connectivity_next, renumbering_next);
914
915 // Renumber the indices group by group according to the incoming
916 // ordering for the groups
917 for (types::global_dof_index row = 0, c = 0; row < n_groups; ++row)
918 for (types::global_dof_index index =
919 group_starts[renumbering_next[row]];
920 index < group_starts[renumbering_next[row] + 1];
921 ++index, ++c)
922 renumbering[c] = group_indices[index];
923 }
924 else
925 {
926 // All groups should have size one and no more recursion is possible,
927 // so use the numbering of the groups
928 unsigned int c = 0;
929 for (const types::global_dof_index i : group_indices)
930 renumbering[c++] = i;
931 }
932 }
933 } // namespace internal
934
935 void
937 const DynamicSparsityPattern &connectivity,
938 std::vector<DynamicSparsityPattern::size_type> &renumbering)
939 {
940 // the internal renumbering keeps the numbering the wrong way around (but
941 // we cannot invert the numbering inside that method because it is used
942 // recursively), so invert it here
943 internal::reorder_hierarchical(connectivity, renumbering);
944 renumbering = Utilities::invert_permutation(renumbering);
945 }
946
947
948
949#ifdef DEAL_II_WITH_MPI
950
951 void
953 const IndexSet &locally_owned_rows,
954 const MPI_Comm mpi_comm,
955 const IndexSet &locally_relevant_rows)
956 {
957 using map_vec_t =
958 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
959
960 // 1. limit rows to non owned:
961 IndexSet requested_rows(locally_relevant_rows);
962 requested_rows.subtract_set(locally_owned_rows);
963
964 std::vector<unsigned int> index_owner =
965 Utilities::MPI::compute_index_owner(locally_owned_rows,
966 requested_rows,
967 mpi_comm);
968
969 // 2. go through requested_rows, figure out the owner and add the row to
970 // request
971 map_vec_t rows_data;
973 i < requested_rows.n_elements();
974 ++i)
975 {
977 requested_rows.nth_index_in_set(i);
978
979 rows_data[index_owner[i]].push_back(row);
980 }
981
982 // 3. get what others ask us to send
983 const auto rows_data_received =
984 Utilities::MPI::some_to_some(mpi_comm, rows_data);
985
986 // 4. now prepare data to be sent in the same format as in
987 // distribute_sparsity_pattern() below, i.e.
988 // rX,num_rX,cols_rX
989 map_vec_t send_data;
990 for (const auto &data : rows_data_received)
991 {
992 for (const auto &row : data.second)
993 {
994 const auto rlen = dsp.row_length(row);
995
996 // skip empty lines
997 if (rlen == 0)
998 continue;
999
1000 // save entries
1001 send_data[data.first].push_back(row); // row index
1002 send_data[data.first].push_back(rlen); // number of entries
1003 for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1004 send_data[data.first].push_back(
1005 dsp.column_number(row, c)); // columns
1006 } // loop over rows
1007 } // loop over received data
1008
1009 // 5. communicate rows
1010 const auto received_data =
1011 Utilities::MPI::some_to_some(mpi_comm, send_data);
1012
1013 // 6. add result to our sparsity
1014 for (const auto &data : received_data)
1015 {
1016 const auto &recv_buf = data.second;
1017 auto ptr = recv_buf.begin();
1018 const auto end = recv_buf.end();
1019 while (ptr != end)
1020 {
1021 const auto row = *(ptr++);
1022 Assert(ptr != end, ExcInternalError());
1023
1024 const auto n_entries = *(ptr++);
1025 Assert(n_entries > 0, ExcInternalError());
1026 Assert(ptr != end, ExcInternalError());
1027
1028 // make sure we clear whatever was previously stored
1029 // in these rows. Otherwise we can't guarantee that the
1030 // data is consistent across MPI communicator.
1031 dsp.clear_row(row);
1032
1033 Assert(ptr + (n_entries - 1) != end, ExcInternalError());
1034 dsp.add_entries(row, ptr, ptr + n_entries, true);
1035 ptr += n_entries;
1036 }
1037 Assert(ptr == end, ExcInternalError());
1038 }
1039 }
1040
1041
1042
1043 void
1046 const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu,
1047 const MPI_Comm mpi_comm,
1048 const IndexSet &myrange)
1049 {
1050 const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1051 std::vector<DynamicSparsityPattern::size_type> start_index(
1052 rows_per_cpu.size() + 1);
1053 start_index[0] = 0;
1054 for (DynamicSparsityPattern::size_type i = 0; i < rows_per_cpu.size(); ++i)
1055 start_index[i + 1] = start_index[i] + rows_per_cpu[i];
1056
1057 IndexSet owned(start_index.back());
1058 owned.add_range(start_index[myid], start_index[myid] + rows_per_cpu[myid]);
1059
1060 distribute_sparsity_pattern(dsp, owned, mpi_comm, myrange);
1061 }
1062
1063
1064
1065 void
1067 const IndexSet &locally_owned_rows,
1068 const MPI_Comm mpi_comm,
1069 const IndexSet &locally_relevant_rows)
1070 {
1072 dsp.row_index_set() == locally_relevant_rows,
1073 ExcMessage(
1074 "The DynamicSparsityPattern must be initialized with an IndexSet that contains locally relevant indices."));
1075
1076 IndexSet requested_rows(locally_relevant_rows);
1077 requested_rows.subtract_set(locally_owned_rows);
1078
1079 std::vector<unsigned int> index_owner =
1080 Utilities::MPI::compute_index_owner(locally_owned_rows,
1081 requested_rows,
1082 mpi_comm);
1083
1084 using map_vec_t =
1085 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
1086
1087 map_vec_t send_data;
1088
1090 i < requested_rows.n_elements();
1091 ++i)
1092 {
1094 requested_rows.nth_index_in_set(i);
1095
1096 const auto rlen = dsp.row_length(row);
1097
1098 // skip empty lines
1099 if (rlen == 0)
1100 continue;
1101
1102 // save entries
1103 send_data[index_owner[i]].push_back(row); // row index
1104 send_data[index_owner[i]].push_back(rlen); // number of entries
1105 for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1106 {
1107 // columns
1108 const auto column = dsp.column_number(row, c);
1109 send_data[index_owner[i]].push_back(column);
1110 }
1111 }
1112
1113 const auto receive_data = Utilities::MPI::some_to_some(mpi_comm, send_data);
1114
1115 // add what we received
1116 for (const auto &data : receive_data)
1117 {
1118 const auto &recv_buf = data.second;
1119 auto ptr = recv_buf.begin();
1120 const auto end = recv_buf.end();
1121 while (ptr != end)
1122 {
1123 const auto row = *(ptr++);
1124 Assert(ptr != end, ExcInternalError());
1125 const auto n_entries = *(ptr++);
1126
1127 Assert(ptr + (n_entries - 1) != end, ExcInternalError());
1128 dsp.add_entries(row, ptr, ptr + n_entries, true);
1129 ptr += n_entries;
1130 }
1131 Assert(ptr == end, ExcInternalError());
1132 }
1133 }
1134
1135
1136
1137 void
1139 const std::vector<IndexSet> &owned_set_per_cpu,
1140 const MPI_Comm mpi_comm,
1141 const IndexSet &myrange)
1142 {
1143 const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1145 owned_set_per_cpu[myid],
1146 mpi_comm,
1147 myrange);
1148 }
1149
1150
1151
1152 void
1154 const IndexSet &locally_owned_rows,
1155 const MPI_Comm mpi_comm,
1156 const IndexSet &locally_relevant_rows)
1157 {
1158 using map_vec_t =
1160 std::vector<BlockDynamicSparsityPattern::size_type>>;
1161 map_vec_t send_data;
1162
1163 IndexSet requested_rows(locally_relevant_rows);
1164 requested_rows.subtract_set(locally_owned_rows);
1165
1166 std::vector<unsigned int> index_owner =
1167 Utilities::MPI::compute_index_owner(locally_owned_rows,
1168 requested_rows,
1169 mpi_comm);
1170
1172 i < requested_rows.n_elements();
1173 ++i)
1174 {
1176 requested_rows.nth_index_in_set(i);
1177
1179
1180 // skip empty lines
1181 if (rlen == 0)
1182 continue;
1183
1184 // save entries
1185 std::vector<BlockDynamicSparsityPattern::size_type> &dst =
1186 send_data[index_owner[i]];
1187
1188 dst.push_back(rlen); // number of entries
1189 dst.push_back(row); // row index
1190 for (BlockDynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1191 {
1192 // columns
1194 dsp.column_number(row, c);
1195 dst.push_back(column);
1196 }
1197 }
1198
1199 unsigned int num_receive = 0;
1200 {
1201 std::vector<unsigned int> send_to;
1202 send_to.reserve(send_data.size());
1203 for (const auto &sparsity_line : send_data)
1204 send_to.push_back(sparsity_line.first);
1205
1206 num_receive =
1208 send_to);
1209 }
1210
1211 std::vector<MPI_Request> requests(send_data.size());
1212
1213
1214 // send data
1215
1217 Utilities::MPI::CollectiveMutex::ScopedLock lock(mutex, mpi_comm);
1218
1219 const int mpi_tag = Utilities::MPI::internal::Tags::
1221
1222 {
1223 unsigned int idx = 0;
1224 for (const auto &sparsity_line : send_data)
1225 {
1226 const int ierr = MPI_Isend(sparsity_line.second.data(),
1227 sparsity_line.second.size(),
1229 sparsity_line.first,
1230 mpi_tag,
1231 mpi_comm,
1232 &requests[idx++]);
1233 AssertThrowMPI(ierr);
1234 }
1235 }
1236
1237 {
1238 // receive
1239 std::vector<BlockDynamicSparsityPattern::size_type> recv_buf;
1240 for (unsigned int index = 0; index < num_receive; ++index)
1241 {
1242 MPI_Status status;
1243 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, mpi_comm, &status);
1244 AssertThrowMPI(ierr);
1245
1246 int len;
1247 ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
1248 AssertThrowMPI(ierr);
1249
1250 recv_buf.resize(len);
1251 ierr = MPI_Recv(recv_buf.data(),
1252 len,
1254 status.MPI_SOURCE,
1255 status.MPI_TAG,
1256 mpi_comm,
1257 &status);
1258 AssertThrowMPI(ierr);
1259
1260 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1261 ptr = recv_buf.begin();
1262 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1263 end = recv_buf.end();
1264 while (ptr != end)
1265 {
1267 Assert(ptr != end, ExcInternalError());
1269 for (unsigned int c = 0; c < num; ++c)
1270 {
1271 Assert(ptr != end, ExcInternalError());
1272 dsp.add(row, *ptr);
1273 ++ptr;
1274 }
1275 }
1276 Assert(ptr == end, ExcInternalError());
1277 }
1278 }
1279
1280 // complete all sends, so that we can safely destroy the buffers.
1281 if (requests.size() > 0)
1282 {
1283 const int ierr =
1284 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1285 AssertThrowMPI(ierr);
1286 }
1287 }
1288#endif
1289} // namespace SparsityTools
1290
size_type column_number(const size_type row, const unsigned int index) const
unsigned int row_length(const size_type row) const
void add(const size_type i, const size_type j)
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
const IndexSet & row_index_set() const
size_type row_length(const size_type row) const
size_type column_number(const size_type row, const size_type index) const
void clear_row(const size_type row)
void add(const size_type i, const size_type j)
size_type size() const
Definition index_set.h:1776
size_type n_elements() const
Definition index_set.h:1934
void subtract_set(const IndexSet &other)
Definition index_set.cc:473
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1803
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1982
types::global_dof_index size_type
size_type n_rows() const
size_type n_cols() const
bool is_compressed() const
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
iterator begin() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
iterator end() const
unsigned int row_length(const size_type row) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcMETISNotInstalled()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMETISError(int arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
static ::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcZOLTANNotInstalled()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
DynamicSparsityPattern::size_type find_unnumbered_starting_index(const DynamicSparsityPattern &sparsity, const std::vector< DynamicSparsityPattern::size_type > &new_indices)
void reorder_hierarchical(const DynamicSparsityPattern &connectivity, std::vector< DynamicSparsityPattern::size_type > &renumbering)
unsigned int color_sparsity_pattern(const SparsityPattern &sparsity_pattern, std::vector< unsigned int > &color_indices)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
void reorder_hierarchical(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
void gather_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
void reorder_Cuthill_McKee(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices, const std::vector< DynamicSparsityPattern::size_type > &starting_indices=std::vector< DynamicSparsityPattern::size_type >())
@ sparsity_tools_distribute_sparsity_pattern
SparsityTools::sparsity_tools_distribute_sparsity_pattern()
Definition mpi_tags.h:73
std::map< unsigned int, T > some_to_some(const MPI_Comm comm, const std::map< unsigned int, T > &objects_to_send)
std::vector< unsigned int > compute_index_owner(const IndexSet &owned_indices, const IndexSet &indices_to_look_up, const MPI_Comm comm)
Definition mpi.cc:697
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
unsigned int compute_n_point_to_point_communications(const MPI_Comm mpi_comm, const std::vector< unsigned int > &destinations)
Definition mpi.cc:362
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1699
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
const types::global_dof_index invalid_size_type
Definition types.h:233
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition types.h:102