deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+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
task_info.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2018 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#include <deal.II/base/conditional_ostream.h>
15#include <deal.II/base/memory_consumption.h>
16#include <deal.II/base/mpi.h>
17#include <deal.II/base/multithread_info.h>
18#include <deal.II/base/parallel.h>
19#include <deal.II/base/utilities.h>
20
21#include <deal.II/lac/dynamic_sparsity_pattern.h>
22
23#include <deal.II/matrix_free/task_info.h>
24#include <deal.II/matrix_free/util.h>
25
26
27#ifdef DEAL_II_WITH_TBB
28# include <tbb/blocked_range.h>
29# include <tbb/parallel_for.h>
30# include <tbb/task.h>
31# ifndef DEAL_II_TBB_WITH_ONEAPI
32# include <tbb/task_scheduler_init.h>
33# endif
34#endif
35
36#include <iostream>
37#include <set>
38
39//
40// TBB with oneAPI API has deprecated and removed the
41// <code>tbb::tasks</code> backend. With this it is no longer possible to
42// compile the following code that builds a directed acyclic graph (DAG) of
43// (thread parallel) tasks without a major porting effort. It turned out
44// that such a dynamic handling of dependencies and structures is not as
45// competitive as initially assumed. Consequently, this part of the matrix
46// free infrastructure has seen less attention than the rest over the last
47// years and is (presumably) not used that often.
48//
49// In case of detected oneAPI backend we simply disable threading in the
50// matrix free backend for now.
51//
52// Matthias Maier, Martin Kronbichler, 2021
53//
54
55DEAL_II_NAMESPACE_OPEN
56
57
58
59/*-------------------- Implementation of the matrix-free loop --------------*/
60namespace internal
61{
62 namespace MatrixFreeFunctions
63 {
64#if defined(DEAL_II_WITH_TBB) && !defined(DEAL_II_TBB_WITH_ONEAPI)
65
66 // This defines the TBB data structures that are needed to schedule the
67 // partition-partition variant
68
69 namespace partition
70 {
71 class ActualCellWork
72 {
73 public:
74 ActualCellWork(MFWorkerInterface **worker_pointer,
75 const unsigned int partition,
76 const TaskInfo &task_info)
77 : worker(nullptr)
78 , worker_pointer(worker_pointer)
79 , partition(partition)
80 , task_info(task_info)
81 {}
82
83 ActualCellWork(MFWorkerInterface &worker,
84 const unsigned int partition,
85 const TaskInfo &task_info)
86 : worker(&worker)
87 , worker_pointer(nullptr)
88 , partition(partition)
89 , task_info(task_info)
90 {}
91
92 void
93 operator()() const
94 {
95 MFWorkerInterface *used_worker =
96 worker != nullptr ? worker : *worker_pointer;
97 Assert(used_worker != nullptr, ExcInternalError());
98 used_worker->cell(partition);
99
100 if (task_info.face_partition_data.empty() == false)
101 {
102 used_worker->face(partition);
103 used_worker->boundary(partition);
104 }
105 }
106
107 private:
108 MFWorkerInterface *worker;
109 MFWorkerInterface **worker_pointer;
110 const unsigned int partition;
111 const TaskInfo &task_info;
112 };
113
114 class CellWork : public tbb::task
115 {
116 public:
117 CellWork(MFWorkerInterface &worker,
118 const unsigned int partition,
119 const TaskInfo &task_info,
120 const bool is_blocked)
121 : dummy(nullptr)
122 , work(worker, partition, task_info)
123 , is_blocked(is_blocked)
124 {}
125
126 tbb::task *
127 execute() override
128 {
129 work();
130
131 if (is_blocked == true)
132 tbb::empty_task::spawn(*dummy);
133 return nullptr;
134 }
135
136 tbb::empty_task *dummy;
137
138 private:
139 ActualCellWork work;
140 const bool is_blocked;
141 };
142
143
144
145 class PartitionWork : public tbb::task
146 {
147 public:
148 PartitionWork(MFWorkerInterface &function_in,
149 const unsigned int partition_in,
150 const TaskInfo &task_info_in,
151 const bool is_blocked_in = false)
152 : dummy(nullptr)
153 , function(function_in)
154 , partition(partition_in)
155 , task_info(task_info_in)
156 , is_blocked(is_blocked_in)
157 {}
158
159 tbb::task *
160 execute() override
161 {
162 tbb::empty_task *root =
163 new (tbb::task::allocate_root()) tbb::empty_task;
164 const unsigned int evens = task_info.partition_evens[partition];
165 const unsigned int odds = task_info.partition_odds[partition];
166 const unsigned int n_blocked_workers =
167 task_info.partition_n_blocked_workers[partition];
168 const unsigned int n_workers =
169 task_info.partition_n_workers[partition];
170 std::vector<CellWork *> worker(n_workers);
171 std::vector<CellWork *> blocked_worker(n_blocked_workers);
172
173 root->set_ref_count(evens + 1);
174 for (unsigned int j = 0; j < evens; ++j)
175 {
176 worker[j] = new (root->allocate_child())
177 CellWork(function,
178 task_info.partition_row_index[partition] + 2 * j,
179 task_info,
180 false);
181 if (j > 0)
182 {
183 worker[j]->set_ref_count(2);
184 blocked_worker[j - 1]->dummy =
185 new (worker[j]->allocate_child()) tbb::empty_task;
186 tbb::task::spawn(*blocked_worker[j - 1]);
187 }
188 else
189 worker[j]->set_ref_count(1);
190 if (j < evens - 1)
191 {
192 blocked_worker[j] = new (worker[j]->allocate_child())
193 CellWork(function,
194 task_info.partition_row_index[partition] + 2 * j +
195 1,
196 task_info,
197 true);
198 }
199 else
200 {
201 if (odds == evens)
202 {
203 worker[evens] = new (worker[j]->allocate_child())
204 CellWork(function,
205 task_info.partition_row_index[partition] +
206 2 * j + 1,
207 task_info,
208 false);
209 tbb::task::spawn(*worker[evens]);
210 }
211 else
212 {
213 tbb::empty_task *child =
214 new (worker[j]->allocate_child()) tbb::empty_task();
215 tbb::task::spawn(*child);
216 }
217 }
218 }
219
220 root->wait_for_all();
221 root->destroy(*root);
222 if (is_blocked == true)
223 tbb::empty_task::spawn(*dummy);
224 return nullptr;
225 }
226
227 tbb::empty_task *dummy;
228
229 private:
230 MFWorkerInterface &function;
231 const unsigned int partition;
232 const TaskInfo &task_info;
233 const bool is_blocked;
234 };
235
236 } // end of namespace partition
237
238
239
240 namespace color
241 {
242 class CellWork
243 {
244 public:
245 CellWork(MFWorkerInterface &worker_in,
246 const TaskInfo &task_info_in,
247 const unsigned int partition_in)
248 : worker(worker_in)
249 , task_info(task_info_in)
250 , partition(partition_in)
251 {}
252
253 void
254 operator()(const tbb::blocked_range<unsigned int> &r) const
255 {
256 const unsigned int start_index =
257 task_info.cell_partition_data[partition] +
258 task_info.block_size * r.begin();
259 const unsigned int end_index =
260 std::min(start_index + task_info.block_size * (r.end() - r.begin()),
261 task_info.cell_partition_data[partition + 1]);
262 worker.cell(std::make_pair(start_index, end_index));
263
264 if (task_info.face_partition_data.empty() == false)
265 {
266 AssertThrow(false, ExcNotImplemented());
267 }
268 }
269
270 private:
271 MFWorkerInterface &worker;
272 const TaskInfo &task_info;
273 const unsigned int partition;
274 };
275
276
277
278 class PartitionWork : public tbb::task
279 {
280 public:
281 PartitionWork(MFWorkerInterface &worker_in,
282 const unsigned int partition_in,
283 const TaskInfo &task_info_in,
284 const bool is_blocked_in)
285 : dummy(nullptr)
286 , worker(worker_in)
287 , partition(partition_in)
288 , task_info(task_info_in)
289 , is_blocked(is_blocked_in)
290 {}
291
292 tbb::task *
293 execute() override
294 {
295 const unsigned int n_chunks =
296 (task_info.cell_partition_data[partition + 1] -
297 task_info.cell_partition_data[partition] + task_info.block_size -
298 1) /
299 task_info.block_size;
300 parallel_for(tbb::blocked_range<unsigned int>(0, n_chunks, 1),
301 CellWork(worker, task_info, partition));
302 if (is_blocked == true)
303 tbb::empty_task::spawn(*dummy);
304 return nullptr;
305 }
306
307 tbb::empty_task *dummy;
308
309 private:
310 MFWorkerInterface &worker;
311 const unsigned int partition;
312 const TaskInfo &task_info;
313 const bool is_blocked;
314 };
315
316 } // end of namespace color
317
318
319
320 class MPICommunication : public tbb::task
321 {
322 public:
323 MPICommunication(MFWorkerInterface &worker_in, const bool do_compress)
324 : worker(worker_in)
325 , do_compress(do_compress)
326 {}
327
328 tbb::task *
329 execute() override
330 {
331 if (do_compress == false)
332 worker.vector_update_ghosts_finish();
333 else
334 worker.vector_compress_start();
335 return nullptr;
336 }
337
338 private:
339 MFWorkerInterface &worker;
340 const bool do_compress;
341 };
342
343#endif // DEAL_II_WITH_TBB
344
345
346
347 void
348 TaskInfo::loop(MFWorkerInterface &funct) const
349 {
350 // If we use thread parallelism, we do not currently support to schedule
351 // pieces of updates within the loop, so this index will collect all
352 // calls in that case and work like a single complete loop over all
353 // cells
354 if (scheme != none)
355 funct.cell_loop_pre_range(numbers::invalid_unsigned_int);
356 else
357 funct.cell_loop_pre_range(
358 partition_row_index[partition_row_index.size() - 2]);
359
360 funct.vector_update_ghosts_start();
361
362#if defined(DEAL_II_WITH_TBB) && !defined(DEAL_II_TBB_WITH_ONEAPI)
363
364 if (scheme != none)
365 {
366 funct.zero_dst_vector_range(numbers::invalid_unsigned_int);
367 if (scheme == partition_partition && evens > 0)
368 {
369 tbb::empty_task *root =
370 new (tbb::task::allocate_root()) tbb::empty_task;
371 root->set_ref_count(evens + 1);
372 std::vector<partition::PartitionWork *> worker(n_workers);
373 std::vector<partition::PartitionWork *> blocked_worker(
374 n_blocked_workers);
375 MPICommunication *worker_compr =
376 new (root->allocate_child()) MPICommunication(funct, true);
377 worker_compr->set_ref_count(1);
378 for (unsigned int j = 0; j < evens; ++j)
379 {
380 if (j > 0)
381 {
382 worker[j] = new (root->allocate_child())
383 partition::PartitionWork(funct, 2 * j, *this, false);
384 worker[j]->set_ref_count(2);
385 blocked_worker[j - 1]->dummy =
386 new (worker[j]->allocate_child()) tbb::empty_task;
387 tbb::task::spawn(*blocked_worker[j - 1]);
388 }
389 else
390 {
391 worker[j] = new (worker_compr->allocate_child())
392 partition::PartitionWork(funct, 2 * j, *this, false);
393 worker[j]->set_ref_count(2);
394 MPICommunication *worker_dist =
395 new (worker[j]->allocate_child())
396 MPICommunication(funct, false);
397 tbb::task::spawn(*worker_dist);
398 }
399 if (j < evens - 1)
400 {
401 blocked_worker[j] = new (worker[j]->allocate_child())
402 partition::PartitionWork(funct, 2 * j + 1, *this, true);
403 }
404 else
405 {
406 if (odds == evens)
407 {
408 worker[evens] = new (worker[j]->allocate_child())
409 partition::PartitionWork(funct,
410 2 * j + 1,
411 *this,
412 false);
413 tbb::task::spawn(*worker[evens]);
414 }
415 else
416 {
417 tbb::empty_task *child =
418 new (worker[j]->allocate_child()) tbb::empty_task();
419 tbb::task::spawn(*child);
420 }
421 }
422 }
423
424 root->wait_for_all();
425 root->destroy(*root);
426 }
427 else if (scheme == partition_partition)
428 {
429 // catch the case of empty partition list: we still need to call
430 // the vector communication routines to clean up and initiate
431 // things
432 funct.vector_update_ghosts_finish();
433 funct.vector_compress_start();
434 }
435 else // end of partition-partition, start of partition-color
436 {
437 // check whether there is only one partition. if not, build up the
438 // tree of partitions
439 if (odds > 0)
440 {
441 tbb::empty_task *root =
442 new (tbb::task::allocate_root()) tbb::empty_task;
443 root->set_ref_count(evens + 1);
444 const unsigned int n_blocked_workers =
445 odds - (odds + evens + 1) % 2;
446 const unsigned int n_workers =
447 cell_partition_data.size() - 1 - n_blocked_workers;
448 std::vector<color::PartitionWork *> worker(n_workers);
449 std::vector<color::PartitionWork *> blocked_worker(
450 n_blocked_workers);
451 unsigned int worker_index = 0, slice_index = 0;
452 int spawn_index_child = -2;
453 MPICommunication *worker_compr =
454 new (root->allocate_child()) MPICommunication(funct, true);
455 worker_compr->set_ref_count(1);
456 for (unsigned int part = 0;
457 part < partition_row_index.size() - 1;
458 part++)
459 {
460 if (part == 0)
461 worker[worker_index] =
462 new (worker_compr->allocate_child())
463 color::PartitionWork(funct,
464 slice_index,
465 *this,
466 false);
467 else
468 worker[worker_index] = new (root->allocate_child())
469 color::PartitionWork(funct,
470 slice_index,
471 *this,
472 false);
473 ++slice_index;
474 for (; slice_index < partition_row_index[part + 1];
475 slice_index++)
476 {
477 worker[worker_index]->set_ref_count(1);
478 ++worker_index;
479 worker[worker_index] =
480 new (worker[worker_index - 1]->allocate_child())
481 color::PartitionWork(funct,
482 slice_index,
483 *this,
484 false);
485 }
486 worker[worker_index]->set_ref_count(2);
487 if (part > 0)
488 {
489 blocked_worker[(part - 1) / 2]->dummy =
490 new (worker[worker_index]->allocate_child())
491 tbb::empty_task;
492 ++worker_index;
493 if (spawn_index_child == -1)
494 tbb::task::spawn(*blocked_worker[(part - 1) / 2]);
495 else
496 {
497 Assert(spawn_index_child >= 0,
498 ExcInternalError());
499 tbb::task::spawn(*worker[spawn_index_child]);
500 }
501 spawn_index_child = -2;
502 }
503 else
504 {
505 MPICommunication *worker_dist =
506 new (worker[worker_index]->allocate_child())
507 MPICommunication(funct, false);
508 tbb::task::spawn(*worker_dist);
509 ++worker_index;
510 }
511 part += 1;
512 if (part < partition_row_index.size() - 1)
513 {
514 if (part < partition_row_index.size() - 2)
515 {
516 blocked_worker[part / 2] =
517 new (worker[worker_index - 1]->allocate_child())
518 color::PartitionWork(funct,
519 slice_index,
520 *this,
521 true);
522 ++slice_index;
523 if (slice_index < partition_row_index[part + 1])
524 {
525 blocked_worker[part / 2]->set_ref_count(1);
526 worker[worker_index] = new (
527 blocked_worker[part / 2]->allocate_child())
528 color::PartitionWork(funct,
529 slice_index,
530 *this,
531 false);
532 ++slice_index;
533 }
534 else
535 {
536 spawn_index_child = -1;
537 continue;
538 }
539 }
540 for (; slice_index < partition_row_index[part + 1];
541 slice_index++)
542 {
543 if (slice_index > partition_row_index[part])
544 {
545 worker[worker_index]->set_ref_count(1);
546 ++worker_index;
547 }
548 worker[worker_index] =
549 new (worker[worker_index - 1]->allocate_child())
550 color::PartitionWork(funct,
551 slice_index,
552 *this,
553 false);
554 }
555 spawn_index_child = worker_index;
556 ++worker_index;
557 }
558 else
559 {
560 tbb::empty_task *final =
561 new (worker[worker_index - 1]->allocate_child())
562 tbb::empty_task;
563 tbb::task::spawn(*final);
564 spawn_index_child = worker_index - 1;
565 }
566 }
567 if (evens == odds)
568 {
569 Assert(spawn_index_child >= 0, ExcInternalError());
570 tbb::task::spawn(*worker[spawn_index_child]);
571 }
572 root->wait_for_all();
573 root->destroy(*root);
574 }
575 // case when we only have one partition: this is the usual
576 // coloring scheme, and we just schedule a parallel for loop for
577 // each color
578 else
579 {
580 Assert(evens <= 1, ExcInternalError());
581 funct.vector_update_ghosts_finish();
582
583 for (unsigned int color = 0; color < partition_row_index[1];
584 ++color)
585 {
586 tbb::empty_task *root =
587 new (tbb::task::allocate_root()) tbb::empty_task;
588 root->set_ref_count(2);
589 color::PartitionWork *worker =
590 new (root->allocate_child())
591 color::PartitionWork(funct, color, *this, false);
592 tbb::empty_task::spawn(*worker);
593 root->wait_for_all();
594 root->destroy(*root);
595 }
596
597 funct.vector_compress_start();
598 }
599 }
600 }
601 else
602#endif
603 // serial loop, go through up to three times and do the MPI transfer at
604 // the beginning/end of the second part
605 {
606 for (unsigned int part = 0; part < partition_row_index.size() - 2;
607 ++part)
608 {
609 if (part == 1)
610 funct.vector_update_ghosts_finish();
611
612 for (unsigned int i = partition_row_index[part];
613 i < partition_row_index[part + 1];
614 ++i)
615 {
616 funct.cell_loop_pre_range(i);
617 funct.zero_dst_vector_range(i);
618 AssertIndexRange(i + 1, cell_partition_data.size());
619 if (cell_partition_data[i + 1] > cell_partition_data[i])
620 {
621 funct.cell(i);
622 }
623
624 if (face_partition_data.empty() == false)
625 {
626 if (face_partition_data[i + 1] > face_partition_data[i])
627 funct.face(i);
628 if (boundary_partition_data[i + 1] >
629 boundary_partition_data[i])
630 funct.boundary(i);
631 }
632 funct.cell_loop_post_range(i);
633 }
634
635 if (part == 1)
636 funct.vector_compress_start();
637 }
638 }
639 funct.vector_compress_finish();
640
641 if (scheme != none)
642 funct.cell_loop_post_range(numbers::invalid_unsigned_int);
643 else
644 funct.cell_loop_post_range(
645 partition_row_index[partition_row_index.size() - 2]);
646 }
647
648
649
650 TaskInfo::TaskInfo()
651 {
652 clear();
653 }
654
655
656
657 void
658 TaskInfo::clear()
659 {
660 n_active_cells = 0;
661 n_ghost_cells = 0;
662 vectorization_length = 1;
663 block_size = 0;
664 n_blocks = 0;
665 scheme = none;
666 partition_row_index.clear();
667 partition_row_index.resize(2);
668 cell_partition_data.clear();
669 face_partition_data.clear();
670 boundary_partition_data.clear();
671 evens = 0;
672 odds = 0;
673 n_blocked_workers = 0;
674 n_workers = 0;
675 partition_evens.clear();
676 partition_odds.clear();
677 partition_n_blocked_workers.clear();
678 partition_n_workers.clear();
679 communicator = MPI_COMM_SELF;
680 my_pid = 0;
681 n_procs = 1;
682 }
683
684
685
686 template <typename StreamType>
687 void
688 TaskInfo::print_memory_statistics(StreamType &out,
689 const std::size_t data_length) const
690 {
691 Utilities::MPI::MinMaxAvg memory_c =
692 Utilities::MPI::min_max_avg(1e-6 * data_length, communicator);
693 if (n_procs < 2)
694 out << memory_c.min;
695 else
696 out << memory_c.min << "/" << memory_c.avg << "/" << memory_c.max;
697 out << " MB" << std::endl;
698 }
699
700
701
702 std::size_t
716
717
718
719 void
720 TaskInfo::make_boundary_cells_divisible(
721 std::vector<unsigned int> &boundary_cells)
722 {
723 // try to make the number of boundary cells divisible by the number of
724 // vectors in vectorization
725 unsigned int fillup_needed =
726 (vectorization_length - boundary_cells.size() % vectorization_length) %
727 vectorization_length;
728 if (fillup_needed > 0 && boundary_cells.size() < n_active_cells)
729 {
730 // fill additional cells into the list of boundary cells to get a
731 // balanced number. Go through the indices successively until we
732 // found enough indices
733 std::vector<unsigned int> new_boundary_cells;
734 new_boundary_cells.reserve(boundary_cells.size());
735
736 unsigned int next_free_slot = 0, bound_index = 0;
737 while (fillup_needed > 0 && bound_index < boundary_cells.size())
738 {
739 if (next_free_slot < boundary_cells[bound_index])
740 {
741 // check if there are enough cells to fill with in the
742 // current slot
743 if (next_free_slot + fillup_needed <=
744 boundary_cells[bound_index])
745 {
746 for (unsigned int j =
747 boundary_cells[bound_index] - fillup_needed;
748 j < boundary_cells[bound_index];
749 ++j)
750 new_boundary_cells.push_back(j);
751 fillup_needed = 0;
752 }
753 // ok, not enough indices, so just take them all up to the
754 // next boundary cell
755 else
756 {
757 for (unsigned int j = next_free_slot;
758 j < boundary_cells[bound_index];
759 ++j)
760 new_boundary_cells.push_back(j);
761 fillup_needed -=
762 boundary_cells[bound_index] - next_free_slot;
763 }
764 }
765 new_boundary_cells.push_back(boundary_cells[bound_index]);
766 next_free_slot = boundary_cells[bound_index] + 1;
767 ++bound_index;
768 }
769 while (fillup_needed > 0 &&
770 (new_boundary_cells.empty() ||
771 new_boundary_cells.back() < n_active_cells - 1))
772 new_boundary_cells.push_back(new_boundary_cells.back() + 1);
773 while (bound_index < boundary_cells.size())
774 new_boundary_cells.push_back(boundary_cells[bound_index++]);
775
776 boundary_cells.swap(new_boundary_cells);
777 }
778
779 // set the number of cells
780 std::sort(boundary_cells.begin(), boundary_cells.end());
781
782 // check that number of boundary cells is divisible by
783 // vectorization_length or that it contains all cells
784 Assert(boundary_cells.size() % vectorization_length == 0 ||
785 boundary_cells.size() == n_active_cells,
786 ExcInternalError());
787 }
788
789
790
791 void
792 TaskInfo::create_blocks_serial(
793 const std::vector<unsigned int> &cells_with_comm,
794 const unsigned int dofs_per_cell,
795 const bool categories_are_hp,
796 const std::vector<unsigned int> &cell_vectorization_categories,
797 const bool cell_vectorization_categories_strict,
798 const std::vector<unsigned int> &parent_relation,
799 std::vector<unsigned int> &renumbering,
800 std::vector<unsigned char> &incompletely_filled_vectorization)
801 {
802 Assert(dofs_per_cell > 0, ExcInternalError());
803 // This function is decomposed into several steps to determine a good
804 // ordering that satisfies the following constraints:
805 // a. Only cells belonging to the same category (or next higher if the
806 // cell_vectorization_categories_strict is false) can be grouped into
807 // the same SIMD batch
808 // b. hp-adaptive computations must form contiguous ranges for the same
809 // degree (category) in cell_partition_data
810 // c. We want to group the cells with the same parent in the same SIMD
811 // lane if possible
812 // d. The cell order should be similar to the initial one
813 // e. Form sets without MPI communication and those with to overlap
814 // communication with computation
815 //
816 // These constraints are satisfied by first grouping by the categories
817 // and, within the groups, to distinguish between cells with a parent
818 // and those without. All of this is set up with batches of cells (with
819 // padding if the size does not match). Then we define a vector of
820 // arrays where we define sorting criteria for the cell batches to
821 // satisfy the items b and d together, split by different parts to
822 // satisfy item e.
823
824 // Give the compiler a chance to detect that vectorization_length is a
825 // power of two, which allows it to replace integer divisions by shifts
826 const unsigned int n_lanes = indicate_power_of_two(vectorization_length);
827
828 // Step 1: find tight map of categories for not taking exceeding amounts
829 // of memory below. Sort the new categories by the numbers in the
830 // old one to ensure we respect the given rules
831 unsigned int n_categories = 1;
832 std::vector<unsigned int> tight_category_map(n_active_cells, 0);
833 if (cell_vectorization_categories.empty() == false)
834 {
835 AssertDimension(cell_vectorization_categories.size(),
836 n_active_cells + n_ghost_cells);
837
838 std::vector<unsigned int> used_categories_vector(
839 cell_vectorization_categories);
840 std::sort(used_categories_vector.begin(),
841 used_categories_vector.end());
842 used_categories_vector.erase(
843 std::unique(used_categories_vector.begin(),
844 used_categories_vector.end()),
845 used_categories_vector.end());
846
847 for (unsigned int i = 0; i < n_active_cells; ++i)
848 {
849 const unsigned int index =
850 std::lower_bound(used_categories_vector.begin(),
851 used_categories_vector.end(),
852 cell_vectorization_categories[i]) -
853 used_categories_vector.begin();
854 AssertIndexRange(index, used_categories_vector.size());
855 tight_category_map[i] = index;
856 }
857 n_categories = used_categories_vector.size();
858 }
859
860 std::vector<unsigned int> temporary_numbering;
861 temporary_numbering.reserve(n_active_cells +
862 (n_lanes - 1) * n_categories);
863 std::vector<unsigned int> category_size;
864
865 // Step 2: Sort the cells by the category
866 //
867 // If we have many categories, there is no point in trying to be clever
868 // to group things, so then we just sort the cells by the order given by
869 // the categories
870 bool do_advanced_reordering = false;
871 if (cell_vectorization_categories_strict == true &&
872 n_categories >= n_active_cells / n_lanes / 4)
873 {
874 std::vector<std::pair<unsigned int, unsigned int>> renumbered(
875 n_active_cells);
876 for (unsigned int i = 0; i < n_active_cells; ++i)
877 renumbered[i] = std::make_pair(tight_category_map[i], i);
878 std::sort(renumbered.begin(), renumbered.end());
879 for (unsigned int i = 0; i < n_active_cells;)
880 {
881 unsigned int j = 1;
882 while (i + j < n_active_cells &&
883 renumbered[i + j].first == renumbered[i].first)
884 ++j;
885 for (unsigned int k = 0; k < j; ++k)
886 temporary_numbering.push_back(renumbered[i + k].second);
887 while (temporary_numbering.size() % n_lanes != 0)
888 temporary_numbering.push_back(numbers::invalid_unsigned_int);
889 i += j;
890 }
891 }
892 // if the number of categories is the same as the number of cells and
893 // are allowed to merge categories, we simply pick the order specified
894 // by the categories.
895 else if (cell_vectorization_categories_strict == false &&
896 n_categories >= n_active_cells)
897 {
898 temporary_numbering.resize(n_active_cells);
899 for (unsigned int i = 0; i < n_active_cells; ++i)
900 temporary_numbering[tight_category_map[i]] = i;
901 }
902 // In all other cases, we perform the sorting by categories by trying to
903 // fill up the ranges in vectorization through promotion of cells to a
904 // higher category if possible, otherwise some lanes remain empty.
905 else
906 {
907 do_advanced_reordering = true;
908 std::vector<std::vector<unsigned int>> renumbering_category(
909 n_categories);
910 for (unsigned int i = 0; i < n_active_cells; ++i)
911 renumbering_category[tight_category_map[i]].push_back(i);
912
913 if (cell_vectorization_categories_strict == false && n_categories > 1)
914 for (unsigned int j = n_categories - 1; j > 0; --j)
915 {
916 unsigned int lower_index = j - 1;
917 while ((renumbering_category[j].size() % n_lanes) != 0u)
918 {
919 while (((renumbering_category[j].size() % n_lanes) != 0u) &&
920 !renumbering_category[lower_index].empty())
921 {
922 renumbering_category[j].push_back(
923 renumbering_category[lower_index].back());
924 renumbering_category[lower_index].pop_back();
925 }
926 if (lower_index == 0)
927 break;
928 else
929 --lower_index;
930 }
931 }
932
933 // Step 3: Use the parent relation to find a good grouping of cells.
934 // To do this, we first put cells of each category defined above into
935 // two bins, those which we know can be grouped together by the given
936 // parent relation and those which cannot
937 const unsigned int n_cells_per_parent =
938 std::count(parent_relation.begin(), parent_relation.end(), 0);
939 for (unsigned int j = 0; j < n_categories; ++j)
940 {
941 std::vector<std::pair<unsigned int, unsigned int>> grouped_cells;
942 std::vector<unsigned int> other_cells;
943 for (const unsigned int cell : renumbering_category[j])
944 if (parent_relation.empty() ||
945 parent_relation[cell] == numbers::invalid_unsigned_int)
946 other_cells.push_back(cell);
947 else
948 grouped_cells.emplace_back(parent_relation[cell], cell);
949
950 // Compute the number of cells per group
951 std::sort(grouped_cells.begin(), grouped_cells.end());
952 std::vector<unsigned int> n_cells_per_group;
953 unsigned int length = 0;
954 for (unsigned int i = 0; i < grouped_cells.size(); ++i, ++length)
955 if (i > 0 &&
956 grouped_cells[i].first != grouped_cells[i - 1].first)
957 {
958 n_cells_per_group.push_back(length);
959 length = 0;
960 }
961 if (length > 0)
962 n_cells_per_group.push_back(length);
963
964 // Move groups that do not have the complete size (due to
965 // categories) to the 'other_cells'. The cells with correct group
966 // size are immediately appended to the temporary cell numbering
967 auto group_it = grouped_cells.begin();
968 for (const unsigned int length : n_cells_per_group)
969 if (length < n_cells_per_parent)
970 for (unsigned int j = 0; j < length; ++j)
971 other_cells.push_back((group_it++)->second);
972 else
973 {
974 // we should not have more cells in a group than in the
975 // first check we did above
976 AssertDimension(length, n_cells_per_parent);
977 for (unsigned int j = 0; j < length; ++j)
978 temporary_numbering.push_back((group_it++)->second);
979 }
980
981 // Sort the remaining cells and append them as well
982 std::sort(other_cells.begin(), other_cells.end());
983 temporary_numbering.insert(temporary_numbering.end(),
984 other_cells.begin(),
985 other_cells.end());
986
987 while (temporary_numbering.size() % n_lanes != 0)
988 temporary_numbering.push_back(numbers::invalid_unsigned_int);
989
990 category_size.push_back(temporary_numbering.size());
991 }
992 }
993
994 while (temporary_numbering.size() % n_lanes != 0)
995 temporary_numbering.push_back(numbers::invalid_unsigned_int);
996
997 // Step 4: Identify the batches with cells marked as "comm"
998 std::vector<unsigned int> temporary_numbering_inverse(n_active_cells);
999 for (unsigned int i = 0; i < temporary_numbering.size(); ++i)
1000 if (temporary_numbering[i] != numbers::invalid_unsigned_int)
1001 temporary_numbering_inverse[temporary_numbering[i]] = i;
1002 std::vector<bool> batch_with_comm(temporary_numbering.size() / n_lanes,
1003 false);
1004 for (const unsigned int cell : cells_with_comm)
1005 batch_with_comm[temporary_numbering_inverse[cell] / n_lanes] = true;
1006
1007 // Step 5: Sort the batches of cells. We do this by two different cases:
1008 // If we have many categories according to step 2 above, where no
1009 // additional numbering is performed, we go by the numbering of the
1010 // categories in ascending order. In the other case, we sort cells by
1011 // their last cell index to get good locality, assuming that the initial
1012 // cell order is of good locality. In case we have hp-calculations with
1013 // categories, we prioritize the order given by the category.
1014 std::vector<std::array<unsigned int, 3>> batch_order;
1015 std::vector<std::array<unsigned int, 3>> batch_order_comm;
1016 for (unsigned int i = 0; i < temporary_numbering.size(); i += n_lanes)
1017 {
1018 unsigned int max_index = 0;
1019 if (do_advanced_reordering)
1020 for (unsigned int j = 0; j < n_lanes; ++j)
1021 if (temporary_numbering[i + j] < numbers::invalid_unsigned_int)
1022 max_index = std::max(temporary_numbering[i + j], max_index);
1023
1024 const unsigned int category_hp =
1025 categories_are_hp ?
1026 std::upper_bound(category_size.begin(), category_size.end(), i) -
1027 category_size.begin() :
1028 0;
1029 const std::array<unsigned int, 3> next{{category_hp, max_index, i}};
1030 if (batch_with_comm[i / n_lanes] || !do_advanced_reordering)
1031 batch_order_comm.emplace_back(next);
1032 else
1033 batch_order.emplace_back(next);
1034 }
1035
1036 std::sort(batch_order.begin(), batch_order.end());
1037 std::sort(batch_order_comm.begin(), batch_order_comm.end());
1038
1039 // Step 6: Put the cells with communication in the middle of the cell
1040 // range. For the MPI case, we need three groups to enable overlap for
1041 // communication and computation (part before comm, part with comm, part
1042 // after comm), whereas we need one for the other case. And in each
1043 // case, we allow for a slot of "ghosted" cells.
1044 std::vector<unsigned int> blocks;
1045 if (n_procs == 1)
1046 {
1047 if (batch_order.empty())
1048 std::swap(batch_order_comm, batch_order);
1049 Assert(batch_order_comm.empty(), ExcInternalError());
1050 partition_row_index.resize(3);
1051 blocks = {0, static_cast<unsigned int>(batch_order.size())};
1052 }
1053 else
1054 {
1055 partition_row_index.resize(5);
1056 const unsigned int comm_begin = batch_order.size() / 2;
1057 batch_order.insert(batch_order.begin() + comm_begin,
1058 batch_order_comm.begin(),
1059 batch_order_comm.end());
1060 const unsigned int comm_end = comm_begin + batch_order_comm.size();
1061 const unsigned int end = batch_order.size();
1062 blocks = {0, comm_begin, comm_end, end};
1063 }
1064
1065 // Step 7: sort ghost cells according to the category
1066 std::vector<std::array<unsigned int, 2>> tight_category_map_ghost;
1067
1068 if (cell_vectorization_categories.empty() == false)
1069 {
1070 tight_category_map_ghost.reserve(n_ghost_cells);
1071
1072 std::set<unsigned int> used_categories;
1073 for (unsigned int i = 0; i < n_ghost_cells; ++i)
1074 used_categories.insert(
1075 cell_vectorization_categories[i + n_active_cells]);
1076
1077 std::vector<unsigned int> used_categories_vector(
1078 used_categories.size());
1079 n_categories = 0;
1080 for (const auto &it : used_categories)
1081 used_categories_vector[n_categories++] = it;
1082
1083 std::vector<unsigned int> counters(n_categories, 0);
1084
1085 for (unsigned int i = 0; i < n_ghost_cells; ++i)
1086 {
1087 const unsigned int index =
1088 std::lower_bound(
1089 used_categories_vector.begin(),
1090 used_categories_vector.end(),
1091 cell_vectorization_categories[i + n_active_cells]) -
1092 used_categories_vector.begin();
1093 AssertIndexRange(index, used_categories_vector.size());
1094 tight_category_map_ghost.emplace_back(
1095 std::array<unsigned int, 2>{{index, i}});
1096
1097 // account for padding in the hp and strict case
1098 if (categories_are_hp || cell_vectorization_categories_strict)
1099 counters[index]++;
1100 }
1101
1102 // insert padding
1103 for (unsigned int i = 0; i < counters.size(); ++i)
1104 if (counters[i] % n_lanes != 0)
1105 for (unsigned int j = counters[i] % n_lanes; j < n_lanes; ++j)
1106 tight_category_map_ghost.emplace_back(
1107 std::array<unsigned int, 2>{
1108 {i, numbers::invalid_unsigned_int}});
1109
1110 std::sort(tight_category_map_ghost.begin(),
1111 tight_category_map_ghost.end());
1112 }
1113
1114 // Step 8: Fill in the data by batches for the locally owned cells.
1115 const unsigned int n_cell_batches = batch_order.size();
1116 const unsigned int n_ghost_batches =
1117 ((tight_category_map_ghost.empty() ? n_ghost_cells :
1118 tight_category_map_ghost.size()) +
1119 n_lanes - 1) /
1120 n_lanes;
1121 incompletely_filled_vectorization.resize(n_cell_batches +
1122 n_ghost_batches);
1123
1124 cell_partition_data.clear();
1125 cell_partition_data.resize(1, 0);
1126
1127 renumbering.clear();
1128 renumbering.resize(n_active_cells + n_ghost_cells,
1129 numbers::invalid_unsigned_int);
1130
1131 unsigned int counter = 0;
1132 for (unsigned int block = 0; block < blocks.size() - 1; ++block)
1133 {
1134 const unsigned int grain_size =
1135 std::max((2048U / dofs_per_cell) / 8 * 4, 2U);
1136 for (unsigned int k = blocks[block]; k < blocks[block + 1];
1137 k += grain_size)
1138 cell_partition_data.push_back(
1139 std::min(k + grain_size, blocks[block + 1]));
1140 partition_row_index[block + 1] = cell_partition_data.size() - 1;
1141
1142 // Set the numbering according to the reordered temporary one
1143 for (unsigned int k = blocks[block]; k < blocks[block + 1]; ++k)
1144 {
1145 const unsigned int pos = batch_order[k][2];
1146 unsigned int j = 0;
1147 for (; j < n_lanes && temporary_numbering[pos + j] !=
1148 numbers::invalid_unsigned_int;
1149 ++j)
1150 renumbering[counter++] = temporary_numbering[pos + j];
1151 if (j < n_lanes)
1152 incompletely_filled_vectorization[k] = j;
1153 }
1154 }
1155 AssertDimension(counter, n_active_cells);
1156
1157 // Step 9: Treat the ghost cells
1158 if (tight_category_map_ghost.empty())
1159 {
1160 for (unsigned int cell = 0; cell < n_ghost_cells; ++cell)
1161 renumbering[n_active_cells + cell] = n_active_cells + cell;
1162
1163 if ((n_ghost_cells % n_lanes) != 0u)
1164 incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes;
1165 }
1166 else
1167 {
1168 for (unsigned int k = 0, ptr = 0; k < n_ghost_batches;
1169 ++k, ptr += n_lanes)
1170 {
1171 unsigned int j = 0;
1172
1173 for (;
1174 j < n_lanes && (ptr + j < tight_category_map_ghost.size()) &&
1175 (tight_category_map_ghost[ptr + j][1] !=
1176 numbers::invalid_unsigned_int);
1177 ++j)
1178 renumbering[counter++] =
1179 n_active_cells + tight_category_map_ghost[ptr + j][1];
1180
1181 if (j < n_lanes)
1182 incompletely_filled_vectorization[n_cell_batches + k] = j;
1183 }
1184
1185 AssertDimension(counter, n_active_cells + n_ghost_cells);
1186 }
1187
1188 cell_partition_data.push_back(n_cell_batches + n_ghost_batches);
1189 partition_row_index.back() = cell_partition_data.size() - 1;
1190
1191 if constexpr (running_in_debug_mode())
1192 {
1193 std::vector<unsigned int> renumber_cpy(renumbering);
1194 std::sort(renumber_cpy.begin(), renumber_cpy.end());
1195 for (unsigned int i = 0; i < renumber_cpy.size(); ++i)
1196 AssertDimension(i, renumber_cpy[i]);
1197 }
1198 }
1199
1200
1201
1202 void
1203 TaskInfo::initial_setup_blocks_tasks(
1204 const std::vector<unsigned int> &boundary_cells,
1205 std::vector<unsigned int> &renumbering,
1206 std::vector<unsigned char> &incompletely_filled_vectorization)
1207 {
1208 const unsigned int n_cell_batches =
1209 (n_active_cells + vectorization_length - 1) / vectorization_length;
1210 const unsigned int n_ghost_slots =
1211 (n_ghost_cells + vectorization_length - 1) / vectorization_length;
1212 incompletely_filled_vectorization.resize(n_cell_batches + n_ghost_slots);
1213 if (n_cell_batches * vectorization_length > n_active_cells)
1214 incompletely_filled_vectorization[n_cell_batches - 1] =
1215 vectorization_length -
1216 (n_cell_batches * vectorization_length - n_active_cells);
1217 if (n_ghost_slots * vectorization_length > n_ghost_cells)
1218 incompletely_filled_vectorization[n_cell_batches + n_ghost_slots - 1] =
1219 vectorization_length -
1220 (n_ghost_slots * vectorization_length - n_ghost_cells);
1221
1222 std::vector<unsigned int> reverse_numbering(
1223 n_active_cells, numbers::invalid_unsigned_int);
1224 for (unsigned int j = 0; j < boundary_cells.size(); ++j)
1225 reverse_numbering[boundary_cells[j]] = j;
1226 unsigned int counter = boundary_cells.size();
1227 for (unsigned int j = 0; j < n_active_cells; ++j)
1228 if (reverse_numbering[j] == numbers::invalid_unsigned_int)
1229 reverse_numbering[j] = counter++;
1230
1231 AssertDimension(counter, n_active_cells);
1232 renumbering = Utilities::invert_permutation(reverse_numbering);
1233
1234 for (unsigned int j = n_active_cells; j < n_active_cells + n_ghost_cells;
1235 ++j)
1236 renumbering.push_back(j);
1237
1238 // TODO: might be able to simplify this code by not relying on the cell
1239 // partition data while computing the thread graph
1240 cell_partition_data.clear();
1241 cell_partition_data.push_back(0);
1242 if (n_procs > 1)
1243 {
1244 const unsigned int n_macro_boundary_cells =
1245 (boundary_cells.size() + vectorization_length - 1) /
1246 vectorization_length;
1247 cell_partition_data.push_back(
1248 (n_cell_batches - n_macro_boundary_cells) / 2);
1249 cell_partition_data.push_back(cell_partition_data[1] +
1250 n_macro_boundary_cells);
1251 }
1252 else
1253 AssertDimension(boundary_cells.size(), 0);
1254 cell_partition_data.push_back(n_cell_batches);
1255 cell_partition_data.push_back(cell_partition_data.back() + n_ghost_slots);
1256 partition_row_index.resize(n_procs > 1 ? 4 : 2);
1257 partition_row_index[0] = 0;
1258 partition_row_index[1] = 1;
1259 if (n_procs > 1)
1260 {
1261 partition_row_index[2] = 2;
1262 partition_row_index[3] = 3;
1263 }
1264 }
1265
1266
1267
1268 void
1269 TaskInfo::guess_block_size(const unsigned int dofs_per_cell)
1270 {
1271 // user did not say a positive number, so we have to guess
1272 if (block_size == 0)
1273 {
1274 // we would like to have enough work to do, so as first guess, try
1275 // to get 16 times as many chunks as we have threads on the system.
1276 block_size = n_active_cells / (MultithreadInfo::n_threads() * 16 *
1277 vectorization_length);
1278
1279 // if there are too few degrees of freedom per cell, need to
1280 // increase the block size
1281 const unsigned int minimum_parallel_grain_size = 200;
1282 if (dofs_per_cell * block_size < minimum_parallel_grain_size)
1283 block_size = (minimum_parallel_grain_size / dofs_per_cell + 1);
1284 if (dofs_per_cell * block_size > 10000)
1285 block_size /= 4;
1286
1287 block_size =
1288 1 << static_cast<unsigned int>(std::log2(block_size + 1));
1289 }
1290 if (block_size > n_active_cells)
1291 block_size = std::max(1U, n_active_cells);
1292 }
1293
1294
1295
1296 void
1297 TaskInfo::make_thread_graph_partition_color(
1298 DynamicSparsityPattern &connectivity_large,
1299 std::vector<unsigned int> &renumbering,
1300 std::vector<unsigned char> &irregular_cells,
1301 const bool)
1302 {
1303 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1304 if (n_cell_batches == 0)
1305 return;
1306
1307 Assert(vectorization_length > 0, ExcInternalError());
1308
1309 unsigned int partition = 0, counter = 0;
1310
1311 // Create connectivity graph for blocks based on connectivity graph for
1312 // cells.
1313 DynamicSparsityPattern connectivity(n_blocks, n_blocks);
1314 make_connectivity_cells_to_blocks(irregular_cells,
1315 connectivity_large,
1316 connectivity);
1317
1318 // Create cell-block partitioning.
1319
1320 // For each block of cells, this variable saves to which partitions the
1321 // block belongs. Initialize all to -1 to mark them as not yet assigned
1322 // a partition.
1323 std::vector<unsigned int> cell_partition(n_blocks,
1324 numbers::invalid_unsigned_int);
1325
1326 // In element j of this variable, one puts the old number of the block
1327 // that should be the jth block in the new numeration.
1328 std::vector<unsigned int> partition_list(n_blocks, 0);
1329 std::vector<unsigned int> partition_color_list(n_blocks, 0);
1330
1331 // This vector points to the start of each partition.
1332 std::vector<unsigned int> partition_size(2, 0);
1333
1334 // blocking_connectivity = true;
1335
1336 // The cluster_size in make_partitioning defines that the no. of cells
1337 // in each partition should be a multiple of cluster_size.
1338 unsigned int cluster_size = 1;
1339
1340 // Make the partitioning of the first layer of the blocks of cells.
1341 make_partitioning(connectivity,
1342 cluster_size,
1343 cell_partition,
1344 partition_list,
1345 partition_size,
1346 partition);
1347
1348 // Color the cells within each partition
1349 make_coloring_within_partitions_pre_blocked(connectivity,
1350 partition,
1351 cell_partition,
1352 partition_list,
1353 partition_size,
1354 partition_color_list);
1355
1356 partition_list = renumbering;
1357
1358 if constexpr (running_in_debug_mode())
1359 {
1360 // in debug mode, check that the partition color list is one-to-one
1361 {
1362 std::vector<unsigned int> sorted_pc_list(partition_color_list);
1363 std::sort(sorted_pc_list.begin(), sorted_pc_list.end());
1364 for (unsigned int i = 0; i < sorted_pc_list.size(); ++i)
1365 Assert(sorted_pc_list[i] == i, ExcInternalError());
1366 }
1367 }
1368
1369 // set the start list for each block and compute the renumbering of
1370 // cells
1371 std::vector<unsigned int> block_start(n_cell_batches + 1);
1372 std::vector<unsigned char> irregular(n_cell_batches);
1373
1374 unsigned int mcell_start = 0;
1375 block_start[0] = 0;
1376 for (unsigned int block = 0; block < n_blocks; ++block)
1377 {
1378 block_start[block + 1] = block_start[block];
1379 for (unsigned int mcell = mcell_start;
1380 mcell < std::min(mcell_start + block_size, n_cell_batches);
1381 ++mcell)
1382 {
1383 unsigned int n_comp = (irregular_cells[mcell] > 0) ?
1384 irregular_cells[mcell] :
1385 vectorization_length;
1386 block_start[block + 1] += n_comp;
1387 ++counter;
1388 }
1389 mcell_start += block_size;
1390 }
1391 counter = 0;
1392 unsigned int counter_macro = 0;
1393 unsigned int block_size_last =
1394 n_cell_batches - block_size * (n_blocks - 1);
1395 if (block_size_last == 0)
1396 block_size_last = block_size;
1397
1398 unsigned int tick = 0;
1399 for (unsigned int block = 0; block < n_blocks; ++block)
1400 {
1401 unsigned int present_block = partition_color_list[block];
1402 for (unsigned int cell = block_start[present_block];
1403 cell < block_start[present_block + 1];
1404 ++cell)
1405 renumbering[counter++] = partition_list[cell];
1406 unsigned int this_block_size =
1407 (present_block == n_blocks - 1) ? block_size_last : block_size;
1408
1409 // Also re-compute the content of cell_partition_data to
1410 // contain the numbers of cells, not blocks
1411 if (cell_partition_data[tick] == block)
1412 cell_partition_data[tick++] = counter_macro;
1413
1414 for (unsigned int j = 0; j < this_block_size; ++j)
1415 irregular[counter_macro++] =
1416 irregular_cells[present_block * block_size + j];
1417 }
1418 AssertDimension(tick + 1, cell_partition_data.size());
1419 cell_partition_data.back() = counter_macro;
1420
1421 irregular_cells.swap(irregular);
1422 AssertDimension(counter, n_active_cells);
1423 AssertDimension(counter_macro, n_cell_batches);
1424
1425 // check that the renumbering is one-to-one
1426 if constexpr (running_in_debug_mode())
1427 {
1428 {
1429 std::vector<unsigned int> sorted_renumbering(renumbering);
1430 std::sort(sorted_renumbering.begin(), sorted_renumbering.end());
1431 for (unsigned int i = 0; i < sorted_renumbering.size(); ++i)
1432 Assert(sorted_renumbering[i] == i, ExcInternalError());
1433 }
1434 }
1435
1436
1437 update_task_info(
1438 partition); // Actually sets too much for partition color case
1439
1440 AssertDimension(cell_partition_data.back(), n_cell_batches);
1441 }
1442
1443
1444
1445 void
1446 TaskInfo::make_thread_graph(
1447 const std::vector<unsigned int> &cell_active_fe_index,
1448 DynamicSparsityPattern &connectivity,
1449 std::vector<unsigned int> &renumbering,
1450 std::vector<unsigned char> &irregular_cells,
1451 const bool hp_bool)
1452 {
1453 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1454 if (n_cell_batches == 0)
1455 return;
1456
1457 Assert(vectorization_length > 0, ExcInternalError());
1458
1459 // if we want to block before partitioning, create connectivity graph
1460 // for blocks based on connectivity graph for cells.
1461 DynamicSparsityPattern connectivity_blocks(n_blocks, n_blocks);
1462 make_connectivity_cells_to_blocks(irregular_cells,
1463 connectivity,
1464 connectivity_blocks);
1465
1466 unsigned int n_blocks = 0;
1467 if (scheme == partition_color ||
1468 scheme == color) // blocking_connectivity == true
1469 n_blocks = this->n_blocks;
1470 else
1471 n_blocks = n_active_cells;
1472
1473 // For each block of cells, this variable saves to which partitions the
1474 // block belongs. Initialize all to -1 to mark them as not yet assigned
1475 // a partition.
1476 std::vector<unsigned int> cell_partition(n_blocks,
1477 numbers::invalid_unsigned_int);
1478
1479 // In element j of this variable, one puts the old number (but after
1480 // renumbering according to the input renumbering) of the block that
1481 // should be the jth block in the new numeration.
1482 std::vector<unsigned int> partition_list(n_blocks, 0);
1483 std::vector<unsigned int> partition_2layers_list(n_blocks, 0);
1484
1485 // This vector points to the start of each partition.
1486 std::vector<unsigned int> partition_size(2, 0);
1487
1488 unsigned int partition = 0;
1489
1490 // Within the partitions we want to be able to block for the case that
1491 // we do not block already in the connectivity. The cluster_size in
1492 // make_partitioning defines that the no. of cells in each partition
1493 // should be a multiple of cluster_size.
1494 unsigned int cluster_size = 1;
1495 if (scheme == partition_partition)
1496 cluster_size = block_size * vectorization_length;
1497
1498 // Make the partitioning of the first layer of the blocks of cells.
1499 if (scheme == partition_color || scheme == color)
1500 make_partitioning(connectivity_blocks,
1501 cluster_size,
1502 cell_partition,
1503 partition_list,
1504 partition_size,
1505 partition);
1506 else
1507 make_partitioning(connectivity,
1508 cluster_size,
1509 cell_partition,
1510 partition_list,
1511 partition_size,
1512 partition);
1513
1514 // Partition or color second layer
1515 if (scheme == partition_partition)
1516
1517 {
1518 // Partition within partitions.
1519 make_partitioning_within_partitions_post_blocked(
1520 connectivity,
1521 cell_active_fe_index,
1522 partition,
1523 cluster_size,
1524 hp_bool,
1525 cell_partition,
1526 partition_list,
1527 partition_size,
1528 partition_2layers_list,
1529 irregular_cells);
1530 }
1531 else if (scheme == partition_color || scheme == color)
1532 {
1533 make_coloring_within_partitions_pre_blocked(connectivity_blocks,
1534 partition,
1535 cell_partition,
1536 partition_list,
1537 partition_size,
1538 partition_2layers_list);
1539 }
1540
1541 // in debug mode, check that the partition_2layers_list is one-to-one
1542 if constexpr (running_in_debug_mode())
1543 {
1544 {
1545 std::vector<unsigned int> sorted_pc_list(partition_2layers_list);
1546 std::sort(sorted_pc_list.begin(), sorted_pc_list.end());
1547 for (unsigned int i = 0; i < sorted_pc_list.size(); ++i)
1548 Assert(sorted_pc_list[i] == i, ExcInternalError());
1549 }
1550 }
1551
1552 // Set the new renumbering
1553 std::vector<unsigned int> renumbering_in(n_active_cells, 0);
1554 renumbering_in.swap(renumbering);
1555 if (scheme == partition_partition) // blocking_connectivity == false
1556 {
1557 // This is the simple case. The renumbering is just a combination of
1558 // the renumbering that we were given as an input and the
1559 // renumbering of partition/coloring given in partition_2layers_list
1560 for (unsigned int j = 0; j < renumbering.size(); ++j)
1561 renumbering[j] = renumbering_in[partition_2layers_list[j]];
1562 // Account for the ghost cells, finally.
1563 for (unsigned int i = 0; i < n_ghost_cells; ++i)
1564 renumbering.push_back(i + n_active_cells);
1565 }
1566 else
1567 {
1568 // set the start list for each block and compute the renumbering of
1569 // cells
1570 std::vector<unsigned int> block_start(n_cell_batches + 1);
1571 std::vector<unsigned char> irregular(n_cell_batches);
1572
1573 unsigned int counter = 0;
1574 unsigned int mcell_start = 0;
1575 block_start[0] = 0;
1576 for (unsigned int block = 0; block < n_blocks; ++block)
1577 {
1578 block_start[block + 1] = block_start[block];
1579 for (unsigned int mcell = mcell_start;
1580 mcell < std::min(mcell_start + block_size, n_cell_batches);
1581 ++mcell)
1582 {
1583 unsigned int n_comp = (irregular_cells[mcell] > 0) ?
1584 irregular_cells[mcell] :
1585 vectorization_length;
1586 block_start[block + 1] += n_comp;
1587 ++counter;
1588 }
1589 mcell_start += block_size;
1590 }
1591 counter = 0;
1592 unsigned int counter_macro = 0;
1593 unsigned int block_size_last =
1594 n_cell_batches - block_size * (n_blocks - 1);
1595 if (block_size_last == 0)
1596 block_size_last = block_size;
1597
1598 unsigned int tick = 0;
1599 for (unsigned int block = 0; block < n_blocks; ++block)
1600 {
1601 unsigned int present_block = partition_2layers_list[block];
1602 for (unsigned int cell = block_start[present_block];
1603 cell < block_start[present_block + 1];
1604 ++cell)
1605 renumbering[counter++] = renumbering_in[cell];
1606 unsigned int this_block_size =
1607 (present_block == n_blocks - 1) ? block_size_last : block_size;
1608
1609 // Also re-compute the content of cell_partition_data to
1610 // contain the numbers of cells, not blocks
1611 if (cell_partition_data[tick] == block)
1612 cell_partition_data[tick++] = counter_macro;
1613
1614 for (unsigned int j = 0; j < this_block_size; ++j)
1615 irregular[counter_macro++] =
1616 irregular_cells[present_block * block_size + j];
1617 }
1618 AssertDimension(tick + 1, cell_partition_data.size());
1619 cell_partition_data.back() = counter_macro;
1620
1621 irregular_cells.swap(irregular);
1622 AssertDimension(counter, n_active_cells);
1623 AssertDimension(counter_macro, n_cell_batches);
1624 // check that the renumbering is one-to-one
1625 if constexpr (running_in_debug_mode())
1626 {
1627 {
1628 std::vector<unsigned int> sorted_renumbering(renumbering);
1629 std::sort(sorted_renumbering.begin(), sorted_renumbering.end());
1630 for (unsigned int i = 0; i < sorted_renumbering.size(); ++i)
1631 Assert(sorted_renumbering[i] == i, ExcInternalError());
1632 }
1633 }
1634 }
1635
1636 // Update the task_info with the more information for the thread graph.
1637 update_task_info(partition);
1638 }
1639
1640
1641
1642 void
1643 TaskInfo::make_thread_graph_partition_partition(
1644 const std::vector<unsigned int> &cell_active_fe_index,
1645 DynamicSparsityPattern &connectivity,
1646 std::vector<unsigned int> &renumbering,
1647 std::vector<unsigned char> &irregular_cells,
1648 const bool hp_bool)
1649 {
1650 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1651 if (n_cell_batches == 0)
1652 return;
1653
1654 const unsigned int cluster_size = block_size * vectorization_length;
1655
1656 // Create cell-block partitioning.
1657
1658 // For each block of cells, this variable saves to which partitions the
1659 // block belongs. Initialize all to n_cell_batches to mark them as not
1660 // yet assigned a partition.
1661 std::vector<unsigned int> cell_partition(n_active_cells,
1662 numbers::invalid_unsigned_int);
1663
1664
1665 // In element j of this variable, one puts the old number of the block
1666 // that should be the jth block in the new numeration.
1667 std::vector<unsigned int> partition_list(n_active_cells, 0);
1668 std::vector<unsigned int> partition_partition_list(n_active_cells, 0);
1669
1670 // This vector points to the start of each partition.
1671 std::vector<unsigned int> partition_size(2, 0);
1672
1673 unsigned int partition = 0;
1674 // Here, we do not block inside the connectivity graph
1675 // blocking_connectivity = false;
1676
1677 // Make the partitioning of the first layer of the blocks of cells.
1678 make_partitioning(connectivity,
1679 cluster_size,
1680 cell_partition,
1681 partition_list,
1682 partition_size,
1683 partition);
1684
1685 // Partition within partitions.
1686 make_partitioning_within_partitions_post_blocked(connectivity,
1687 cell_active_fe_index,
1688 partition,
1689 cluster_size,
1690 hp_bool,
1691 cell_partition,
1692 partition_list,
1693 partition_size,
1694 partition_partition_list,
1695 irregular_cells);
1696
1697 partition_list.swap(renumbering);
1698
1699 for (unsigned int j = 0; j < renumbering.size(); ++j)
1700 renumbering[j] = partition_list[partition_partition_list[j]];
1701
1702 for (unsigned int i = 0; i < n_ghost_cells; ++i)
1703 renumbering.push_back(i + n_active_cells);
1704
1705 update_task_info(partition);
1706 }
1707
1708
1709
1710 void
1711 TaskInfo::make_connectivity_cells_to_blocks(
1712 const std::vector<unsigned char> &irregular_cells,
1713 const DynamicSparsityPattern &connectivity_cells,
1714 DynamicSparsityPattern &connectivity_blocks) const
1715 {
1716 std::vector<std::vector<unsigned int>> cell_blocks(n_blocks);
1717 std::vector<unsigned int> touched_cells(n_active_cells);
1718 unsigned int cell = 0;
1719 for (unsigned int i = 0, mcell = 0; i < n_blocks; ++i)
1720 {
1721 for (unsigned int c = 0;
1722 c < block_size && mcell < *(cell_partition_data.end() - 2);
1723 ++c, ++mcell)
1724 {
1725 unsigned int ncomp = (irregular_cells[mcell] > 0) ?
1726 irregular_cells[mcell] :
1727 vectorization_length;
1728 for (unsigned int c = 0; c < ncomp; ++c, ++cell)
1729 {
1730 cell_blocks[i].push_back(cell);
1731 touched_cells[cell] = i;
1732 }
1733 }
1734 }
1735 AssertDimension(cell, n_active_cells);
1736 for (unsigned int i = 0; i < cell_blocks.size(); ++i)
1737 for (unsigned int col = 0; col < cell_blocks[i].size(); ++col)
1738 {
1739 for (DynamicSparsityPattern::iterator it =
1740 connectivity_cells.begin(cell_blocks[i][col]);
1741 it != connectivity_cells.end(cell_blocks[i][col]);
1742 ++it)
1743 {
1744 if (touched_cells[it->column()] != i)
1745 connectivity_blocks.add(i, touched_cells[it->column()]);
1746 }
1747 }
1748 }
1749
1750
1751
1752 // Function to create partitioning on the second layer within each
1753 // partition. Version without preblocking.
1754 void
1755 TaskInfo::make_partitioning_within_partitions_post_blocked(
1756 const DynamicSparsityPattern &connectivity,
1757 const std::vector<unsigned int> &cell_active_fe_index,
1758 const unsigned int partition,
1759 const unsigned int cluster_size,
1760 const bool hp_bool,
1761 const std::vector<unsigned int> &cell_partition,
1762 const std::vector<unsigned int> &partition_list,
1763 const std::vector<unsigned int> &partition_size,
1764 std::vector<unsigned int> &partition_partition_list,
1765 std::vector<unsigned char> &irregular_cells)
1766 {
1767 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1768 const unsigned int n_ghost_slots =
1769 *(cell_partition_data.end() - 1) - n_cell_batches;
1770
1771 // List of cells in previous partition
1772 std::vector<unsigned int> neighbor_list;
1773 // List of cells in current partition for use as neighbors in next
1774 // partition
1775 std::vector<unsigned int> neighbor_neighbor_list;
1776
1777 std::vector<unsigned int> renumbering(n_active_cells);
1778
1779 irregular_cells.back() = 0;
1780 irregular_cells.resize(n_active_cells + n_ghost_slots);
1781
1782 unsigned int max_fe_index = 0;
1783 for (const unsigned int fe_index : cell_active_fe_index)
1784 max_fe_index = std::max(fe_index, max_fe_index);
1785
1786 Assert(!hp_bool || cell_active_fe_index.size() == n_active_cells,
1787 ExcInternalError());
1788
1789 {
1790 unsigned int n_cell_batches_before = 0;
1791 // Create partitioning within partitions.
1792
1793 // For each block of cells, this variable saves to which partitions
1794 // the block belongs. Initialize all to n_cell_batches to mark them as
1795 // not yet assigned a partition.
1796 std::vector<unsigned int> cell_partition_l2(
1797 n_active_cells, numbers::invalid_unsigned_int);
1798 partition_row_index.clear();
1799 partition_row_index.resize(partition + 1, 0);
1800 cell_partition_data.resize(1, 0);
1801
1802 unsigned int counter = 0;
1803 unsigned int missing_macros;
1804 for (unsigned int part = 0; part < partition; ++part)
1805 {
1806 neighbor_neighbor_list.resize(0);
1807 neighbor_list.resize(0);
1808 bool work = true;
1809 unsigned int partition_l2 = 0;
1810 unsigned int start_up = partition_size[part];
1811 unsigned int partition_counter = 0;
1812 while (work)
1813 {
1814 if (neighbor_list.empty())
1815 {
1816 work = false;
1817 partition_counter = 0;
1818 for (unsigned int j = start_up;
1819 j < partition_size[part + 1];
1820 ++j)
1821 if (cell_partition[partition_list[j]] == part &&
1822 cell_partition_l2[partition_list[j]] ==
1823 numbers::invalid_unsigned_int)
1824 {
1825 start_up = j;
1826 work = true;
1827 partition_counter = 1;
1828 // To start up, set the start_up cell to partition
1829 // and list all its neighbors.
1830 AssertIndexRange(start_up, partition_size[part + 1]);
1831 cell_partition_l2[partition_list[start_up]] =
1832 partition_l2;
1833 neighbor_neighbor_list.push_back(
1834 partition_list[start_up]);
1835 partition_partition_list[counter++] =
1836 partition_list[start_up];
1837 ++start_up;
1838 break;
1839 }
1840 }
1841 else
1842 {
1843 partition_counter = 0;
1844 for (const unsigned int neighbor : neighbor_list)
1845 {
1846 Assert(cell_partition[neighbor] == part,
1847 ExcInternalError());
1848 Assert(cell_partition_l2[neighbor] == partition_l2 - 1,
1849 ExcInternalError());
1850 auto neighbor_it = connectivity.begin(neighbor);
1851 const auto end_it = connectivity.end(neighbor);
1852 for (; neighbor_it != end_it; ++neighbor_it)
1853 {
1854 if (cell_partition[neighbor_it->column()] == part &&
1855 cell_partition_l2[neighbor_it->column()] ==
1856 numbers::invalid_unsigned_int)
1857 {
1858 cell_partition_l2[neighbor_it->column()] =
1859 partition_l2;
1860 neighbor_neighbor_list.push_back(
1861 neighbor_it->column());
1862 partition_partition_list[counter++] =
1863 neighbor_it->column();
1864 ++partition_counter;
1865 }
1866 }
1867 }
1868 }
1869 if (partition_counter > 0)
1870 {
1871 int index_before = neighbor_neighbor_list.size(),
1872 index = index_before;
1873 {
1874 // put the cells into separate lists for each FE index
1875 // within one partition-partition
1876 missing_macros = 0;
1877 std::vector<unsigned int> remaining_per_cell_batch(
1878 max_fe_index + 1);
1879 std::vector<std::vector<unsigned int>>
1880 renumbering_fe_index;
1881 unsigned int cell;
1882 bool filled = true;
1883 if (hp_bool == true)
1884 {
1885 renumbering_fe_index.resize(max_fe_index + 1);
1886 for (cell = counter - partition_counter;
1887 cell < counter;
1888 ++cell)
1889 {
1890 renumbering_fe_index
1891 [cell_active_fe_index.empty() ?
1892 0 :
1893 cell_active_fe_index
1894 [partition_partition_list[cell]]]
1895 .push_back(partition_partition_list[cell]);
1896 }
1897 // check how many more cells are needed in the lists
1898 for (unsigned int j = 0; j < max_fe_index + 1; ++j)
1899 {
1900 remaining_per_cell_batch[j] =
1901 renumbering_fe_index[j].size() %
1902 vectorization_length;
1903 if (remaining_per_cell_batch[j] != 0)
1904 filled = false;
1905 missing_macros +=
1906 ((renumbering_fe_index[j].size() +
1907 vectorization_length - 1) /
1908 vectorization_length);
1909 }
1910 }
1911 else
1912 {
1913 remaining_per_cell_batch.resize(1);
1914 remaining_per_cell_batch[0] =
1915 partition_counter % vectorization_length;
1916 missing_macros =
1917 partition_counter / vectorization_length;
1918 if (remaining_per_cell_batch[0] != 0)
1919 {
1920 filled = false;
1921 ++missing_macros;
1922 }
1923 }
1924 missing_macros =
1925 cluster_size - (missing_macros % cluster_size);
1926
1927 // now we realized that there are some cells missing.
1928 while (missing_macros > 0 || filled == false)
1929 {
1930 if (index == 0)
1931 {
1932 index = neighbor_neighbor_list.size();
1933 if (index == index_before)
1934 {
1935 if (missing_macros != 0)
1936 {
1937 neighbor_neighbor_list.resize(0);
1938 }
1939 start_up--;
1940 break; // not connected - start again
1941 }
1942 index_before = index;
1943 }
1944 index--;
1945 unsigned int additional =
1946 neighbor_neighbor_list[index];
1947
1948 // go through the neighbors of the last cell in the
1949 // current partition and check if we find some to
1950 // fill up with.
1951 DynamicSparsityPattern::iterator neighbor =
1952 connectivity.begin(
1953 additional),
1954 end =
1955 connectivity.end(
1956 additional);
1957 for (; neighbor != end; ++neighbor)
1958 {
1959 if (cell_partition[neighbor->column()] == part &&
1960 cell_partition_l2[neighbor->column()] ==
1961 numbers::invalid_unsigned_int)
1962 {
1963 unsigned int this_index = 0;
1964 if (hp_bool == true)
1965 this_index =
1966 cell_active_fe_index.empty() ?
1967 0 :
1968 cell_active_fe_index[neighbor
1969 ->column()];
1970
1971 // Only add this cell if we need more macro
1972 // cells in the current block or if there is
1973 // a macro cell with the FE index that is
1974 // not yet fully populated
1975 if (missing_macros > 0 ||
1976 remaining_per_cell_batch[this_index] > 0)
1977 {
1978 cell_partition_l2[neighbor->column()] =
1979 partition_l2;
1980 neighbor_neighbor_list.push_back(
1981 neighbor->column());
1982 if (hp_bool == true)
1983 renumbering_fe_index[this_index]
1984 .push_back(neighbor->column());
1985 partition_partition_list[counter] =
1986 neighbor->column();
1987 ++counter;
1988 ++partition_counter;
1989 if (remaining_per_cell_batch
1990 [this_index] == 0 &&
1991 missing_macros > 0)
1992 missing_macros--;
1993 remaining_per_cell_batch[this_index]++;
1994 if (remaining_per_cell_batch
1995 [this_index] ==
1996 vectorization_length)
1997 {
1998 remaining_per_cell_batch[this_index] =
1999 0;
2000 }
2001 if (missing_macros == 0)
2002 {
2003 filled = true;
2004 for (unsigned int fe_ind = 0;
2005 fe_ind < max_fe_index + 1;
2006 ++fe_ind)
2007 if (remaining_per_cell_batch
2008 [fe_ind] != 0)
2009 filled = false;
2010 }
2011 if (filled == true)
2012 break;
2013 }
2014 }
2015 }
2016 }
2017 if (hp_bool == true)
2018 {
2019 // set the renumbering according to their active FE
2020 // index within one partition-partition which was
2021 // implicitly assumed above
2022 cell = counter - partition_counter;
2023 for (unsigned int j = 0; j < max_fe_index + 1; ++j)
2024 {
2025 for (const unsigned int jj :
2026 renumbering_fe_index[j])
2027 renumbering[cell++] = jj;
2028 if (renumbering_fe_index[j].size() %
2029 vectorization_length !=
2030 0)
2031 irregular_cells[renumbering_fe_index[j].size() /
2032 vectorization_length +
2033 n_cell_batches_before] =
2034 renumbering_fe_index[j].size() %
2035 vectorization_length;
2036 n_cell_batches_before +=
2037 (renumbering_fe_index[j].size() +
2038 vectorization_length - 1) /
2039 vectorization_length;
2040 renumbering_fe_index[j].resize(0);
2041 }
2042 }
2043 else
2044 {
2045 n_cell_batches_before +=
2046 partition_counter / vectorization_length;
2047 if (partition_counter % vectorization_length != 0)
2048 {
2049 irregular_cells[n_cell_batches_before] =
2050 partition_counter % vectorization_length;
2051 ++n_cell_batches_before;
2052 }
2053 }
2054 }
2055 cell_partition_data.push_back(n_cell_batches_before);
2056 partition_l2++;
2057 }
2058 neighbor_list = neighbor_neighbor_list;
2059 neighbor_neighbor_list.resize(0);
2060 }
2061 partition_row_index[part + 1] =
2062 partition_row_index[part] + partition_l2;
2063 }
2064 }
2065 if (hp_bool == true)
2066 {
2067 partition_partition_list.swap(renumbering);
2068 }
2069 }
2070
2071
2072
2073 // Function to create coloring on the second layer within each partition.
2074 // Version assumes preblocking.
2075 void
2076 TaskInfo::make_coloring_within_partitions_pre_blocked(
2077 const DynamicSparsityPattern &connectivity,
2078 const unsigned int partition,
2079 const std::vector<unsigned int> &cell_partition,
2080 const std::vector<unsigned int> &partition_list,
2081 const std::vector<unsigned int> &partition_size,
2082 std::vector<unsigned int> &partition_color_list)
2083 {
2084 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
2085 std::vector<unsigned int> cell_color(n_blocks, n_cell_batches);
2086 std::vector<bool> color_finder;
2087
2088 partition_row_index.resize(partition + 1);
2089 cell_partition_data.clear();
2090 unsigned int color_counter = 0, index_counter = 0;
2091 for (unsigned int part = 0; part < partition; ++part)
2092 {
2093 partition_row_index[part] = index_counter;
2094 unsigned int max_color = 0;
2095 for (unsigned int k = partition_size[part];
2096 k < partition_size[part + 1];
2097 k++)
2098 {
2099 unsigned int cell = partition_list[k];
2100 unsigned int n_neighbors = connectivity.row_length(cell);
2101
2102 // In the worst case, each neighbor has a different color. So we
2103 // find at least one available color between 0 and n_neighbors.
2104 color_finder.resize(n_neighbors + 1);
2105 for (unsigned int j = 0; j <= n_neighbors; ++j)
2106 color_finder[j] = true;
2107 DynamicSparsityPattern::iterator neighbor =
2108 connectivity.begin(cell),
2109 end = connectivity.end(cell);
2110 for (; neighbor != end; ++neighbor)
2111 {
2112 // Mark the color that a neighbor within the partition has
2113 // as taken
2114 if (cell_partition[neighbor->column()] == part &&
2115 cell_color[neighbor->column()] <= n_neighbors)
2116 color_finder[cell_color[neighbor->column()]] = false;
2117 }
2118 // Choose the smallest color that is not taken for the block
2119 cell_color[cell] = 0;
2120 while (color_finder[cell_color[cell]] == false)
2121 cell_color[cell]++;
2122 if (cell_color[cell] > max_color)
2123 max_color = cell_color[cell];
2124 }
2125 // Reorder within partition: First, all blocks that belong the 0 and
2126 // then so on until those with color max (Note that the smaller the
2127 // number the larger the partition)
2128 for (unsigned int color = 0; color <= max_color; ++color)
2129 {
2130 cell_partition_data.push_back(color_counter);
2131 ++index_counter;
2132 for (unsigned int k = partition_size[part];
2133 k < partition_size[part + 1];
2134 k++)
2135 {
2136 unsigned int cell = partition_list[k];
2137 if (cell_color[cell] == color)
2138 {
2139 partition_color_list[color_counter++] = cell;
2140 }
2141 }
2142 }
2143 }
2144 cell_partition_data.push_back(n_blocks);
2145 partition_row_index[partition] = index_counter;
2146 AssertDimension(color_counter, n_blocks);
2147 }
2148
2149
2150 // Function to create partitioning on the first layer.
2151 void
2152 TaskInfo::make_partitioning(const DynamicSparsityPattern &connectivity,
2153 const unsigned int cluster_size,
2154 std::vector<unsigned int> &cell_partition,
2155 std::vector<unsigned int> &partition_list,
2156 std::vector<unsigned int> &partition_size,
2157 unsigned int &partition) const
2158
2159 {
2160 // For each block of cells, this variable saves to which partitions the
2161 // block belongs. Initialize all to n_cell_batches to mark them as not
2162 // yet assigned a partition.
2163 // std::vector<unsigned int> cell_partition (n_active_cells,
2164 // numbers::invalid_unsigned_int);
2165 // List of cells in previous partition
2166 std::vector<unsigned int> neighbor_list;
2167 // List of cells in current partition for use as neighbors in next
2168 // partition
2169 std::vector<unsigned int> neighbor_neighbor_list;
2170
2171 // In element j of this variable, one puts the old number of the block
2172 // that should be the jth block in the new numeration.
2173 // std::vector<unsigned int> partition_list(n_active_cells,0);
2174
2175 // This vector points to the start of each partition.
2176 // std::vector<unsigned int> partition_size(2,0);
2177
2178 partition = 0;
2179 unsigned int counter = 0;
2180 unsigned int start_nonboundary =
2181 cell_partition_data.size() == 5 ?
2182 vectorization_length *
2183 (cell_partition_data[2] - cell_partition_data[1]) :
2184 0;
2185
2186 const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
2187 if (n_cell_batches == 0)
2188 return;
2189 if (scheme == color)
2190 start_nonboundary = n_cell_batches;
2191 if (scheme == partition_color ||
2192 scheme == color) // blocking_connectivity == true
2193 start_nonboundary = ((start_nonboundary + block_size - 1) / block_size);
2194 unsigned int n_blocks;
2195 if (scheme == partition_color ||
2196 scheme == color) // blocking_connectivity == true
2197 n_blocks = this->n_blocks;
2198 else
2199 n_blocks = n_active_cells;
2200
2201 if (start_nonboundary > n_blocks)
2202 start_nonboundary = n_blocks;
2203
2204
2205 unsigned int start_up = 0;
2206 bool work = true;
2207 unsigned int remainder = cluster_size;
2208
2209 // this performs a classical breath-first search in the connectivity
2210 // graph of the cells under the restriction that the size of the
2211 // partitions should be a multiple of the given block size
2212 while (work)
2213 {
2214 // put the cells with neighbors on remote MPI processes up front
2215 if (start_nonboundary > 0)
2216 {
2217 for (unsigned int cell = 0; cell < start_nonboundary; ++cell)
2218 {
2219 const unsigned int cell_nn = cell;
2220 cell_partition[cell_nn] = partition;
2221 neighbor_list.push_back(cell_nn);
2222 partition_list[counter++] = cell_nn;
2223 partition_size.back()++;
2224 }
2225 start_nonboundary = 0;
2226 remainder -= (start_nonboundary % cluster_size);
2227 if (remainder == cluster_size)
2228 remainder = 0;
2229 }
2230 else
2231 {
2232 // To start up, set the start_up cell to partition and list all
2233 // its neighbors.
2234 cell_partition[start_up] = partition;
2235 neighbor_list.push_back(start_up);
2236 partition_list[counter++] = start_up;
2237 partition_size.back()++;
2238 ++start_up;
2239 remainder--;
2240 if (remainder == cluster_size)
2241 remainder = 0;
2242 }
2243 int index_before = neighbor_list.size(), index = index_before,
2244 index_stop = 0;
2245 while (remainder > 0)
2246 {
2247 if (index == index_stop)
2248 {
2249 index = neighbor_list.size();
2250 if (index == index_before)
2251 {
2252 neighbor_list.resize(0);
2253 goto not_connect;
2254 }
2255 index_stop = index_before;
2256 index_before = index;
2257 }
2258 index--;
2259 unsigned int additional = neighbor_list[index];
2260 DynamicSparsityPattern::iterator neighbor =
2261 connectivity.begin(additional),
2262 end =
2263 connectivity.end(additional);
2264 for (; neighbor != end; ++neighbor)
2265 {
2266 if (cell_partition[neighbor->column()] ==
2267 numbers::invalid_unsigned_int)
2268 {
2269 partition_size.back()++;
2270 cell_partition[neighbor->column()] = partition;
2271 neighbor_list.push_back(neighbor->column());
2272 partition_list[counter++] = neighbor->column();
2273 remainder--;
2274 if (remainder == 0)
2275 break;
2276 }
2277 }
2278 }
2279
2280 while (neighbor_list.size() > 0)
2281 {
2282 ++partition;
2283
2284 // counter for number of cells so far in current partition
2285 unsigned int partition_counter = 0;
2286
2287 // Mark the start of the new partition
2288 partition_size.push_back(partition_size.back());
2289
2290 // Loop through the list of cells in previous partition and put
2291 // all their neighbors in current partition
2292 for (const unsigned int cell : neighbor_list)
2293 {
2294 Assert(cell_partition[cell] == partition - 1,
2295 ExcInternalError());
2296 auto neighbor = connectivity.begin(cell);
2297 const auto end = connectivity.end(cell);
2298 for (; neighbor != end; ++neighbor)
2299 {
2300 if (cell_partition[neighbor->column()] ==
2301 numbers::invalid_unsigned_int)
2302 {
2303 partition_size.back()++;
2304 cell_partition[neighbor->column()] = partition;
2305
2306 // collect the cells of the current partition for
2307 // use as neighbors in next partition
2308 neighbor_neighbor_list.push_back(neighbor->column());
2309 partition_list[counter++] = neighbor->column();
2310 ++partition_counter;
2311 }
2312 }
2313 }
2314 remainder = cluster_size - (partition_counter % cluster_size);
2315 if (remainder == cluster_size)
2316 remainder = 0;
2317 int index_stop = 0;
2318 int index_before = neighbor_neighbor_list.size(),
2319 index = index_before;
2320 while (remainder > 0)
2321 {
2322 if (index == index_stop)
2323 {
2324 index = neighbor_neighbor_list.size();
2325 if (index == index_before)
2326 {
2327 neighbor_neighbor_list.resize(0);
2328 break;
2329 }
2330 index_stop = index_before;
2331 index_before = index;
2332 }
2333 index--;
2334 unsigned int additional = neighbor_neighbor_list[index];
2335 DynamicSparsityPattern::iterator neighbor =
2336 connectivity.begin(
2337 additional),
2338 end = connectivity.end(
2339 additional);
2340 for (; neighbor != end; ++neighbor)
2341 {
2342 if (cell_partition[neighbor->column()] ==
2343 numbers::invalid_unsigned_int)
2344 {
2345 partition_size.back()++;
2346 cell_partition[neighbor->column()] = partition;
2347 neighbor_neighbor_list.push_back(neighbor->column());
2348 partition_list[counter++] = neighbor->column();
2349 remainder--;
2350 if (remainder == 0)
2351 break;
2352 }
2353 }
2354 }
2355
2356 neighbor_list = neighbor_neighbor_list;
2357 neighbor_neighbor_list.resize(0);
2358 }
2359 not_connect:
2360 // One has to check if the graph is not connected so we have to find
2361 // another partition.
2362 work = false;
2363 for (unsigned int j = start_up; j < n_blocks; ++j)
2364 if (cell_partition[j] == numbers::invalid_unsigned_int)
2365 {
2366 start_up = j;
2367 work = true;
2368 if (remainder == 0)
2369 remainder = cluster_size;
2370 break;
2371 }
2372 }
2373 if (remainder != 0)
2374 ++partition;
2375
2376 AssertDimension(partition_size[partition], n_blocks);
2377 }
2378
2379
2380 void
2381 TaskInfo::update_task_info(const unsigned int partition)
2382 {
2383 evens = (partition + 1) / 2;
2384 odds = partition / 2;
2385 n_blocked_workers = odds - (odds + evens + 1) % 2;
2386 n_workers = evens + odds - n_blocked_workers;
2387 // From here only used for partition partition option.
2388 partition_evens.resize(partition);
2389 partition_odds.resize(partition);
2390 partition_n_blocked_workers.resize(partition);
2391 partition_n_workers.resize(partition);
2392 for (unsigned int part = 0; part < partition; ++part)
2393 {
2394 partition_evens[part] =
2395 (partition_row_index[part + 1] - partition_row_index[part] + 1) / 2;
2396 partition_odds[part] =
2397 (partition_row_index[part + 1] - partition_row_index[part]) / 2;
2398 partition_n_blocked_workers[part] =
2399 partition_odds[part] -
2400 (partition_odds[part] + partition_evens[part] + 1) % 2;
2401 partition_n_workers[part] = partition_evens[part] +
2402 partition_odds[part] -
2403 partition_n_blocked_workers[part];
2404 }
2405 }
2406 } // namespace MatrixFreeFunctions
2407} // namespace internal
2408
2409
2410
2411// explicit instantiations of template functions
2412template void
2413internal::MatrixFreeFunctions::TaskInfo::print_memory_statistics<std::ostream>(
2414 std::ostream &,
2415 const std::size_t) const;
2416template void
2417internal::MatrixFreeFunctions::TaskInfo::print_memory_statistics<
2418 ConditionalOStream>(ConditionalOStream &, const std::size_t) const;
2419
2420
2421DEAL_II_NAMESPACE_CLOSE
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10398
operator()
*  *  Point< dim > operator()(const Point< dim > &p) const * 
Definition ALE_Finite_Strain_Plasticity.h:7901
push_back
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
DynamicSparsityPattern::row_length
size_type row_length(const size_type row) const
Definition dynamic_sparsity_pattern.h:1074
DynamicSparsityPattern::add
void add(const size_type i, const size_type j)
Definition dynamic_sparsity_pattern.h:1031
MultithreadInfo::n_threads
static unsigned int n_threads()
Definition multithread_info.cc:127
conditional_ostream.h
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
running_in_debug_mode
constexpr bool running_in_debug_mode()
Definition config.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
dynamic_sparsity_pattern.h
first
Point< 2 > first
Definition grid_out.cc:4639
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
memory_consumption.h
size
std::size_t size
Definition mpi.cc:733
n_procs
const unsigned int n_procs
Definition mpi.cc:923
multithread_info.h
MemoryConsumption::memory_consumption
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Definition memory_consumption.h:271
SparsityTools::partition
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
Definition sparsity_tools.cc:394
Utilities::MPI::min_max_avg
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm mpi_communicator)
Definition mpi.cc:77
Utilities::invert_permutation
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1670
internal::MatrixFreeFunctions::indicate_power_of_two
unsigned int indicate_power_of_two(const unsigned int vectorization_length)
Definition util.h:190
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
parallel::internal::parallel_for
void parallel_for(Iterator x_begin, Iterator x_end, const Functor &functor, const unsigned int grainsize)
Definition parallel.h:127
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
internal::MFWorkerInterface::face
virtual void face(const unsigned int range_index)=0
internal::MFWorkerInterface::zero_dst_vector_range
virtual void zero_dst_vector_range(const unsigned int range_index)=0
internal::MFWorkerInterface::cell
virtual void cell(const std::pair< unsigned int, unsigned int > &cell_range)=0
internal::MFWorkerInterface::boundary
virtual void boundary(const unsigned int range_index)=0
internal::MFWorkerInterface::vector_compress_start
virtual void vector_compress_start()=0
Starts the communication for the vector compress operation.
internal::MFWorkerInterface::cell_loop_post_range
virtual void cell_loop_post_range(const unsigned int range_index)=0
internal::MFWorkerInterface::vector_update_ghosts_start
virtual void vector_update_ghosts_start()=0
Starts the communication for the update ghost values operation.
internal::MFWorkerInterface::cell_loop_pre_range
virtual void cell_loop_pre_range(const unsigned int range_index)=0
internal::MFWorkerInterface::vector_update_ghosts_finish
virtual void vector_update_ghosts_finish()=0
Finishes the communication for the update ghost values operation.
internal::MFWorkerInterface::vector_compress_finish
virtual void vector_compress_finish()=0
Finishes the communication for the vector compress operation.
internal::MatrixFreeFunctions::TaskInfo::boundary_partition_data
std::vector< unsigned int > boundary_partition_data
Definition task_info.h:516
internal::MatrixFreeFunctions::TaskInfo::loop
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:348
internal::MatrixFreeFunctions::TaskInfo::partition_n_workers
std::vector< unsigned int > partition_n_workers
Definition task_info.h:574
internal::MatrixFreeFunctions::TaskInfo::create_blocks_serial
void create_blocks_serial(const std::vector< unsigned int > &cells_with_comm, const unsigned int dofs_per_cell, const bool categories_are_hp, const std::vector< unsigned int > &cell_vectorization_categories, const bool cell_vectorization_categories_strict, const std::vector< unsigned int > &parent_relation, std::vector< unsigned int > &renumbering, std::vector< unsigned char > &incompletely_filled_vectorization)
Definition task_info.cc:792
internal::MatrixFreeFunctions::TaskInfo::print_memory_statistics
void print_memory_statistics(StreamType &out, std::size_t data_length) const
Definition task_info.cc:688
internal::MatrixFreeFunctions::TaskInfo::partition_row_index
std::vector< unsigned int > partition_row_index
Definition task_info.h:464
internal::MatrixFreeFunctions::TaskInfo::partition_evens
std::vector< unsigned int > partition_evens
Definition task_info.h:556
internal::MatrixFreeFunctions::TaskInfo::partition_n_blocked_workers
std::vector< unsigned int > partition_n_blocked_workers
Definition task_info.h:568
internal::MatrixFreeFunctions::TaskInfo::cell_partition_data
std::vector< unsigned int > cell_partition_data
Definition task_info.h:472
internal::MatrixFreeFunctions::TaskInfo::make_boundary_cells_divisible
void make_boundary_cells_divisible(std::vector< unsigned int > &boundary_cells)
Definition task_info.cc:720
internal::MatrixFreeFunctions::TaskInfo::partition_odds
std::vector< unsigned int > partition_odds
Definition task_info.h:562
internal::MatrixFreeFunctions::TaskInfo::face_partition_data
std::vector< unsigned int > face_partition_data
Definition task_info.h:494
task_info.h
Generated by doxygen 1.9.8