Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+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 - 2023 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 std::vector<types::global_dof_index> touched_nodes(
730 connectivity.n_rows(), numbers::invalid_dof_index);
731 std::vector<unsigned int> row_lengths(connectivity.n_rows());
732 std::set<types::global_dof_index> current_neighbors;
733 std::vector<std::vector<types::global_dof_index>> groups;
734
735 // First collect the number of neighbors for each node. We use this
736 // field to find next nodes with the minimum number of non-touched
737 // neighbors in the field n_remaining_neighbors, so we will count down
738 // on this field. We also cache the row lengths because we need this
739 // data frequently and getting it from the sparsity pattern is more
740 // expensive.
741 for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row)
742 {
743 row_lengths[row] = connectivity.row_length(row);
744 Assert(row_lengths[row] > 0, ExcInternalError());
745 }
746 std::vector<unsigned int> n_remaining_neighbors(row_lengths);
747
748 // This outer loop is typically traversed only once, unless the global
749 // graph is not connected
750 while (true)
751 {
752 // Find cell with the minimal number of neighbors (typically a
753 // corner node when based on FEM meshes). If no cell is left, we are
754 // done. Together with the outer while loop, this loop can possibly
755 // be of quadratic complexity in the number of disconnected
756 // partitions, i.e. up to connectivity.n_rows() in the worst case,
757 // but that is not the usual use case of this loop and thus not
758 // optimized for.
759 std::pair<types::global_dof_index, types::global_dof_index>
760 min_neighbors(numbers::invalid_dof_index,
762 for (types::global_dof_index i = 0; i < touched_nodes.size(); ++i)
763 if (touched_nodes[i] == numbers::invalid_dof_index)
764 if (row_lengths[i] < min_neighbors.second)
765 {
766 min_neighbors = std::make_pair(i, n_remaining_neighbors[i]);
767 if (n_remaining_neighbors[i] <= 1)
768 break;
769 }
770 if (min_neighbors.first == numbers::invalid_dof_index)
771 break;
772
773 Assert(min_neighbors.second > 0, ExcInternalError());
774
775 current_neighbors.clear();
776 current_neighbors.insert(min_neighbors.first);
777 while (!current_neighbors.empty())
778 {
779 // Find node with minimum number of untouched neighbors among the
780 // next set of possible neighbors
781 min_neighbors = std::make_pair(numbers::invalid_dof_index,
783 for (const auto current_neighbor : current_neighbors)
784 {
785 Assert(touched_nodes[current_neighbor] ==
788 if (n_remaining_neighbors[current_neighbor] <
789 min_neighbors.second)
790 min_neighbors =
791 std::make_pair(current_neighbor,
792 n_remaining_neighbors[current_neighbor]);
793 }
794
795 // Among the set of nodes with the minimal number of neighbors,
796 // choose the one with the largest number of touched neighbors,
797 // i.e., the one with the largest row length
798 const types::global_dof_index best_row_length =
799 min_neighbors.second;
800 for (const auto current_neighbor : current_neighbors)
801 if (n_remaining_neighbors[current_neighbor] == best_row_length)
802 if (row_lengths[current_neighbor] > min_neighbors.second)
803 min_neighbors =
804 std::make_pair(current_neighbor,
805 row_lengths[current_neighbor]);
806
807 // Add the pivot and all direct neighbors of the pivot node not
808 // yet touched to the list of new entries.
809 groups.emplace_back();
810 std::vector<types::global_dof_index> &next_group = groups.back();
811
812 next_group.push_back(min_neighbors.first);
813 touched_nodes[min_neighbors.first] = groups.size() - 1;
815 connectivity.begin(min_neighbors.first);
816 it != connectivity.end(min_neighbors.first);
817 ++it)
818 if (touched_nodes[it->column()] == numbers::invalid_dof_index)
819 {
820 next_group.push_back(it->column());
821 touched_nodes[it->column()] = groups.size() - 1;
822 }
823
824 // Add all neighbors of the current list not yet touched to the
825 // set of possible next pivots. The added node is no longer a
826 // valid neighbor (here we assume symmetry of the
827 // connectivity). Delete the entries of the current list from
828 // the set of possible next pivots.
829 for (const auto index : next_group)
830 {
832 connectivity.begin(index);
833 it != connectivity.end(index);
834 ++it)
835 {
836 if (touched_nodes[it->column()] ==
838 current_neighbors.insert(it->column());
839 n_remaining_neighbors[it->column()]--;
840 }
841 current_neighbors.erase(index);
842 }
843 }
844 }
845
846 // Sanity check: for all nodes, there should not be any neighbors left
847 for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row)
848 Assert(n_remaining_neighbors[row] == 0, ExcInternalError());
849
850 // If the number of groups is smaller than the number of nodes, we
851 // continue by recursively calling this method
852 if (groups.size() < connectivity.n_rows())
853 {
854 // Form the connectivity of the groups
855 DynamicSparsityPattern connectivity_next(groups.size(),
856 groups.size());
857 for (types::global_dof_index i = 0; i < groups.size(); ++i)
858 for (types::global_dof_index col = 0; col < groups[i].size(); ++col)
860 connectivity.begin(groups[i][col]);
861 it != connectivity.end(groups[i][col]);
862 ++it)
863 connectivity_next.add(i, touched_nodes[it->column()]);
864
865 // Recursively call the reordering
866 std::vector<types::global_dof_index> renumbering_next(groups.size());
867 reorder_hierarchical(connectivity_next, renumbering_next);
868
869 // Renumber the indices group by group according to the incoming
870 // ordering for the groups
871 for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i)
872 for (types::global_dof_index col = 0;
873 col < groups[renumbering_next[i]].size();
874 ++col, ++count)
875 renumbering[count] = groups[renumbering_next[i]][col];
876 }
877 else
878 {
879 // All groups should have size one and no more recursion is possible,
880 // so use the numbering of the groups
881 for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i)
882 for (types::global_dof_index col = 0; col < groups[i].size();
883 ++col, ++count)
884 renumbering[count] = groups[i][col];
885 }
886 }
887 } // namespace internal
888
889 void
891 const DynamicSparsityPattern &connectivity,
892 std::vector<DynamicSparsityPattern::size_type> &renumbering)
893 {
894 // the internal renumbering keeps the numbering the wrong way around (but
895 // we cannot invert the numbering inside that method because it is used
896 // recursively), so invert it here
897 internal::reorder_hierarchical(connectivity, renumbering);
898 renumbering = Utilities::invert_permutation(renumbering);
899 }
900
901
902
903#ifdef DEAL_II_WITH_MPI
904
905 void
907 const IndexSet &locally_owned_rows,
908 const MPI_Comm mpi_comm,
909 const IndexSet &locally_relevant_rows)
910 {
911 using map_vec_t =
912 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
913
914 // 1. limit rows to non owned:
915 IndexSet requested_rows(locally_relevant_rows);
916 requested_rows.subtract_set(locally_owned_rows);
917
918 std::vector<unsigned int> index_owner =
919 Utilities::MPI::compute_index_owner(locally_owned_rows,
920 requested_rows,
921 mpi_comm);
922
923 // 2. go through requested_rows, figure out the owner and add the row to
924 // request
925 map_vec_t rows_data;
927 i < requested_rows.n_elements();
928 ++i)
929 {
931 requested_rows.nth_index_in_set(i);
932
933 rows_data[index_owner[i]].push_back(row);
934 }
935
936 // 3. get what others ask us to send
937 const auto rows_data_received =
938 Utilities::MPI::some_to_some(mpi_comm, rows_data);
939
940 // 4. now prepare data to be sent in the same format as in
941 // distribute_sparsity_pattern() below, i.e.
942 // rX,num_rX,cols_rX
943 map_vec_t send_data;
944 for (const auto &data : rows_data_received)
945 {
946 for (const auto &row : data.second)
947 {
948 const auto rlen = dsp.row_length(row);
949
950 // skip empty lines
951 if (rlen == 0)
952 continue;
953
954 // save entries
955 send_data[data.first].push_back(row); // row index
956 send_data[data.first].push_back(rlen); // number of entries
957 for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
958 send_data[data.first].push_back(
959 dsp.column_number(row, c)); // columns
960 } // loop over rows
961 } // loop over received data
962
963 // 5. communicate rows
964 const auto received_data =
965 Utilities::MPI::some_to_some(mpi_comm, send_data);
966
967 // 6. add result to our sparsity
968 for (const auto &data : received_data)
969 {
970 const auto &recv_buf = data.second;
971 auto ptr = recv_buf.begin();
972 const auto end = recv_buf.end();
973 while (ptr != end)
974 {
975 const auto row = *(ptr++);
976 Assert(ptr != end, ExcInternalError());
977
978 const auto n_entries = *(ptr++);
979 Assert(n_entries > 0, ExcInternalError());
980 Assert(ptr != end, ExcInternalError());
981
982 // make sure we clear whatever was previously stored
983 // in these rows. Otherwise we can't guarantee that the
984 // data is consistent across MPI communicator.
985 dsp.clear_row(row);
986
987 Assert(ptr + (n_entries - 1) != end, ExcInternalError());
988 dsp.add_entries(row, ptr, ptr + n_entries, true);
989 ptr += n_entries;
990 }
991 Assert(ptr == end, ExcInternalError());
992 }
993 }
994
995
996
997 void
1000 const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu,
1001 const MPI_Comm mpi_comm,
1002 const IndexSet &myrange)
1003 {
1004 const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1005 std::vector<DynamicSparsityPattern::size_type> start_index(
1006 rows_per_cpu.size() + 1);
1007 start_index[0] = 0;
1008 for (DynamicSparsityPattern::size_type i = 0; i < rows_per_cpu.size(); ++i)
1009 start_index[i + 1] = start_index[i] + rows_per_cpu[i];
1010
1011 IndexSet owned(start_index.back());
1012 owned.add_range(start_index[myid], start_index[myid] + rows_per_cpu[myid]);
1013
1014 distribute_sparsity_pattern(dsp, owned, mpi_comm, myrange);
1015 }
1016
1017
1018
1019 void
1021 const IndexSet &locally_owned_rows,
1022 const MPI_Comm mpi_comm,
1023 const IndexSet &locally_relevant_rows)
1024 {
1025 IndexSet requested_rows(locally_relevant_rows);
1026 requested_rows.subtract_set(locally_owned_rows);
1027
1028 std::vector<unsigned int> index_owner =
1029 Utilities::MPI::compute_index_owner(locally_owned_rows,
1030 requested_rows,
1031 mpi_comm);
1032
1033 using map_vec_t =
1034 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
1035
1036 map_vec_t send_data;
1037
1039 i < requested_rows.n_elements();
1040 ++i)
1041 {
1043 requested_rows.nth_index_in_set(i);
1044
1045 const auto rlen = dsp.row_length(row);
1046
1047 // skip empty lines
1048 if (rlen == 0)
1049 continue;
1050
1051 // save entries
1052 send_data[index_owner[i]].push_back(row); // row index
1053 send_data[index_owner[i]].push_back(rlen); // number of entries
1054 for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1055 {
1056 // columns
1057 const auto column = dsp.column_number(row, c);
1058 send_data[index_owner[i]].push_back(column);
1059 }
1060 }
1061
1062 const auto receive_data = Utilities::MPI::some_to_some(mpi_comm, send_data);
1063
1064 // add what we received
1065 for (const auto &data : receive_data)
1066 {
1067 const auto &recv_buf = data.second;
1068 auto ptr = recv_buf.begin();
1069 const auto end = recv_buf.end();
1070 while (ptr != end)
1071 {
1072 const auto row = *(ptr++);
1073 Assert(ptr != end, ExcInternalError());
1074 const auto n_entries = *(ptr++);
1075
1076 Assert(ptr + (n_entries - 1) != end, ExcInternalError());
1077 dsp.add_entries(row, ptr, ptr + n_entries, true);
1078 ptr += n_entries;
1079 }
1080 Assert(ptr == end, ExcInternalError());
1081 }
1082 }
1083
1084
1085
1086 void
1088 const std::vector<IndexSet> &owned_set_per_cpu,
1089 const MPI_Comm mpi_comm,
1090 const IndexSet &myrange)
1091 {
1092 const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1094 owned_set_per_cpu[myid],
1095 mpi_comm,
1096 myrange);
1097 }
1098
1099
1100
1101 void
1103 const IndexSet &locally_owned_rows,
1104 const MPI_Comm mpi_comm,
1105 const IndexSet &locally_relevant_rows)
1106 {
1107 using map_vec_t =
1109 std::vector<BlockDynamicSparsityPattern::size_type>>;
1110 map_vec_t send_data;
1111
1112 IndexSet requested_rows(locally_relevant_rows);
1113 requested_rows.subtract_set(locally_owned_rows);
1114
1115 std::vector<unsigned int> index_owner =
1116 Utilities::MPI::compute_index_owner(locally_owned_rows,
1117 requested_rows,
1118 mpi_comm);
1119
1121 i < requested_rows.n_elements();
1122 ++i)
1123 {
1125 requested_rows.nth_index_in_set(i);
1126
1128
1129 // skip empty lines
1130 if (rlen == 0)
1131 continue;
1132
1133 // save entries
1134 std::vector<BlockDynamicSparsityPattern::size_type> &dst =
1135 send_data[index_owner[i]];
1136
1137 dst.push_back(rlen); // number of entries
1138 dst.push_back(row); // row index
1139 for (BlockDynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1140 {
1141 // columns
1143 dsp.column_number(row, c);
1144 dst.push_back(column);
1145 }
1146 }
1147
1148 unsigned int num_receive = 0;
1149 {
1150 std::vector<unsigned int> send_to;
1151 send_to.reserve(send_data.size());
1152 for (const auto &sparsity_line : send_data)
1153 send_to.push_back(sparsity_line.first);
1154
1155 num_receive =
1157 send_to);
1158 }
1159
1160 std::vector<MPI_Request> requests(send_data.size());
1161
1162
1163 // send data
1164
1166 Utilities::MPI::CollectiveMutex::ScopedLock lock(mutex, mpi_comm);
1167
1168 const int mpi_tag = Utilities::MPI::internal::Tags::
1170
1171 {
1172 unsigned int idx = 0;
1173 for (const auto &sparsity_line : send_data)
1174 {
1175 const int ierr = MPI_Isend(sparsity_line.second.data(),
1176 sparsity_line.second.size(),
1178 sparsity_line.first,
1179 mpi_tag,
1180 mpi_comm,
1181 &requests[idx++]);
1182 AssertThrowMPI(ierr);
1183 }
1184 }
1185
1186 {
1187 // receive
1188 std::vector<BlockDynamicSparsityPattern::size_type> recv_buf;
1189 for (unsigned int index = 0; index < num_receive; ++index)
1190 {
1191 MPI_Status status;
1192 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, mpi_comm, &status);
1193 AssertThrowMPI(ierr);
1194
1195 int len;
1196 ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
1197 AssertThrowMPI(ierr);
1198
1199 recv_buf.resize(len);
1200 ierr = MPI_Recv(recv_buf.data(),
1201 len,
1203 status.MPI_SOURCE,
1204 status.MPI_TAG,
1205 mpi_comm,
1206 &status);
1207 AssertThrowMPI(ierr);
1208
1209 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1210 ptr = recv_buf.begin();
1211 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1212 end = recv_buf.end();
1213 while (ptr != end)
1214 {
1216 Assert(ptr != end, ExcInternalError());
1218 for (unsigned int c = 0; c < num; ++c)
1219 {
1220 Assert(ptr != end, ExcInternalError());
1221 dsp.add(row, *ptr);
1222 ++ptr;
1223 }
1224 }
1225 Assert(ptr == end, ExcInternalError());
1226 }
1227 }
1228
1229 // complete all sends, so that we can safely destroy the buffers.
1230 if (requests.size() > 0)
1231 {
1232 const int ierr =
1233 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1234 AssertThrowMPI(ierr);
1235 }
1236 }
1237#endif
1238} // namespace SparsityTools
1239
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:1765
size_type n_elements() const
Definition index_set.h:1923
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1792
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1971
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:74
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:1068
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:143
unsigned int compute_n_point_to_point_communications(const MPI_Comm mpi_comm, const std::vector< unsigned int > &destinations)
Definition mpi.cc:415
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1696
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