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