deal.II version GIT relicensing-3509-g790ffced1c 2025-06-14 07: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
mpi.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 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
18#include <deal.II/base/mpi.h>
19#include <deal.II/base/mpi.templates.h>
24
27
28#include <boost/serialization/utility.hpp>
29
30// In this file, we use offsetof, which is a macro. When compiling
31// with C++20 modules, this presents a problem because we wrap all of
32// namespace std -- and then don't have access to macros. As a
33// consequence, we really do need the following #include, even when
34// building modules:
35#include <cstddef> // Do not convert for module purposes
36#include <iostream>
37#include <limits>
38#include <numeric>
39#include <set>
40#include <vector>
41
42#if defined(DEAL_II_WITH_MPI)
44# include <mpi.h>
46#endif
47
48
50
51
52namespace Utilities
53{
56 const unsigned int my_partition_id,
57 const unsigned int n_partitions,
58 const types::global_dof_index total_size)
59 {
60 static_assert(std::is_same_v<types::global_dof_index, IndexSet::size_type>,
61 "IndexSet::size_type must match types::global_dof_index for "
62 "using this function");
63 const unsigned int remain = total_size % n_partitions;
64
65 const IndexSet::size_type min_size = total_size / n_partitions;
66
67 const IndexSet::size_type begin =
68 min_size * my_partition_id + std::min(my_partition_id, remain);
69 const IndexSet::size_type end =
70 min_size * (my_partition_id + 1) + std::min(my_partition_id + 1, remain);
71 IndexSet result(total_size);
72 result.add_range(begin, end);
73 return result;
74 }
75
76 namespace MPI
77 {
78 MinMaxAvg
79 min_max_avg(const double my_value, const MPI_Comm mpi_communicator)
80 {
81 MinMaxAvg result;
84 mpi_communicator);
85
86 return result;
87 }
88
89
90
91 std::vector<MinMaxAvg>
92 min_max_avg(const std::vector<double> &my_values,
93 const MPI_Comm mpi_communicator)
94 {
95 std::vector<MinMaxAvg> results(my_values.size());
96 min_max_avg(my_values, results, mpi_communicator);
97
98 return results;
99 }
100
101
102
103#ifdef DEAL_II_WITH_MPI
104 unsigned int
105 n_mpi_processes(const MPI_Comm mpi_communicator)
106 {
107 if (job_supports_mpi())
108 {
109 int n_jobs = 1;
110 const int ierr = MPI_Comm_size(mpi_communicator, &n_jobs);
111 AssertThrowMPI(ierr);
112 return n_jobs;
113 }
114 else
115 return 1;
116 }
117
118
119 unsigned int
120 this_mpi_process(const MPI_Comm mpi_communicator)
121 {
122 if (job_supports_mpi())
123 {
124 int rank = 0;
125 const int ierr = MPI_Comm_rank(mpi_communicator, &rank);
126 AssertThrowMPI(ierr);
127 return rank;
128 }
129 else
130 return 0;
131 }
132
133
134
135 std::vector<unsigned int>
137 const MPI_Comm comm_small)
138 {
140 return std::vector<unsigned int>{0};
141
142 const unsigned int rank = Utilities::MPI::this_mpi_process(comm_large);
143 const unsigned int size = Utilities::MPI::n_mpi_processes(comm_small);
144
145 std::vector<unsigned int> ranks(size);
146 const int ierr = MPI_Allgather(
147 &rank, 1, MPI_UNSIGNED, ranks.data(), 1, MPI_UNSIGNED, comm_small);
148 AssertThrowMPI(ierr);
149
150 return ranks;
151 }
152
153
154
156 duplicate_communicator(const MPI_Comm mpi_communicator)
157 {
158 MPI_Comm new_communicator;
159 const int ierr = MPI_Comm_dup(mpi_communicator, &new_communicator);
160 AssertThrowMPI(ierr);
161 return new_communicator;
162 }
163
164
165
166 void
167 free_communicator(MPI_Comm mpi_communicator)
168 {
169 // MPI_Comm_free will set the argument to MPI_COMM_NULL automatically.
170 const int ierr = MPI_Comm_free(&mpi_communicator);
171 AssertThrowMPI(ierr);
172 }
173
174
175
176 std::vector<IndexSet>
178 const MPI_Comm comm,
180 {
181 static_assert(
182 std::is_same_v<types::global_dof_index, IndexSet::size_type>,
183 "IndexSet::size_type must match types::global_dof_index for "
184 "using this function");
185 const unsigned int n_proc = n_mpi_processes(comm);
186 const std::vector<IndexSet::size_type> sizes =
188 const auto total_size =
189 std::accumulate(sizes.begin(), sizes.end(), IndexSet::size_type(0));
190
191 std::vector<IndexSet> res(n_proc, IndexSet(total_size));
192
193 IndexSet::size_type begin = 0;
194 for (unsigned int i = 0; i < n_proc; ++i)
195 {
196 res[i].add_range(begin, begin + sizes[i]);
197 begin = begin + sizes[i];
198 }
199
200 return res;
201 }
202
203
204
207 const MPI_Comm comm,
208 const types::global_dof_index total_size)
209 {
210 const unsigned int this_proc = this_mpi_process(comm);
211 const unsigned int n_proc = n_mpi_processes(comm);
212
214 n_proc,
215 total_size);
216 }
217
218
219
220 std::unique_ptr<MPI_Datatype, void (*)(MPI_Datatype *)>
221 create_mpi_data_type_n_bytes(const std::size_t n_bytes)
222 {
223 MPI_Datatype result;
224 int ierr = LargeCount::Type_contiguous_c(n_bytes, MPI_BYTE, &result);
225 AssertThrowMPI(ierr);
226 ierr = MPI_Type_commit(&result);
227 AssertThrowMPI(ierr);
228
229 if constexpr (running_in_debug_mode())
230 {
231 MPI_Count size64;
232 ierr = MPI_Type_size_x(result, &size64);
233 AssertThrowMPI(ierr);
234
235 Assert(size64 == static_cast<MPI_Count>(n_bytes), ExcInternalError());
236 }
237
238 // Now put the new data type into a std::unique_ptr with a custom
239 // deleter. We call the std::unique_ptr constructor that as first
240 // argument takes a pointer (here, a pointer to a copy of the `result`
241 // object, and as second argument a pointer-to-function, for which
242 // we here use a lambda function without captures that acts as the
243 // 'deleter' object: it calls `MPI_Type_free` and then deletes the
244 // pointer. To avoid a compiler warning about a null this pointer
245 // in the lambda (which don't make sense: the lambda doesn't store
246 // anything), we create the deleter first.
247 auto deleter = [](MPI_Datatype *p) {
248 if (p != nullptr)
249 {
250 const int ierr = MPI_Type_free(p);
251 AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
252 delete p;
253 }
254 };
255
256 return std::unique_ptr<MPI_Datatype, void (*)(MPI_Datatype *)>(
257 new MPI_Datatype(result), deleter);
258 }
259
260
261
262 std::vector<unsigned int>
264 const MPI_Comm mpi_comm,
265 const std::vector<unsigned int> &destinations)
266 {
267 const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
268 const unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_comm);
269
270 if constexpr (running_in_debug_mode())
271 {
272 for (const unsigned int destination : destinations)
273 AssertIndexRange(destination, n_procs);
274 }
275
276 // Have a little function that checks if destinations provided
277 // to the current process are unique. The way it does this is
278 // to create a sorted list of destinations and then walk through
279 // the list and look at successive elements -- if we find the
280 // same number twice, we know that the destinations were not
281 // unique
282 const bool my_destinations_are_unique = [destinations]() {
283 if (destinations.empty())
284 return true;
285 else
286 {
287 std::vector<unsigned int> my_destinations = destinations;
288 std::sort(my_destinations.begin(), my_destinations.end());
289 return (std::adjacent_find(my_destinations.begin(),
290 my_destinations.end()) ==
291 my_destinations.end());
292 }
293 }();
294
295 // If all processes report that they have unique destinations,
296 // then we can short-cut the process using a consensus algorithm (which
297 // is implemented only for the case of unique destinations):
298 if (Utilities::MPI::min((my_destinations_are_unique ? 1 : 0), mpi_comm) ==
299 1)
300 {
301 return ConsensusAlgorithms::nbx<char, char>(
302 destinations, {}, {}, {}, mpi_comm);
303 }
304
305 // So we need to run a different algorithm, specifically one that
306 // requires more memory -- MPI_Reduce_scatter_block will require memory
307 // proportional to the number of processes involved; that function is
308 // available for MPI 2.2 or later:
309 static CollectiveMutex mutex;
310 CollectiveMutex::ScopedLock lock(mutex, mpi_comm);
311
312 const int mpi_tag =
314
315 // Calculate the number of messages to send to each process
316 std::vector<unsigned int> dest_vector(n_procs);
317 for (const auto &el : destinations)
318 ++dest_vector[el];
319
320 // Find how many processes will send to this one
321 // by reducing with sum and then scattering the
322 // results over all processes
323 unsigned int n_recv_from;
324 const int ierr = MPI_Reduce_scatter_block(
325 dest_vector.data(), &n_recv_from, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm);
326
327 AssertThrowMPI(ierr);
328
329 // Send myid to every process in `destinations` vector...
330 std::vector<MPI_Request> send_requests(destinations.size());
331 for (const auto &el : destinations)
332 {
333 const int ierr =
334 MPI_Isend(&myid,
335 1,
336 MPI_UNSIGNED,
337 el,
338 mpi_tag,
339 mpi_comm,
340 send_requests.data() + (&el - destinations.data()));
341 AssertThrowMPI(ierr);
342 }
343
344
345 // Receive `n_recv_from` times from the processes
346 // who communicate with this one. Store the obtained id's
347 // in the resulting vector
348 std::vector<unsigned int> origins(n_recv_from);
349 for (auto &el : origins)
350 {
351 const int ierr = MPI_Recv(&el,
352 1,
353 MPI_UNSIGNED,
354 MPI_ANY_SOURCE,
355 mpi_tag,
356 mpi_comm,
357 MPI_STATUS_IGNORE);
358 AssertThrowMPI(ierr);
359 }
360
361 if (destinations.size() > 0)
362 {
363 const int ierr = MPI_Waitall(destinations.size(),
364 send_requests.data(),
365 MPI_STATUSES_IGNORE);
366 AssertThrowMPI(ierr);
367 }
368
369 return origins;
370 }
371
372
373
374 unsigned int
376 const MPI_Comm mpi_comm,
377 const std::vector<unsigned int> &destinations)
378 {
379 // Have a little function that checks if destinations provided
380 // to the current process are unique:
381 const bool my_destinations_are_unique = [destinations]() {
382 std::vector<unsigned int> my_destinations = destinations;
383 const unsigned int n_destinations = my_destinations.size();
384 std::sort(my_destinations.begin(), my_destinations.end());
385 my_destinations.erase(std::unique(my_destinations.begin(),
386 my_destinations.end()),
387 my_destinations.end());
388 return (my_destinations.size() == n_destinations);
389 }();
390
391 // If all processes report that they have unique destinations,
392 // then we can short-cut the process using a consensus algorithm:
393
394 if (Utilities::MPI::min((my_destinations_are_unique ? 1 : 0), mpi_comm) ==
395 1)
396 {
397 return ConsensusAlgorithms::nbx<char, char>(
398 destinations, {}, {}, {}, mpi_comm)
399 .size();
400 }
401 else
402 {
403 const unsigned int n_procs =
405
406 if constexpr (running_in_debug_mode())
407 {
408 for (const unsigned int destination : destinations)
409 {
410 AssertIndexRange(destination, n_procs);
411 Assert(
412 destination != Utilities::MPI::this_mpi_process(mpi_comm),
414 "There is no point in communicating with ourselves."));
415 }
416 }
417
418 // Calculate the number of messages to send to each process
419 std::vector<unsigned int> dest_vector(n_procs);
420 for (const auto &el : destinations)
421 ++dest_vector[el];
422
423 // Find out how many processes will send to this one
424 // MPI_Reduce_scatter(_block) does exactly this
425 unsigned int n_recv_from = 0;
426
427 const int ierr = MPI_Reduce_scatter_block(dest_vector.data(),
428 &n_recv_from,
429 1,
430 MPI_UNSIGNED,
431 MPI_SUM,
432 mpi_comm);
433
434 AssertThrowMPI(ierr);
435
436 return n_recv_from;
437 }
438 }
439
440
441
442 namespace
443 {
444 // custom MIP_Op for calculate_collective_mpi_min_max_avg
445 void
446 max_reduce(const void *in_lhs_,
447 void *inout_rhs_,
448 int *len,
449 MPI_Datatype *)
450 {
451 const MinMaxAvg *in_lhs = static_cast<const MinMaxAvg *>(in_lhs_);
452 MinMaxAvg *inout_rhs = static_cast<MinMaxAvg *>(inout_rhs_);
453
454 for (int i = 0; i < *len; ++i)
455 {
456 inout_rhs[i].sum += in_lhs[i].sum;
457 if (inout_rhs[i].min > in_lhs[i].min)
458 {
459 inout_rhs[i].min = in_lhs[i].min;
460 inout_rhs[i].min_index = in_lhs[i].min_index;
461 }
462 else if (inout_rhs[i].min == in_lhs[i].min)
463 {
464 // choose lower cpu index when tied to make operator commutative
465 if (inout_rhs[i].min_index > in_lhs[i].min_index)
466 inout_rhs[i].min_index = in_lhs[i].min_index;
467 }
468
469 if (inout_rhs[i].max < in_lhs[i].max)
470 {
471 inout_rhs[i].max = in_lhs[i].max;
472 inout_rhs[i].max_index = in_lhs[i].max_index;
473 }
474 else if (inout_rhs[i].max == in_lhs[i].max)
475 {
476 // choose lower cpu index when tied to make operator commutative
477 if (inout_rhs[i].max_index > in_lhs[i].max_index)
478 inout_rhs[i].max_index = in_lhs[i].max_index;
479 }
480 }
481 }
482 } // namespace
483
484
485
486 void
488 const ArrayView<MinMaxAvg> &result,
489 const MPI_Comm mpi_communicator)
490 {
491 // If MPI was not started, we have a serial computation and cannot run
492 // the other MPI commands
493 if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 1)
494 {
495 for (unsigned int i = 0; i < my_values.size(); ++i)
496 {
497 result[i].sum = my_values[i];
498 result[i].avg = my_values[i];
499 result[i].min = my_values[i];
500 result[i].max = my_values[i];
501 result[i].min_index = 0;
502 result[i].max_index = 0;
503 }
504 return;
505 }
506
507 /*
508 * A custom MPI datatype handle describing the memory layout of the
509 * MinMaxAvg struct. Initialized on first pass control reaches the
510 * static variable. So hopefully not initialized too early.
511 */
512 static MPI_Datatype type = []() {
513 MPI_Datatype type;
514
515 int lengths[] = {3, 2, 1};
516
517 MPI_Aint displacements[] = {0,
518 offsetof(MinMaxAvg, min_index),
519 offsetof(MinMaxAvg, avg)};
520
521 MPI_Datatype types[] = {MPI_DOUBLE, MPI_INT, MPI_DOUBLE};
522
523 int ierr =
524 MPI_Type_create_struct(3, lengths, displacements, types, &type);
525 AssertThrowMPI(ierr);
526
527 ierr = MPI_Type_commit(&type);
528 AssertThrowMPI(ierr);
529
530 /* Ensure that we free the allocated datatype again at the end of
531 * the program run just before we call MPI_Finalize():*/
532 InitFinalize::signals.at_mpi_finalize.connect([type]() mutable {
533 int ierr = MPI_Type_free(&type);
534 AssertThrowMPI(ierr);
535 });
536
537 return type;
538 }();
539
540 /*
541 * A custom MPI op handle for our max_reduce function.
542 * Initialized on first pass control reaches the static variable. So
543 * hopefully not initialized too early.
544 */
545 static MPI_Op op = []() {
546 MPI_Op op;
547
548 int ierr =
549 MPI_Op_create(reinterpret_cast<MPI_User_function *>(&max_reduce),
550 static_cast<int>(true),
551 &op);
552 AssertThrowMPI(ierr);
553
554 /* Ensure that we free the allocated op again at the end of the
555 * program run just before we call MPI_Finalize():*/
556 InitFinalize::signals.at_mpi_finalize.connect([op]() mutable {
557 int ierr = MPI_Op_free(&op);
558 AssertThrowMPI(ierr);
559 });
560
561 return op;
562 }();
563
564 AssertDimension(Utilities::MPI::min(my_values.size(), mpi_communicator),
565 Utilities::MPI::max(my_values.size(), mpi_communicator));
566
567 AssertDimension(my_values.size(), result.size());
568
569 // To avoid uninitialized values on some MPI implementations, provide
570 // result with a default value already...
571 MinMaxAvg dummy = {0.,
572 std::numeric_limits<double>::max(),
573 std::numeric_limits<double>::lowest(),
574 0,
575 0,
576 0.};
577
578 for (auto &i : result)
579 i = dummy;
580
581 const unsigned int my_id =
582 ::Utilities::MPI::this_mpi_process(mpi_communicator);
583 const unsigned int numproc =
584 ::Utilities::MPI::n_mpi_processes(mpi_communicator);
585
586 std::vector<MinMaxAvg> in(my_values.size());
587
588 for (unsigned int i = 0; i < my_values.size(); ++i)
589 {
590 in[i].sum = in[i].min = in[i].max = my_values[i];
591 in[i].min_index = in[i].max_index = my_id;
592 }
593
594 int ierr = MPI_Allreduce(
595 in.data(), result.data(), my_values.size(), type, op, mpi_communicator);
596 AssertThrowMPI(ierr);
597
598 for (auto &r : result)
599 r.avg = r.sum / numproc;
600 }
601
602
603#else
604
605 unsigned int
607 {
608 return 1;
609 }
610
611
612
613 unsigned int
615 {
616 return 0;
617 }
618
619
620
621 std::vector<unsigned int>
623 {
624 return std::vector<unsigned int>{0};
625 }
626
627
628
629 std::vector<IndexSet>
631 const MPI_Comm /*comm*/,
633 {
634 return std::vector<IndexSet>(1, complete_index_set(locally_owned_size));
635 }
636
639 const MPI_Comm /*comm*/,
640 const types::global_dof_index total_size)
641 {
642 return complete_index_set(total_size);
643 }
644
645
646
648 duplicate_communicator(const MPI_Comm mpi_communicator)
649 {
650 return mpi_communicator;
651 }
652
653
654
655 void
656 free_communicator(MPI_Comm /*mpi_communicator*/)
657 {}
658
659
660
661 void
662 min_max_avg(const ArrayView<const double> &my_values,
663 const ArrayView<MinMaxAvg> &result,
664 const MPI_Comm)
665 {
666 AssertDimension(my_values.size(), result.size());
667
668 for (unsigned int i = 0; i < my_values.size(); ++i)
669 {
670 result[i].sum = my_values[i];
671 result[i].avg = my_values[i];
672 result[i].min = my_values[i];
673 result[i].max = my_values[i];
674 result[i].min_index = 0;
675 result[i].max_index = 0;
676 }
677 }
678
679#endif
680
681
683 char **&argv,
684 const unsigned int max_num_threads)
685 : InitFinalize(argc,
686 argv,
690 max_num_threads)
691 {}
692
693
694
695 bool
697 {
698#ifdef DEAL_II_WITH_MPI
699 int MPI_has_been_started = 0;
700 const int ierr = MPI_Initialized(&MPI_has_been_started);
701 AssertThrowMPI(ierr);
702
703 return (MPI_has_been_started > 0);
704#else
705 return false;
706#endif
707 }
708
709
710
711 namespace
712 {
717 namespace ComputeIndexOwner
718 {
719 class FlexibleIndexStorage
720 {
721 public:
722 using index_type = unsigned int;
723 static const index_type invalid_index_value =
725
726 FlexibleIndexStorage(const bool use_vector = true);
727
728 void
729 reinit(const bool use_vector,
730 const bool index_range_contiguous,
731 const std::size_t size);
732
733 void
734 fill(const std::size_t start,
735 const std::size_t end,
736 const index_type &value);
737
738 index_type &
739 operator[](const std::size_t index);
740
741 index_type
742 operator[](const std::size_t index) const;
743
744 bool
745 entry_has_been_set(const std::size_t index) const;
746
747 private:
749 std::size_t size;
750 std::vector<index_type> data;
751 std::map<std::size_t, index_type> data_map;
752 };
753
754
755
761 struct Dictionary
762 {
784 static constexpr unsigned int range_minimum_grain_size = 64;
785
792 static constexpr unsigned int sparsity_factor = 4;
793
794
800 Dictionary(const IndexSet &owned_indices, const MPI_Comm comm);
801
808 FlexibleIndexStorage actually_owning_ranks;
809
814 std::vector<unsigned int> actually_owning_rank_list;
815
823
829 std::pair<types::global_dof_index, types::global_dof_index>
831
838
843
849
855 unsigned int stride_small_size;
856
862 unsigned int
863 dof_to_dict_rank(const types::global_dof_index i);
864
870 get_index_offset(const unsigned int rank);
871
877 unsigned int
878 get_owning_rank_index(const unsigned int rank_in_owned_indices,
879 const unsigned int guess = 0);
880
881 private:
886 void
887 partition(const IndexSet &owned_indices, const MPI_Comm comm);
888 };
889
890
891
898 class ConsensusAlgorithmsPayload
899 {
900 public:
901 using RequestType = std::vector<
902 std::pair<types::global_dof_index, types::global_dof_index>>;
903 using AnswerType = std::vector<unsigned int>;
904
908 ConsensusAlgorithmsPayload(const IndexSet &owned_indices,
909 const IndexSet &indices_to_look_up,
910 const MPI_Comm comm,
911 std::vector<unsigned int> &owning_ranks,
912 const bool track_index_requesters = false);
913
918
924
929
933 const unsigned int my_rank;
934
939 const unsigned int n_procs;
940
948
954 std::vector<unsigned int> &owning_ranks;
955
959 Dictionary dict;
960
970 std::vector<std::vector<
971 std::pair<unsigned int,
972 std::vector<std::pair<types::global_dof_index,
975
981 std::map<unsigned int,
982 std::pair<std::vector<types::global_dof_index>,
983 std::vector<unsigned int>>>
985
989 std::vector<unsigned int>
990 compute_targets();
991
995 std::vector<
996 std::pair<types::global_dof_index, types::global_dof_index>>
997 create_request(const unsigned int other_rank);
998
1002 std::vector<unsigned int>
1003 answer_request(
1004 const unsigned int other_rank,
1005 const std::vector<std::pair<types::global_dof_index,
1006 types::global_dof_index>> &buffer_recv);
1007
1012 void
1013 process_answer(const unsigned int other_rank,
1014 const std::vector<unsigned int> &recv_buffer);
1015
1031 std::map<unsigned int, IndexSet>
1032 get_requesters();
1033
1034 private:
1048 void
1049 append_index_origin(
1050 const types::global_dof_index index_within_dictionary,
1051 const unsigned int rank_of_request,
1052 const unsigned int rank_of_owner,
1053 unsigned int &owner_index_guess);
1054 };
1055
1056 /* ------------------------- inline functions ----------------------- */
1057
1058 inline unsigned int
1059 Dictionary::dof_to_dict_rank(const types::global_dof_index i)
1060 {
1061 // note: this formula is also explicitly used in
1062 // get_index_offset(), so keep the two in sync
1063 return (i / dofs_per_process) * stride_small_size;
1064 }
1065
1066
1068 Dictionary::get_index_offset(const unsigned int rank)
1069 {
1070 return std::min(dofs_per_process *
1071 static_cast<types::global_dof_index>(
1072 (rank + stride_small_size - 1) /
1074 size);
1075 }
1076
1077
1078
1079 inline unsigned int
1080 Dictionary::get_owning_rank_index(
1081 const unsigned int rank_in_owned_indices,
1082 const unsigned int guess)
1083 {
1085 if (actually_owning_rank_list[guess] == rank_in_owned_indices)
1086 return guess;
1087 else
1088 {
1089 auto it = std::lower_bound(actually_owning_rank_list.begin(),
1091 rank_in_owned_indices);
1093 Assert(*it == rank_in_owned_indices, ExcInternalError());
1094 return it - actually_owning_rank_list.begin();
1095 }
1096 }
1097
1098
1099 const FlexibleIndexStorage::index_type
1100 FlexibleIndexStorage::invalid_index_value;
1101
1102
1103
1104 FlexibleIndexStorage::FlexibleIndexStorage(const bool use_vector)
1106 , size(0)
1107 {}
1108
1109
1110
1111 void
1112 FlexibleIndexStorage::reinit(const bool use_vector,
1113 const bool index_range_contiguous,
1114 const std::size_t size)
1115 {
1116 this->use_vector = use_vector;
1117 this->size = size;
1118
1119 data = {};
1120 data_map.clear();
1121
1122 // in case we have contiguous indices, only fill the vector upon
1123 // first request in `fill`
1124 if (use_vector && !index_range_contiguous)
1125 data.resize(size, invalid_index_value);
1126 }
1127
1128
1129
1130 void
1131 FlexibleIndexStorage::fill(
1132 const std::size_t start,
1133 const std::size_t end,
1134 const FlexibleIndexStorage::index_type &value)
1135 {
1136 AssertIndexRange(start, size);
1137 AssertIndexRange(end, size + 1);
1138
1139 if (use_vector)
1140 {
1141 if (data.empty() && end > start)
1142 {
1143 // in debug mode, we want to track whether we set all
1144 // indices, so we first fill an invalid index and only later
1145 // the actual ones, whereas we simply assign the given rank
1146 // to the complete vector the first time we pass around in
1147 // this function in release mode to avoid touching data
1148 // unnecessarily (and overwrite the smaller pieces), as the
1149 // locally owned part comes first
1150 if constexpr (running_in_debug_mode())
1151 {
1152 data.resize(size, invalid_index_value);
1153 std::fill(data.begin() + start,
1154 data.begin() + end,
1155 value);
1156 }
1157 else
1158 {
1159 data.resize(size, value);
1160 }
1161 }
1162 else
1163 {
1164 AssertDimension(data.size(), size);
1165 std::fill(data.begin() + start, data.begin() + end, value);
1166 }
1167 }
1168 else
1169 {
1170 for (auto i = start; i < end; ++i)
1171 data_map[i] = value;
1172 }
1173 }
1174
1175
1176
1177 FlexibleIndexStorage::index_type &
1178 FlexibleIndexStorage::operator[](const std::size_t index)
1179 {
1180 AssertIndexRange(index, size);
1181
1182 if (use_vector)
1183 {
1184 AssertDimension(data.size(), size);
1185 return data[index];
1186 }
1187 else
1188 {
1189 return data_map.try_emplace(index, invalid_index_value)
1190 .first->second;
1191 }
1192 }
1193
1194
1195
1196 inline bool
1197 FlexibleIndexStorage::entry_has_been_set(const std::size_t index) const
1198 {
1199 AssertIndexRange(index, size);
1200
1201 if (use_vector)
1202 {
1203 if (data.empty())
1204 return false;
1205
1206 AssertDimension(data.size(), size);
1207 return data[index] != invalid_index_value;
1208 }
1209 else
1210 return data_map.find(index) != data_map.end();
1211 }
1212
1213
1214
1215 Dictionary::Dictionary(const IndexSet &owned_indices,
1216 const MPI_Comm comm)
1217 {
1218 // 1) set up the partition
1220
1221 unsigned int my_rank = this_mpi_process(comm);
1222
1223 types::global_dof_index dic_local_received = 0;
1224 std::map<unsigned int,
1225 std::vector<std::pair<types::global_dof_index,
1227 buffers;
1228
1229 const auto owned_indices_size_actual =
1231
1232 actually_owning_ranks.reinit((owned_indices_size_actual *
1234 owned_indices_size_actual ==
1237
1238 // 2) collect relevant processes and process local dict entries
1239 for (auto interval = owned_indices.begin_intervals();
1240 interval != owned_indices.end_intervals();
1241 ++interval)
1242 {
1243 // Due to the granularity of the dictionary, the interval
1244 // might be split into several ranges of processor owner
1245 // ranks. Here, we process the interval by breaking into
1246 // smaller pieces in terms of the dictionary number.
1247 std::pair<types::global_dof_index, types::global_dof_index>
1248 index_range(*interval->begin(), interval->last() + 1);
1249
1250 AssertThrow(index_range.second <= size, ExcInternalError());
1251
1252 while (index_range.first != index_range.second)
1253 {
1254 Assert(index_range.first < index_range.second,
1256
1257 const unsigned int owner =
1258 dof_to_dict_rank(index_range.first);
1259
1260 // this explicitly picks up the formula of
1261 // dof_to_dict_rank, so the two places must be in sync
1262 const types::global_dof_index next_index =
1263 std::min(get_index_offset(owner + 1), index_range.second);
1264
1265 Assert(next_index > index_range.first, ExcInternalError());
1266
1267 if constexpr (running_in_debug_mode())
1268 {
1269 // make sure that the owner is the same on the current
1270 // interval
1271 for (types::global_dof_index i = index_range.first + 1;
1272 i < next_index;
1273 ++i)
1274 AssertDimension(owner, dof_to_dict_rank(i));
1275 }
1276
1277 // add the interval, either to the local range or into a
1278 // buffer to be sent to another processor
1279 if (owner == my_rank)
1280 {
1281 actually_owning_ranks.fill(index_range.first -
1282 local_range.first,
1283 next_index - local_range.first,
1284 my_rank);
1285 dic_local_received += next_index - index_range.first;
1286 if (actually_owning_rank_list.empty())
1288 }
1289 else
1290 buffers[owner].emplace_back(index_range.first, next_index);
1291
1292 index_range.first = next_index;
1293 }
1294 }
1295
1296#ifdef DEAL_II_WITH_MPI
1297 n_dict_procs_in_owned_indices = buffers.size();
1298 std::vector<MPI_Request> request;
1299
1300 // Check if index set space is partitioned globally without gaps.
1301 if (owned_indices_size_actual == owned_indices.size())
1302 {
1303 // no gaps: setup is simple! Processes send their locally owned
1304 // indices to the dictionary. The dictionary stores the sending
1305 // rank for each index. The dictionary knows exactly
1306 // when it is set up when all indices it is responsible for
1307 // have been processed.
1308
1309 request.reserve(n_dict_procs_in_owned_indices);
1310
1311 // protect the following communication steps using a mutex:
1312 static CollectiveMutex mutex;
1313 CollectiveMutex::ScopedLock lock(mutex, comm);
1314
1315 const int mpi_tag =
1317
1318
1319 // 3) send messages with local dofs to the right dict process
1320 for (const auto &rank_pair : buffers)
1321 {
1322 request.push_back(MPI_Request());
1323 const int ierr =
1324 MPI_Isend(rank_pair.second.data(),
1325 rank_pair.second.size() * 2,
1328 rank_pair.first,
1329 mpi_tag,
1330 comm,
1331 &request.back());
1332 AssertThrowMPI(ierr);
1333 }
1334
1335 // 4) receive messages until all dofs in dict are processed
1336 while (this->locally_owned_size != dic_local_received)
1337 {
1338 // wait for an incoming message
1339 MPI_Status status;
1340 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, comm, &status);
1341 AssertThrowMPI(ierr);
1342
1343 // retrieve size of incoming message
1344 int number_amount;
1345 ierr = MPI_Get_count(&status,
1348 &number_amount);
1349 AssertThrowMPI(ierr);
1350
1351 const auto other_rank = status.MPI_SOURCE;
1352 actually_owning_rank_list.push_back(other_rank);
1353
1354 // receive message
1355 Assert(number_amount % 2 == 0, ExcInternalError());
1356 std::vector<
1357 std::pair<types::global_dof_index, types::global_dof_index>>
1358 buffer(number_amount / 2);
1359 ierr = MPI_Recv(buffer.data(),
1360 number_amount,
1363 status.MPI_SOURCE,
1364 status.MPI_TAG,
1365 comm,
1366 MPI_STATUS_IGNORE);
1367 AssertThrowMPI(ierr);
1368 // process message: loop over all intervals
1369 for (auto interval : buffer)
1370 {
1371 if constexpr (library_build_mode ==
1373 {
1374 for (types::global_dof_index i = interval.first;
1375 i < interval.second;
1376 i++)
1377 Assert(actually_owning_ranks.entry_has_been_set(
1378 i - local_range.first) == false,
1380 Assert(interval.first >= local_range.first &&
1381 interval.first < local_range.second,
1383 Assert(interval.second > local_range.first &&
1384 interval.second <= local_range.second,
1386 }
1387
1388 actually_owning_ranks.fill(interval.first -
1389 local_range.first,
1390 interval.second -
1391 local_range.first,
1392 other_rank);
1393 dic_local_received += interval.second - interval.first;
1394 }
1395 }
1396 }
1397 else
1398 {
1399 // with gap: use a ConsensusAlgorithm to determine when all
1400 // dictionaries have been set up.
1401
1402 // 3/4) use a ConsensusAlgorithm to send messages with local
1403 // dofs to the right dict process
1404
1405 using RequestType = std::vector<
1406 std::pair<types::global_dof_index, types::global_dof_index>>;
1407
1408 ConsensusAlgorithms::selector<RequestType>(
1409 /* targets = */
1410 [&buffers]() {
1411 std::vector<unsigned int> targets;
1412 targets.reserve(buffers.size());
1413 for (const auto &rank_pair : buffers)
1414 targets.emplace_back(rank_pair.first);
1415
1416 return targets;
1417 }(),
1418
1419 /* create_request = */
1420 [&buffers](const unsigned int target_rank) -> RequestType {
1421 return buffers.at(target_rank);
1422 },
1423
1424 /* process_request = */
1425 [&](const unsigned int source_rank,
1426 const RequestType &request) -> void {
1427 // process message: loop over all intervals
1428 for (auto interval : request)
1429 {
1430 if constexpr (library_build_mode ==
1432 {
1433 for (types::global_dof_index i = interval.first;
1434 i < interval.second;
1435 i++)
1436 Assert(
1437 actually_owning_ranks.entry_has_been_set(
1438 i - local_range.first) == false,
1439 ExcMessage(
1440 "Multiple processes seem to own the same global index. "
1441 "A possible reason is that the sets of locally owned "
1442 "indices are not distinct."));
1443 Assert(interval.first < interval.second,
1445 Assert(
1446 local_range.first <= interval.first &&
1447 interval.second <= local_range.second,
1448 ExcMessage(
1449 "The specified interval is not handled by the current process."));
1450 }
1451 actually_owning_ranks.fill(interval.first -
1452 local_range.first,
1453 interval.second -
1454 local_range.first,
1455 source_rank);
1456 }
1457 actually_owning_rank_list.push_back(source_rank);
1458 },
1459
1460 comm);
1461 }
1462
1463 std::sort(actually_owning_rank_list.begin(),
1465
1466 for (unsigned int i = 1; i < actually_owning_rank_list.size(); ++i)
1470
1471 // 5) make sure that all messages have been sent
1472 if (request.size() > 0)
1473 {
1474 const int ierr = MPI_Waitall(request.size(),
1475 request.data(),
1476 MPI_STATUSES_IGNORE);
1477 AssertThrowMPI(ierr);
1478 }
1479
1480#else
1481 Assert(buffers.empty(), ExcInternalError());
1482 (void)comm;
1483 (void)dic_local_received;
1484#endif
1485 }
1486
1487
1488
1489 void
1490 Dictionary::partition(const IndexSet &owned_indices,
1491 const MPI_Comm comm)
1492 {
1493 const unsigned int n_procs = n_mpi_processes(comm);
1494 const unsigned int my_rank = this_mpi_process(comm);
1495
1497
1499
1501 std::max<types::global_dof_index>((size + n_procs - 1) / n_procs,
1503
1505 std::max<unsigned int>(dofs_per_process * n_procs / size, 1);
1506
1507 local_range.first = get_index_offset(my_rank);
1508 local_range.second = get_index_offset(my_rank + 1);
1509
1511 }
1512
1513
1514 ConsensusAlgorithmsPayload::ConsensusAlgorithmsPayload(
1515 const IndexSet &owned_indices,
1517 const MPI_Comm comm,
1518 std::vector<unsigned int> &owning_ranks,
1519 const bool track_index_requesters)
1522 , comm(comm)
1529 {}
1530
1531
1532
1533 std::vector<unsigned int>
1534 ConsensusAlgorithmsPayload::compute_targets()
1535 {
1536 std::vector<unsigned int> targets;
1537
1539 unsigned int index = 0;
1540 unsigned int owner_index_guess = 0;
1541 for (auto i : indices_to_look_up)
1542 {
1543 unsigned int other_rank = dict.dof_to_dict_rank(i);
1544 if (other_rank == my_rank)
1545 {
1547 dict.actually_owning_ranks[i - dict.local_range.first];
1549 append_index_origin(i - dict.local_range.first,
1550 my_rank,
1551 owning_ranks[index],
1552 owner_index_guess);
1553 }
1554 else
1555 {
1556 if (targets.empty() || targets.back() != other_rank)
1557 targets.push_back(other_rank);
1558 auto &indices = indices_to_look_up_by_dict_rank[other_rank];
1559 indices.first.push_back(i);
1560 indices.second.push_back(index);
1561 }
1562 ++index;
1563 }
1564
1565 Assert(targets.size() == indices_to_look_up_by_dict_rank.size(),
1566 ExcMessage("Size does not match!"));
1567
1568 return targets;
1569 }
1570
1571
1572
1573 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1574 ConsensusAlgorithmsPayload::create_request(
1575 const unsigned int other_rank)
1576 {
1577 std::vector<
1578 std::pair<types::global_dof_index, types::global_dof_index>>
1579 send_buffer;
1580
1581 // create index set and compress data to be sent
1582 auto &indices_i = indices_to_look_up_by_dict_rank[other_rank].first;
1583 IndexSet is(dict.size);
1584 is.add_indices(indices_i.begin(), indices_i.end());
1585 is.compress();
1586
1587 for (auto interval = is.begin_intervals();
1588 interval != is.end_intervals();
1589 ++interval)
1590 send_buffer.emplace_back(*interval->begin(), interval->last() + 1);
1591
1592 return send_buffer;
1593 }
1594
1595
1596
1597 std::vector<unsigned int>
1598 ConsensusAlgorithmsPayload::answer_request(
1599 const unsigned int other_rank,
1600 const std::vector<std::pair<types::global_dof_index,
1601 types::global_dof_index>> &buffer_recv)
1602 {
1603 std::vector<unsigned int> request_buffer;
1604
1605 unsigned int owner_index_guess = 0;
1606 for (const auto &interval : buffer_recv)
1607 for (auto i = interval.first; i < interval.second; ++i)
1608 {
1609 const unsigned int actual_owner =
1610 dict.actually_owning_ranks[i - dict.local_range.first];
1611 request_buffer.push_back(actual_owner);
1612
1614 append_index_origin(i - dict.local_range.first,
1615 other_rank,
1616 actual_owner,
1617 owner_index_guess);
1618 }
1619
1620 return request_buffer;
1621 }
1622
1623
1624
1625 void
1626 ConsensusAlgorithmsPayload::process_answer(
1627 const unsigned int other_rank,
1628 const std::vector<unsigned int> &recv_buffer)
1629 {
1630 const auto &recv_indices =
1631 indices_to_look_up_by_dict_rank[other_rank].second;
1632 AssertDimension(recv_indices.size(), recv_buffer.size());
1633 for (unsigned int j = 0; j < recv_indices.size(); ++j)
1634 owning_ranks[recv_indices[j]] = recv_buffer[j];
1635 }
1636
1637
1638
1639 std::map<unsigned int, IndexSet>
1640 ConsensusAlgorithmsPayload::get_requesters()
1641 {
1643 ExcMessage("Must enable index range tracking in "
1644 "constructor of ConsensusAlgorithmProcess"));
1645
1646 std::map<unsigned int, ::IndexSet> requested_indices;
1647
1648#ifdef DEAL_II_WITH_MPI
1649
1650 static CollectiveMutex mutex;
1651 CollectiveMutex::ScopedLock lock(mutex, comm);
1652
1653 const int mpi_tag = Utilities::MPI::internal::Tags::
1655
1656 // reserve enough slots for the requests ahead; depending on
1657 // whether the owning rank is one of the requesters or not, we
1658 // might have one less requests to execute, so fill the requests
1659 // on demand.
1660 std::vector<MPI_Request> send_requests;
1661 send_requests.reserve(requesters.size());
1662
1663 // We use an integer vector for the data exchange. Since we send
1664 // data associated to intervals with different requesters, we will
1665 // need to send (a) the MPI rank of the requester, (b) the number
1666 // of intervals directed to this requester, and (c) a list of
1667 // intervals, i.e., two integers per interval. The number of items
1668 // sent in total can be deduced both via the MPI status message at
1669 // the receiver site as well as be counting the buckets from
1670 // different requesters.
1671 std::vector<std::vector<types::global_dof_index>> send_data(
1672 requesters.size());
1673 for (unsigned int i = 0; i < requesters.size(); ++i)
1674 {
1675 // special code for our own indices
1676 if (dict.actually_owning_rank_list[i] == my_rank)
1677 {
1678 for (const auto &j : requesters[i])
1679 {
1680 const types::global_dof_index index_offset =
1681 dict.get_index_offset(my_rank);
1682 IndexSet &my_index_set = requested_indices[j.first];
1683 my_index_set.set_size(owned_indices.size());
1684 for (const auto &interval : j.second)
1685 my_index_set.add_range(index_offset + interval.first,
1686 index_offset + interval.second);
1687 }
1688 }
1689 else
1690 {
1691 for (const auto &j : requesters[i])
1692 {
1693 send_data[i].push_back(j.first);
1694 send_data[i].push_back(j.second.size());
1695 for (const auto &interval : j.second)
1696 {
1697 send_data[i].push_back(interval.first);
1698 send_data[i].push_back(interval.second);
1699 }
1700 }
1701 send_requests.push_back(MPI_Request());
1702 const int ierr =
1703 MPI_Isend(send_data[i].data(),
1704 send_data[i].size(),
1707 dict.actually_owning_rank_list[i],
1708 mpi_tag,
1709 comm,
1710 &send_requests.back());
1711 AssertThrowMPI(ierr);
1712 }
1713 }
1714
1715 // receive the data
1716 for (unsigned int c = 0; c < dict.n_dict_procs_in_owned_indices; ++c)
1717 {
1718 // wait for an incoming message
1719 MPI_Status status;
1720 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, comm, &status);
1721 AssertThrowMPI(ierr);
1722
1723 // retrieve size of incoming message
1724 int number_amount;
1725 ierr = MPI_Get_count(
1726 &status,
1727 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1728 &number_amount);
1729 AssertThrowMPI(ierr);
1730
1731 // receive message
1732 Assert(number_amount % 2 == 0, ExcInternalError());
1733 std::vector<
1734 std::pair<types::global_dof_index, types::global_dof_index>>
1735 buffer(number_amount / 2);
1736 ierr = MPI_Recv(
1737 buffer.data(),
1738 number_amount,
1739 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1740 status.MPI_SOURCE,
1741 status.MPI_TAG,
1742 comm,
1743 &status);
1744 AssertThrowMPI(ierr);
1745
1746 // unpack the message and translate the dictionary-local
1747 // indices coming via MPI to the global index range
1748 const types::global_dof_index index_offset =
1749 dict.get_index_offset(status.MPI_SOURCE);
1750 unsigned int offset = 0;
1751 while (offset < buffer.size())
1752 {
1753 AssertIndexRange(offset + buffer[offset].second,
1754 buffer.size());
1755
1756 IndexSet my_index_set(owned_indices.size());
1757 for (unsigned int i = offset + 1;
1758 i < offset + buffer[offset].second + 1;
1759 ++i)
1760 my_index_set.add_range(index_offset + buffer[i].first,
1761 index_offset + buffer[i].second);
1762
1763 // the underlying index set is able to merge ranges coming
1764 // from different ranks due to the partitioning in the
1765 // dictionary
1766 IndexSet &index_set = requested_indices[buffer[offset].first];
1767 if (index_set.size() == 0)
1768 index_set.set_size(owned_indices.size());
1769 index_set.add_indices(my_index_set);
1770
1771 offset += buffer[offset].second + 1;
1772 }
1773 AssertDimension(offset, buffer.size());
1774 }
1775
1776 if (send_requests.size() > 0)
1777 {
1778 const auto ierr = MPI_Waitall(send_requests.size(),
1779 send_requests.data(),
1780 MPI_STATUSES_IGNORE);
1781 AssertThrowMPI(ierr);
1782 }
1783
1784
1785 if constexpr (running_in_debug_mode())
1786 {
1787 for (const auto &it : requested_indices)
1788 {
1789 IndexSet copy_set = it.second;
1790 copy_set.subtract_set(owned_indices);
1791 Assert(copy_set.n_elements() == 0,
1793 "The indices requested from the current "
1794 "MPI rank should be locally owned here!"));
1795 }
1796 }
1797
1798#endif // DEAL_II_WITH_MPI
1799
1800 return requested_indices;
1801 }
1802
1803
1804
1805 void
1806 ConsensusAlgorithmsPayload::append_index_origin(
1807 const types::global_dof_index index_within_dict,
1808 const unsigned int rank_of_request,
1809 const unsigned int rank_of_owner,
1810 unsigned int &owner_index_guess)
1811 {
1812 // remember who requested which index. We want to use an
1813 // std::vector with simple addressing, via a good guess from the
1814 // preceding index, rather than std::map, because this is an inner
1815 // loop and it avoids the map lookup in every iteration
1816 owner_index_guess =
1817 dict.get_owning_rank_index(rank_of_owner, owner_index_guess);
1818
1819 auto &request = requesters[owner_index_guess];
1820 if (request.empty() || request.back().first != rank_of_request)
1821 request.emplace_back(
1822 rank_of_request,
1823 std::vector<
1824 std::pair<types::global_dof_index, types::global_dof_index>>());
1825
1826 auto &intervals = request.back().second;
1827 if (intervals.empty() || intervals.back().second != index_within_dict)
1828 intervals.emplace_back(index_within_dict, index_within_dict + 1);
1829 else
1830 ++intervals.back().second;
1831 }
1832
1833 } // namespace ComputeIndexOwner
1834 } // namespace
1835
1836
1837
1838 std::vector<unsigned int>
1841 const MPI_Comm comm)
1842 {
1844 ExcMessage("IndexSets have to have the same sizes."));
1845
1846 Assert(
1848 ExcMessage("IndexSets have to have the same size on all processes."));
1849
1850 std::vector<unsigned int> owning_ranks(indices_to_look_up.n_elements());
1851
1852 // Step 1: setup dictionary
1853 // The input owned_indices can be partitioned arbitrarily. In the
1854 // dictionary, the index set is statically repartitioned among the
1855 // processes again and extended with information with the actual owner
1856 // of that the index.
1857 ComputeIndexOwner::ConsensusAlgorithmsPayload process(
1860 comm,
1862 /* keep track of requesters = */ false);
1863
1864 // Step 2: read dictionary
1865 // Communicate with the process who owns the index in the static
1866 // partition (i.e. in the dictionary). This process returns the actual
1867 // owner of the index.
1868 using RequestType =
1869 ComputeIndexOwner::ConsensusAlgorithmsPayload::RequestType;
1870 using AnswerType =
1871 ComputeIndexOwner::ConsensusAlgorithmsPayload::AnswerType;
1872 ConsensusAlgorithms::selector<RequestType, AnswerType>(
1873 process.compute_targets(),
1874 [&process](const unsigned int other_rank) -> RequestType {
1875 return process.create_request(other_rank);
1876 },
1877 [&process](const unsigned int other_rank, const RequestType &r)
1878 -> AnswerType { return process.answer_request(other_rank, r); },
1879 [&process](const unsigned int other_rank, const AnswerType &a) -> void {
1880 process.process_answer(other_rank, a);
1881 },
1882 comm);
1883
1884 return owning_ranks;
1885 }
1886
1887
1888
1889 std::pair<std::vector<unsigned int>, std::map<unsigned int, IndexSet>>
1892 const MPI_Comm &comm)
1893 {
1895 ExcMessage("IndexSets have to have the same sizes."));
1896
1897 Assert(
1899 ExcMessage("IndexSets have to have the same size on all processes."));
1900
1901 std::vector<unsigned int> owning_ranks(indices_to_look_up.n_elements());
1902
1903 // Step 1: setup dictionary
1904 // The input owned_indices can be partitioned arbitrarily. In the
1905 // dictionary, the index set is statically repartitioned among the
1906 // processes again and extended with information with the actual owner
1907 // of that the index.
1908 ComputeIndexOwner::ConsensusAlgorithmsPayload process(
1910
1911 // Step 2: read dictionary
1912 // Communicate with the process who owns the index in the static
1913 // partition (i.e. in the dictionary). This process returns the actual
1914 // owner of the index.
1915 using RequestType =
1916 ComputeIndexOwner::ConsensusAlgorithmsPayload::RequestType;
1917 using AnswerType =
1918 ComputeIndexOwner::ConsensusAlgorithmsPayload::AnswerType;
1919 ConsensusAlgorithms::selector<RequestType, AnswerType>(
1920 process.compute_targets(),
1921 [&process](const unsigned int other_rank) -> RequestType {
1922 return process.create_request(other_rank);
1923 },
1924 [&process](const unsigned int other_rank, const RequestType &r)
1925 -> AnswerType { return process.answer_request(other_rank, r); },
1926 [&process](const unsigned int other_rank, const AnswerType &a) -> void {
1927 process.process_answer(other_rank, a);
1928 },
1929 comm);
1930
1931 return {owning_ranks, process.get_requesters()};
1932 }
1933
1934
1935
1936 namespace internal
1937 {
1938 namespace CollectiveMutexImplementation
1939 {
1944 void
1945 check_exception()
1946 {
1947#ifdef DEAL_II_WITH_MPI
1948 if (std::uncaught_exceptions() > 0)
1949 {
1950 std::cerr
1951 << "---------------------------------------------------------\n"
1952 << "An exception was thrown inside a section of the program\n"
1953 << "guarded by a CollectiveMutex.\n"
1954 << "Because a CollectiveMutex guards critical communication\n"
1955 << "handling the exception would likely\n"
1956 << "deadlock because only the current process is aware of the\n"
1957 << "exception. To prevent this deadlock, the program will be\n"
1958 << "aborted.\n"
1959 << "---------------------------------------------------------"
1960 << std::endl;
1961
1962 MPI_Abort(MPI_COMM_WORLD, 1);
1963 }
1964#endif
1965 }
1966 } // namespace CollectiveMutexImplementation
1967 } // namespace internal
1968
1969
1970
1971 CollectiveMutex::CollectiveMutex()
1972 : locked(false)
1973 , request(MPI_REQUEST_NULL)
1974 {
1976 }
1977
1978
1979
1981 {
1982 // First check if this destructor is called during exception handling
1983 // if so, abort.
1985
1986 Assert(
1987 !locked,
1988 ExcMessage(
1989 "Error: MPI::CollectiveMutex is still locked while being destroyed!"));
1990
1992 }
1993
1994
1995
1996 void
1998 {
1999 Assert(
2000 !locked,
2001 ExcMessage(
2002 "Error: MPI::CollectiveMutex needs to be unlocked before lock()"));
2003
2004#ifdef DEAL_II_WITH_MPI
2005
2006 if (job_supports_mpi())
2007 {
2008 // TODO: For now, we implement this mutex with a blocking barrier in
2009 // the lock and unlock. It needs to be tested, if we can move to a
2010 // nonblocking barrier (code disabled below).
2011
2012 const int ierr = MPI_Barrier(comm);
2013 AssertThrowMPI(ierr);
2014
2015# if 0
2016 // wait for non-blocking barrier to finish. This is a noop the
2017 // first time we lock().
2018 const int ierr = MPI_Wait(&request, MPI_STATUS_IGNORE);
2019 AssertThrowMPI(ierr);
2020# else
2021 // nothing to do as blocking barrier already completed
2022# endif
2023 }
2024#else
2025 (void)comm;
2026#endif
2027
2028 locked = true;
2029 }
2030
2031
2032
2033 void
2035 {
2036 // First check if this function is called during exception handling
2037 // if so, abort. This can happen if a ScopedLock is destroyed.
2039
2040 Assert(
2041 locked,
2042 ExcMessage(
2043 "Error: MPI::CollectiveMutex needs to be locked before unlock()"));
2044
2045#ifdef DEAL_II_WITH_MPI
2046
2047 if (job_supports_mpi())
2048 {
2049 // TODO: For now, we implement this mutex with a blocking barrier
2050 // in the lock and unlock. It needs to be tested, if we can move
2051 // to a nonblocking barrier (code disabled below):
2052# if 0
2053 const int ierr = MPI_Ibarrier(comm, &request);
2054 AssertThrowMPI(ierr);
2055# else
2056 const int ierr = MPI_Barrier(comm);
2057 AssertThrowMPI(ierr);
2058# endif
2059 }
2060#else
2061 (void)comm;
2062#endif
2063
2064 locked = false;
2065 }
2066
2067
2068#ifndef DOXYGEN
2069 // explicit instantiations
2070
2071 // booleans aren't in MPI_SCALARS
2072 template bool
2073 reduce(const bool &,
2074 const MPI_Comm,
2075 const std::function<bool(const bool &, const bool &)> &,
2076 const unsigned int);
2077
2078 template std::vector<bool>
2079 reduce(const std::vector<bool> &,
2080 const MPI_Comm,
2081 const std::function<std::vector<bool>(const std::vector<bool> &,
2082 const std::vector<bool> &)> &,
2083 const unsigned int);
2084
2085 template bool
2086 all_reduce(const bool &,
2087 const MPI_Comm,
2088 const std::function<bool(const bool &, const bool &)> &);
2089
2090 template std::vector<bool>
2091 all_reduce(
2092 const std::vector<bool> &,
2093 const MPI_Comm,
2094 const std::function<std::vector<bool>(const std::vector<bool> &,
2095 const std::vector<bool> &)> &);
2096
2097 // We need an explicit instantiation of this for the same reason as the
2098 // other types described in mpi.inst.in
2099 template void
2100 internal::all_reduce<bool>(const MPI_Op &,
2101 const ArrayView<const bool> &,
2102 const MPI_Comm,
2103 const ArrayView<bool> &);
2104
2105
2106 template bool
2107 logical_or<bool>(const bool &, const MPI_Comm);
2108
2109
2110 template void
2111 logical_or<bool>(const ArrayView<const bool> &,
2112 const MPI_Comm,
2113 const ArrayView<bool> &);
2114
2115
2116 template std::vector<unsigned int>
2117 compute_set_union(const std::vector<unsigned int> &vec,
2118 const MPI_Comm comm);
2119
2120
2121 template std::set<unsigned int>
2122 compute_set_union(const std::set<unsigned int> &set, const MPI_Comm comm);
2123#endif
2124
2125#include "base/mpi.inst"
2126 } // end of namespace MPI
2127} // end of namespace Utilities
2128
value_type * data() const noexcept
Definition array_view.h:666
std::size_t size() const
Definition array_view.h:689
IntervalIterator end_intervals() const
Definition index_set.h:1754
size_type size() const
Definition index_set.h:1787
size_type n_elements() const
Definition index_set.h:1945
void set_size(const size_type size)
Definition index_set.h:1775
IntervalIterator begin_intervals() const
Definition index_set.h:1742
void subtract_set(const IndexSet &other)
Definition index_set.cc:478
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1814
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
static void unregister_request(MPI_Request &request)
static void register_request(MPI_Request &request)
static Signals signals
void lock(const MPI_Comm comm)
Definition mpi.cc:1997
void unlock(const MPI_Comm comm)
Definition mpi.cc:2034
MPI_InitFinalize(int &argc, char **&argv, const unsigned int max_num_threads=numbers::invalid_unsigned_int)
Definition mpi.cc:682
constexpr LibraryBuildMode library_build_mode
Definition config.h:63
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition config.h:598
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition config.h:641
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1215
InitializeLibrary
const unsigned int my_rank
Definition mpi.cc:933
const IndexSet & indices_to_look_up
Definition mpi.cc:923
std::vector< index_type > data
Definition mpi.cc:750
std::pair< types::global_dof_index, types::global_dof_index > local_range
Definition mpi.cc:830
std::vector< unsigned int > actually_owning_rank_list
Definition mpi.cc:814
static constexpr unsigned int range_minimum_grain_size
Definition mpi.cc:784
static constexpr unsigned int sparsity_factor
Definition mpi.cc:792
unsigned int stride_small_size
Definition mpi.cc:855
bool use_vector
Definition mpi.cc:748
Dictionary dict
Definition mpi.cc:959
std::map< std::size_t, index_type > data_map
Definition mpi.cc:751
std::vector< unsigned int > & owning_ranks
Definition mpi.cc:954
std::vector< std::vector< std::pair< unsigned int, std::vector< std::pair< types::global_dof_index, types::global_dof_index > > > > > requesters
Definition mpi.cc:974
std::size_t size
Definition mpi.cc:749
const MPI_Comm comm
Definition mpi.cc:928
const IndexSet & owned_indices
Definition mpi.cc:917
std::map< unsigned int, std::pair< std::vector< types::global_dof_index >, std::vector< unsigned int > > > indices_to_look_up_by_dict_rank
Definition mpi.cc:984
static const index_type invalid_index_value
Definition mpi.cc:723
const bool track_index_requesters
Definition mpi.cc:947
FlexibleIndexStorage actually_owning_ranks
Definition mpi.cc:808
types::global_dof_index locally_owned_size
Definition mpi.cc:837
unsigned int n_dict_procs_in_owned_indices
Definition mpi.cc:848
const unsigned int n_procs
Definition mpi.cc:939
types::global_dof_index dofs_per_process
Definition mpi.cc:822
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
VectorType::value_type * end(VectorType &V)
int Type_contiguous_c(MPI_Count count, MPI_Datatype oldtype, MPI_Datatype *newtype)
@ compute_point_to_point_communication_pattern
Utilities::MPI::compute_point_to_point_communication_pattern()
Definition mpi_tags.h:57
@ consensus_algorithm_payload_get_requesters
ConsensusAlgorithms::Payload::get_requesters()
Definition mpi_tags.h:81
@ dictionary_reinit
Dictionary::reinit()
Definition mpi_tags.h:78
std::vector< IndexSet > create_ascending_partitioning(const MPI_Comm comm, const types::global_dof_index locally_owned_size)
Definition mpi.cc:177
std::pair< std::vector< unsigned int >, std::map< unsigned int, IndexSet > > compute_index_owner_and_requesters(const IndexSet &owned_indices, const IndexSet &indices_to_look_up, const MPI_Comm &comm)
Definition mpi.cc:1890
T sum(const T &t, const MPI_Comm mpi_communicator)
std::unique_ptr< MPI_Datatype, void(*)(MPI_Datatype *)> create_mpi_data_type_n_bytes(const std::size_t n_bytes)
Definition mpi.cc:221
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:105
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
std::vector< T > all_gather(const MPI_Comm comm, const T &object_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:1839
std::vector< T > compute_set_union(const std::vector< T > &vec, const MPI_Comm comm)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:120
T all_reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner)
IndexSet create_evenly_distributed_partitioning(const MPI_Comm comm, const types::global_dof_index total_size)
Definition mpi.cc:206
std::vector< unsigned int > compute_point_to_point_communication_pattern(const MPI_Comm mpi_comm, const std::vector< unsigned int > &destinations)
Definition mpi.cc:263
MPI_Comm duplicate_communicator(const MPI_Comm mpi_communicator)
Definition mpi.cc:156
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
bool job_supports_mpi()
Definition mpi.cc:696
const MPI_Datatype mpi_type_id_for_type
Definition mpi.h:1634
void free_communicator(MPI_Comm mpi_communicator)
Definition mpi.cc:167
std::vector< unsigned int > mpi_processes_within_communicator(const MPI_Comm comm_large, const MPI_Comm comm_small)
Definition mpi.cc:136
unsigned int compute_n_point_to_point_communications(const MPI_Comm mpi_comm, const std::vector< unsigned int > &destinations)
Definition mpi.cc:375
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm mpi_communicator)
Definition mpi.cc:79
IndexSet create_evenly_distributed_partitioning(const unsigned int my_partition_id, const unsigned int n_partitions, const types::global_dof_index total_size)
Definition mpi.cc:55
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned int global_dof_index
Definition types.h:94
boost::signals2::signal< void()> at_mpi_finalize