Reference documentation for deal.II version Git b927e2c03c 2022-01-17 18:55:39 -0500
\(\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\}}\)
task_info.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
19 #include <deal.II/base/mpi.h>
21 #include <deal.II/base/parallel.h>
22 #include <deal.II/base/utilities.h>
23 
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 
56 
57 
58 
59 /*-------------------- Implementation of the matrix-free loop --------------*/
60 namespace 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  {
72  {
73  public:
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 
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:
110  const unsigned int partition;
112  };
113 
114  class CellWork : public tbb::task
115  {
116  public:
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:
140  const bool is_blocked;
141  };
142 
143 
144 
145  class PartitionWork : public tbb::task
146  {
147  public:
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 =
168  const unsigned int n_workers =
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,
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,
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,
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;
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:
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 =
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()),
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:
273  const unsigned int partition;
274  };
275 
276 
277 
278  class PartitionWork : public tbb::task
279  {
280  public:
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 =
298  1) /
300  parallel_for(tbb::blocked_range<unsigned int>(0, n_chunks, 1),
302  if (is_blocked == true)
303  tbb::empty_task::spawn(*dummy);
304  return nullptr;
305  }
306 
307  tbb::empty_task *dummy;
308 
309  private:
311  const unsigned int partition;
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)
333  else
335  return nullptr;
336  }
337 
338  private:
340  const bool do_compress;
341  };
342 
343 #endif // DEAL_II_WITH_TBB
344 
345 
346 
347  void
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)
356  else
357  funct.cell_loop_pre_range(
358  partition_row_index[partition_row_index.size() - 2]);
359 
361 
362 #if defined(DEAL_II_WITH_TBB) && !defined(DEAL_II_TBB_WITH_ONEAPI)
363 
364  if (scheme != none)
365  {
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())
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
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());
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);
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)
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)
643  else
644  funct.cell_loop_post_range(
645  partition_row_index[partition_row_index.size() - 2]);
646  }
647 
648 
649 
651  {
652  clear();
653  }
654 
655 
656 
657  void
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
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
704  {
705  return (
706  sizeof(*this) +
707  MemoryConsumption::memory_consumption(partition_row_index) +
708  MemoryConsumption::memory_consumption(cell_partition_data) +
709  MemoryConsumption::memory_consumption(face_partition_data) +
710  MemoryConsumption::memory_consumption(boundary_partition_data) +
711  MemoryConsumption::memory_consumption(partition_evens) +
712  MemoryConsumption::memory_consumption(partition_odds) +
713  MemoryConsumption::memory_consumption(partition_n_blocked_workers) +
714  MemoryConsumption::memory_consumption(partition_n_workers));
715  }
716 
717 
718 
719  void
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.size() == 0 ||
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
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  unsigned int vectorization_length_bits = 0;
827  unsigned int my_length = vectorization_length;
828  while ((my_length >>= 1) != 0u)
829  ++vectorization_length_bits;
830  const unsigned int n_lanes = 1 << vectorization_length_bits;
831 
832  // Step 1: find tight map of categories for not taking exceeding amounts
833  // of memory below. Sort the new categories by the numbers in the
834  // old one to ensure we respect the given rules
835  unsigned int n_categories = 1;
836  std::vector<unsigned int> tight_category_map(n_active_cells, 0);
837  if (cell_vectorization_categories.empty() == false)
838  {
839  AssertDimension(cell_vectorization_categories.size(),
840  n_active_cells + n_ghost_cells);
841 
842  std::set<unsigned int> used_categories;
843  for (unsigned int i = 0; i < n_active_cells; ++i)
844  used_categories.insert(cell_vectorization_categories[i]);
845  std::vector<unsigned int> used_categories_vector(
846  used_categories.size());
847  n_categories = 0;
848  for (const auto &it : used_categories)
849  used_categories_vector[n_categories++] = it;
850  for (unsigned int i = 0; i < n_active_cells; ++i)
851  {
852  const unsigned int index =
853  std::lower_bound(used_categories_vector.begin(),
854  used_categories_vector.end(),
855  cell_vectorization_categories[i]) -
856  used_categories_vector.begin();
857  AssertIndexRange(index, used_categories_vector.size());
858  tight_category_map[i] = index;
859  }
860  }
861 
862  // Step 2: Sort the cells by the category. If we want to fill up the
863  // ranges in vectorization, promote some of the cells to a higher
864  // category
865  std::vector<std::vector<unsigned int>> renumbering_category(n_categories);
866  for (unsigned int i = 0; i < n_active_cells; ++i)
867  renumbering_category[tight_category_map[i]].push_back(i);
868 
869  if (cell_vectorization_categories_strict == false && n_categories > 1)
870  for (unsigned int j = n_categories - 1; j > 0; --j)
871  {
872  unsigned int lower_index = j - 1;
873  while ((renumbering_category[j].size() % n_lanes) != 0u)
874  {
875  while (((renumbering_category[j].size() % n_lanes) != 0u) &&
876  !renumbering_category[lower_index].empty())
877  {
878  renumbering_category[j].push_back(
879  renumbering_category[lower_index].back());
880  renumbering_category[lower_index].pop_back();
881  }
882  if (lower_index == 0)
883  break;
884  else
885  --lower_index;
886  }
887  }
888 
889  // Step 3: Use the parent relation to find a good grouping of cells. To
890  // do this, we first put cells of each category defined above into two
891  // bins, those which we know can be grouped together by the given parent
892  // relation and those which cannot
893  std::vector<unsigned int> temporary_numbering;
894  temporary_numbering.reserve(n_active_cells +
895  (n_lanes - 1) * n_categories);
896  const unsigned int n_cells_per_parent =
897  std::count(parent_relation.begin(), parent_relation.end(), 0);
898  std::vector<unsigned int> category_size;
899  for (unsigned int j = 0; j < n_categories; ++j)
900  {
901  std::vector<std::pair<unsigned int, unsigned int>> grouped_cells;
902  std::vector<unsigned int> other_cells;
903  for (const unsigned int cell : renumbering_category[j])
904  if (parent_relation.empty() ||
905  parent_relation[cell] == numbers::invalid_unsigned_int)
906  other_cells.push_back(cell);
907  else
908  grouped_cells.emplace_back(parent_relation[cell], cell);
909 
910  // Compute the number of cells per group
911  std::sort(grouped_cells.begin(), grouped_cells.end());
912  std::vector<unsigned int> n_cells_per_group;
913  unsigned int length = 0;
914  for (unsigned int i = 0; i < grouped_cells.size(); ++i, ++length)
915  if (i > 0 && grouped_cells[i].first != grouped_cells[i - 1].first)
916  {
917  n_cells_per_group.push_back(length);
918  length = 0;
919  }
920  if (length > 0)
921  n_cells_per_group.push_back(length);
922 
923  // Move groups that do not have the complete size (due to
924  // categories) to the 'other_cells'. The cells with correct group
925  // size are immediately appended to the temporary cell numbering
926  auto group_it = grouped_cells.begin();
927  for (unsigned int length : n_cells_per_group)
928  if (length < n_cells_per_parent)
929  for (unsigned int j = 0; j < length; ++j)
930  other_cells.push_back((group_it++)->second);
931  else
932  {
933  // we should not have more cells in a group than in the first
934  // check we did above
935  AssertDimension(length, n_cells_per_parent);
936  for (unsigned int j = 0; j < length; ++j)
937  temporary_numbering.push_back((group_it++)->second);
938  }
939 
940  // Sort the remaining cells and append them as well
941  std::sort(other_cells.begin(), other_cells.end());
942  temporary_numbering.insert(temporary_numbering.end(),
943  other_cells.begin(),
944  other_cells.end());
945 
946  while (temporary_numbering.size() % n_lanes != 0)
947  temporary_numbering.push_back(numbers::invalid_unsigned_int);
948 
949  category_size.push_back(temporary_numbering.size());
950  }
951 
952  // Step 4: Identify the batches with cells marked as "comm"
953  std::vector<bool> batch_with_comm(temporary_numbering.size() / n_lanes,
954  false);
955  std::vector<unsigned int> temporary_numbering_inverse(n_active_cells);
956  for (unsigned int i = 0; i < temporary_numbering.size(); ++i)
957  if (temporary_numbering[i] != numbers::invalid_unsigned_int)
958  temporary_numbering_inverse[temporary_numbering[i]] = i;
959  for (const unsigned int cell : cells_with_comm)
960  batch_with_comm[temporary_numbering_inverse[cell] / n_lanes] = true;
961 
962  // Step 5: Sort the batches of cells by their last cell index to get
963  // good locality, assuming that the initial cell order is of good
964  // locality. In case we have hp-calculations with categories, we need to
965  // sort also by the category.
966  std::vector<std::array<unsigned int, 3>> batch_order;
967  std::vector<std::array<unsigned int, 3>> batch_order_comm;
968  for (unsigned int i = 0; i < temporary_numbering.size(); i += n_lanes)
969  {
970  unsigned int max_index = 0;
971  for (unsigned int j = 0; j < n_lanes; ++j)
972  if (temporary_numbering[i + j] < numbers::invalid_unsigned_int)
973  max_index = std::max(temporary_numbering[i + j], max_index);
974  const unsigned int category_hp =
975  categories_are_hp ?
976  std::upper_bound(category_size.begin(), category_size.end(), i) -
977  category_size.begin() :
978  0;
979  const std::array<unsigned int, 3> next{{category_hp, max_index, i}};
980  if (batch_with_comm[i / n_lanes])
981  batch_order_comm.emplace_back(next);
982  else
983  batch_order.emplace_back(next);
984  }
985 
986  std::sort(batch_order.begin(), batch_order.end());
987  std::sort(batch_order_comm.begin(), batch_order_comm.end());
988 
989  // Step 6: Put the cells with communication in the middle of the cell
990  // range. For the MPI case, we need three groups to enable overlap for
991  // communication and computation (part before comm, part with comm, part
992  // after comm), whereas we need one for the other case. And in each
993  // case, we allow for a slot of "ghosted" cells.
994  std::vector<unsigned int> blocks;
995  if (n_procs == 1)
996  {
997  if (batch_order.empty())
998  std::swap(batch_order_comm, batch_order);
999  Assert(batch_order_comm.empty(), ExcInternalError());
1000  partition_row_index.resize(3);
1001  blocks = {0, static_cast<unsigned int>(batch_order.size())};
1002  }
1003  else
1004  {
1005  partition_row_index.resize(5);
1006  const unsigned int comm_begin = batch_order.size() / 2;
1007  batch_order.insert(batch_order.begin() + comm_begin,
1008  batch_order_comm.begin(),
1009  batch_order_comm.end());
1010  const unsigned int comm_end = comm_begin + batch_order_comm.size();
1011  const unsigned int end = batch_order.size();
1012  blocks = {0, comm_begin, comm_end, end};
1013  }
1014 
1015  // Step 7: Fill in the data by batches for the locally owned cells.
1016  const unsigned int n_cell_batches = batch_order.size();
1017  const unsigned int n_ghost_batches =
1018  (n_ghost_cells + n_lanes - 1) / n_lanes;
1019  incompletely_filled_vectorization.resize(n_cell_batches +
1020  n_ghost_batches);
1021 
1022  cell_partition_data.clear();
1023  cell_partition_data.resize(1, 0);
1024 
1025  renumbering.clear();
1026  renumbering.resize(n_active_cells + n_ghost_cells,
1028 
1029  unsigned int counter = 0;
1030  for (unsigned int block = 0; block < blocks.size() - 1; ++block)
1031  {
1032  const unsigned int grain_size =
1033  std::max((2048U / dofs_per_cell) / 8 * 4, 2U);
1034  for (unsigned int k = blocks[block]; k < blocks[block + 1];
1035  k += grain_size)
1036  cell_partition_data.push_back(
1037  std::min(k + grain_size, blocks[block + 1]));
1038  partition_row_index[block + 1] = cell_partition_data.size() - 1;
1039 
1040  // Set the numbering according to the reordered temporary one
1041  for (unsigned int k = blocks[block]; k < blocks[block + 1]; ++k)
1042  {
1043  const unsigned int pos = batch_order[k][2];
1044  unsigned int j = 0;
1045  for (; j < n_lanes && temporary_numbering[pos + j] !=
1047  ++j)
1048  renumbering[counter++] = temporary_numbering[pos + j];
1049  if (j < n_lanes)
1050  incompletely_filled_vectorization[k] = j;
1051  }
1052  }
1053  AssertDimension(counter, n_active_cells);
1054 
1055  // Step 8: Treat the ghost cells
1056  for (unsigned int cell = n_active_cells;
1057  cell < n_active_cells + n_ghost_cells;
1058  ++cell)
1059  {
1060  if (!cell_vectorization_categories.empty())
1061  AssertDimension(cell_vectorization_categories[cell],
1062  cell_vectorization_categories[n_active_cells]);
1063  renumbering[cell] = cell;
1064  }
1065  if ((n_ghost_cells % n_lanes) != 0u)
1066  incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes;
1067  cell_partition_data.push_back(n_cell_batches + n_ghost_batches);
1068  partition_row_index.back() = cell_partition_data.size() - 1;
1069 
1070 #ifdef DEBUG
1071  std::vector<unsigned int> renumber_cpy(renumbering);
1072  std::sort(renumber_cpy.begin(), renumber_cpy.end());
1073  for (unsigned int i = 0; i < renumber_cpy.size(); ++i)
1074  AssertDimension(i, renumber_cpy[i]);
1075 #endif
1076  }
1077 
1078 
1079 
1080  void
1082  const std::vector<unsigned int> &boundary_cells,
1083  std::vector<unsigned int> & renumbering,
1084  std::vector<unsigned char> & incompletely_filled_vectorization)
1085  {
1086  const unsigned int n_cell_batches =
1087  (n_active_cells + vectorization_length - 1) / vectorization_length;
1088  const unsigned int n_ghost_slots =
1089  (n_ghost_cells + vectorization_length - 1) / vectorization_length;
1090  incompletely_filled_vectorization.resize(n_cell_batches + n_ghost_slots);
1091  if (n_cell_batches * vectorization_length > n_active_cells)
1092  incompletely_filled_vectorization[n_cell_batches - 1] =
1093  vectorization_length -
1094  (n_cell_batches * vectorization_length - n_active_cells);
1095  if (n_ghost_slots * vectorization_length > n_ghost_cells)
1096  incompletely_filled_vectorization[n_cell_batches + n_ghost_slots - 1] =
1097  vectorization_length -
1098  (n_ghost_slots * vectorization_length - n_ghost_cells);
1099 
1100  std::vector<unsigned int> reverse_numbering(
1102  for (unsigned int j = 0; j < boundary_cells.size(); ++j)
1103  reverse_numbering[boundary_cells[j]] = j;
1104  unsigned int counter = boundary_cells.size();
1105  for (unsigned int j = 0; j < n_active_cells; ++j)
1106  if (reverse_numbering[j] == numbers::invalid_unsigned_int)
1107  reverse_numbering[j] = counter++;
1108 
1109  AssertDimension(counter, n_active_cells);
1110  renumbering = Utilities::invert_permutation(reverse_numbering);
1111 
1112  for (unsigned int j = n_active_cells; j < n_active_cells + n_ghost_cells;
1113  ++j)
1114  renumbering.push_back(j);
1115 
1116  // TODO: might be able to simplify this code by not relying on the cell
1117  // partition data while computing the thread graph
1118  cell_partition_data.clear();
1119  cell_partition_data.push_back(0);
1120  if (n_procs > 1)
1121  {
1122  const unsigned int n_macro_boundary_cells =
1123  (boundary_cells.size() + vectorization_length - 1) /
1124  vectorization_length;
1125  cell_partition_data.push_back(
1126  (n_cell_batches - n_macro_boundary_cells) / 2);
1127  cell_partition_data.push_back(cell_partition_data[1] +
1128  n_macro_boundary_cells);
1129  }
1130  else
1131  AssertDimension(boundary_cells.size(), 0);
1132  cell_partition_data.push_back(n_cell_batches);
1133  cell_partition_data.push_back(cell_partition_data.back() + n_ghost_slots);
1134  partition_row_index.resize(n_procs > 1 ? 4 : 2);
1135  partition_row_index[0] = 0;
1136  partition_row_index[1] = 1;
1137  if (n_procs > 1)
1138  {
1139  partition_row_index[2] = 2;
1140  partition_row_index[3] = 3;
1141  }
1142  }
1143 
1144 
1145 
1146  void
1147  TaskInfo::guess_block_size(const unsigned int dofs_per_cell)
1148  {
1149  // user did not say a positive number, so we have to guess
1150  if (block_size == 0)
1151  {
1152  // we would like to have enough work to do, so as first guess, try
1153  // to get 16 times as many chunks as we have threads on the system.
1155  vectorization_length);
1156 
1157  // if there are too few degrees of freedom per cell, need to
1158  // increase the block size
1159  const unsigned int minimum_parallel_grain_size = 200;
1160  if (dofs_per_cell * block_size < minimum_parallel_grain_size)
1161  block_size = (minimum_parallel_grain_size / dofs_per_cell + 1);
1162  if (dofs_per_cell * block_size > 10000)
1163  block_size /= 4;
1164 
1165  block_size =
1166  1 << static_cast<unsigned int>(std::log2(block_size + 1));
1167  }
1168  if (block_size > n_active_cells)
1170  }
1171 
1172 
1173 
1174  void
1176  DynamicSparsityPattern & connectivity_large,
1177  std::vector<unsigned int> & renumbering,
1178  std::vector<unsigned char> &irregular_cells,
1179  const bool)
1180  {
1181  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1182  if (n_cell_batches == 0)
1183  return;
1184 
1185  Assert(vectorization_length > 0, ExcInternalError());
1186 
1187  unsigned int partition = 0, counter = 0;
1188 
1189  // Create connectivity graph for blocks based on connectivity graph for
1190  // cells.
1191  DynamicSparsityPattern connectivity(n_blocks, n_blocks);
1192  make_connectivity_cells_to_blocks(irregular_cells,
1193  connectivity_large,
1194  connectivity);
1195 
1196  // Create cell-block partitioning.
1197 
1198  // For each block of cells, this variable saves to which partitions the
1199  // block belongs. Initialize all to -1 to mark them as not yet assigned
1200  // a partition.
1201  std::vector<unsigned int> cell_partition(n_blocks,
1203 
1204  // In element j of this variable, one puts the old number of the block
1205  // that should be the jth block in the new numeration.
1206  std::vector<unsigned int> partition_list(n_blocks, 0);
1207  std::vector<unsigned int> partition_color_list(n_blocks, 0);
1208 
1209  // This vector points to the start of each partition.
1210  std::vector<unsigned int> partition_size(2, 0);
1211 
1212  // blocking_connectivity = true;
1213 
1214  // The cluster_size in make_partitioning defines that the no. of cells
1215  // in each partition should be a multiple of cluster_size.
1216  unsigned int cluster_size = 1;
1217 
1218  // Make the partitioning of the first layer of the blocks of cells.
1219  make_partitioning(connectivity,
1220  cluster_size,
1221  cell_partition,
1222  partition_list,
1223  partition_size,
1224  partition);
1225 
1226  // Color the cells within each partition
1227  make_coloring_within_partitions_pre_blocked(connectivity,
1228  partition,
1229  cell_partition,
1230  partition_list,
1231  partition_size,
1232  partition_color_list);
1233 
1234  partition_list = renumbering;
1235 
1236 #ifdef DEBUG
1237  // in debug mode, check that the partition color list is one-to-one
1238  {
1239  std::vector<unsigned int> sorted_pc_list(partition_color_list);
1240  std::sort(sorted_pc_list.begin(), sorted_pc_list.end());
1241  for (unsigned int i = 0; i < sorted_pc_list.size(); ++i)
1242  Assert(sorted_pc_list[i] == i, ExcInternalError());
1243  }
1244 #endif
1245 
1246  // set the start list for each block and compute the renumbering of
1247  // cells
1248  std::vector<unsigned int> block_start(n_cell_batches + 1);
1249  std::vector<unsigned char> irregular(n_cell_batches);
1250 
1251  unsigned int mcell_start = 0;
1252  block_start[0] = 0;
1253  for (unsigned int block = 0; block < n_blocks; ++block)
1254  {
1255  block_start[block + 1] = block_start[block];
1256  for (unsigned int mcell = mcell_start;
1257  mcell < std::min(mcell_start + block_size, n_cell_batches);
1258  ++mcell)
1259  {
1260  unsigned int n_comp = (irregular_cells[mcell] > 0) ?
1261  irregular_cells[mcell] :
1262  vectorization_length;
1263  block_start[block + 1] += n_comp;
1264  ++counter;
1265  }
1266  mcell_start += block_size;
1267  }
1268  counter = 0;
1269  unsigned int counter_macro = 0;
1270  unsigned int block_size_last =
1271  n_cell_batches - block_size * (n_blocks - 1);
1272  if (block_size_last == 0)
1273  block_size_last = block_size;
1274 
1275  unsigned int tick = 0;
1276  for (unsigned int block = 0; block < n_blocks; ++block)
1277  {
1278  unsigned int present_block = partition_color_list[block];
1279  for (unsigned int cell = block_start[present_block];
1280  cell < block_start[present_block + 1];
1281  ++cell)
1282  renumbering[counter++] = partition_list[cell];
1283  unsigned int this_block_size =
1284  (present_block == n_blocks - 1) ? block_size_last : block_size;
1285 
1286  // Also re-compute the content of cell_partition_data to
1287  // contain the numbers of cells, not blocks
1288  if (cell_partition_data[tick] == block)
1289  cell_partition_data[tick++] = counter_macro;
1290 
1291  for (unsigned int j = 0; j < this_block_size; ++j)
1292  irregular[counter_macro++] =
1293  irregular_cells[present_block * block_size + j];
1294  }
1295  AssertDimension(tick + 1, cell_partition_data.size());
1296  cell_partition_data.back() = counter_macro;
1297 
1298  irregular_cells.swap(irregular);
1299  AssertDimension(counter, n_active_cells);
1300  AssertDimension(counter_macro, n_cell_batches);
1301 
1302  // check that the renumbering is one-to-one
1303 #ifdef DEBUG
1304  {
1305  std::vector<unsigned int> sorted_renumbering(renumbering);
1306  std::sort(sorted_renumbering.begin(), sorted_renumbering.end());
1307  for (unsigned int i = 0; i < sorted_renumbering.size(); ++i)
1308  Assert(sorted_renumbering[i] == i, ExcInternalError());
1309  }
1310 #endif
1311 
1312 
1313  update_task_info(
1314  partition); // Actually sets too much for partition color case
1315 
1316  AssertDimension(cell_partition_data.back(), n_cell_batches);
1317  }
1318 
1319 
1320 
1321  void
1323  const std::vector<unsigned int> &cell_active_fe_index,
1324  DynamicSparsityPattern & connectivity,
1325  std::vector<unsigned int> & renumbering,
1326  std::vector<unsigned char> & irregular_cells,
1327  const bool hp_bool)
1328  {
1329  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1330  if (n_cell_batches == 0)
1331  return;
1332 
1333  Assert(vectorization_length > 0, ExcInternalError());
1334 
1335  // if we want to block before partitioning, create connectivity graph
1336  // for blocks based on connectivity graph for cells.
1337  DynamicSparsityPattern connectivity_blocks(n_blocks, n_blocks);
1338  make_connectivity_cells_to_blocks(irregular_cells,
1339  connectivity,
1340  connectivity_blocks);
1341 
1342  unsigned int n_blocks = 0;
1343  if (scheme == partition_color ||
1344  scheme == color) // blocking_connectivity == true
1345  n_blocks = this->n_blocks;
1346  else
1347  n_blocks = n_active_cells;
1348 
1349  // For each block of cells, this variable saves to which partitions the
1350  // block belongs. Initialize all to -1 to mark them as not yet assigned
1351  // a partition.
1352  std::vector<unsigned int> cell_partition(n_blocks,
1354 
1355  // In element j of this variable, one puts the old number (but after
1356  // renumbering according to the input renumbering) of the block that
1357  // should be the jth block in the new numeration.
1358  std::vector<unsigned int> partition_list(n_blocks, 0);
1359  std::vector<unsigned int> partition_2layers_list(n_blocks, 0);
1360 
1361  // This vector points to the start of each partition.
1362  std::vector<unsigned int> partition_size(2, 0);
1363 
1364  unsigned int partition = 0;
1365 
1366  // Within the partitions we want to be able to block for the case that
1367  // we do not block already in the connectivity. The cluster_size in
1368  // make_partitioning defines that the no. of cells in each partition
1369  // should be a multiple of cluster_size.
1370  unsigned int cluster_size = 1;
1371  if (scheme == partition_partition)
1372  cluster_size = block_size * vectorization_length;
1373 
1374  // Make the partitioning of the first layer of the blocks of cells.
1375  if (scheme == partition_color || scheme == color)
1376  make_partitioning(connectivity_blocks,
1377  cluster_size,
1378  cell_partition,
1379  partition_list,
1380  partition_size,
1381  partition);
1382  else
1383  make_partitioning(connectivity,
1384  cluster_size,
1385  cell_partition,
1386  partition_list,
1387  partition_size,
1388  partition);
1389 
1390  // Partition or color second layer
1391  if (scheme == partition_partition)
1392 
1393  {
1394  // Partition within partitions.
1395  make_partitioning_within_partitions_post_blocked(
1396  connectivity,
1397  cell_active_fe_index,
1398  partition,
1399  cluster_size,
1400  hp_bool,
1401  cell_partition,
1402  partition_list,
1403  partition_size,
1404  partition_2layers_list,
1405  irregular_cells);
1406  }
1407  else if (scheme == partition_color || scheme == color)
1408  {
1409  make_coloring_within_partitions_pre_blocked(connectivity_blocks,
1410  partition,
1411  cell_partition,
1412  partition_list,
1413  partition_size,
1414  partition_2layers_list);
1415  }
1416 
1417  // in debug mode, check that the partition_2layers_list is one-to-one
1418 #ifdef DEBUG
1419  {
1420  std::vector<unsigned int> sorted_pc_list(partition_2layers_list);
1421  std::sort(sorted_pc_list.begin(), sorted_pc_list.end());
1422  for (unsigned int i = 0; i < sorted_pc_list.size(); ++i)
1423  Assert(sorted_pc_list[i] == i, ExcInternalError());
1424  }
1425 #endif
1426 
1427  // Set the new renumbering
1428  std::vector<unsigned int> renumbering_in(n_active_cells, 0);
1429  renumbering_in.swap(renumbering);
1430  if (scheme == partition_partition) // blocking_connectivity == false
1431  {
1432  // This is the simple case. The renumbering is just a combination of
1433  // the renumbering that we were given as an input and the
1434  // renumbering of partition/coloring given in partition_2layers_list
1435  for (unsigned int j = 0; j < renumbering.size(); ++j)
1436  renumbering[j] = renumbering_in[partition_2layers_list[j]];
1437  // Account for the ghost cells, finally.
1438  for (unsigned int i = 0; i < n_ghost_cells; ++i)
1439  renumbering.push_back(i + n_active_cells);
1440  }
1441  else
1442  {
1443  // set the start list for each block and compute the renumbering of
1444  // cells
1445  std::vector<unsigned int> block_start(n_cell_batches + 1);
1446  std::vector<unsigned char> irregular(n_cell_batches);
1447 
1448  unsigned int counter = 0;
1449  unsigned int mcell_start = 0;
1450  block_start[0] = 0;
1451  for (unsigned int block = 0; block < n_blocks; ++block)
1452  {
1453  block_start[block + 1] = block_start[block];
1454  for (unsigned int mcell = mcell_start;
1455  mcell < std::min(mcell_start + block_size, n_cell_batches);
1456  ++mcell)
1457  {
1458  unsigned int n_comp = (irregular_cells[mcell] > 0) ?
1459  irregular_cells[mcell] :
1460  vectorization_length;
1461  block_start[block + 1] += n_comp;
1462  ++counter;
1463  }
1464  mcell_start += block_size;
1465  }
1466  counter = 0;
1467  unsigned int counter_macro = 0;
1468  unsigned int block_size_last =
1469  n_cell_batches - block_size * (n_blocks - 1);
1470  if (block_size_last == 0)
1471  block_size_last = block_size;
1472 
1473  unsigned int tick = 0;
1474  for (unsigned int block = 0; block < n_blocks; ++block)
1475  {
1476  unsigned int present_block = partition_2layers_list[block];
1477  for (unsigned int cell = block_start[present_block];
1478  cell < block_start[present_block + 1];
1479  ++cell)
1480  renumbering[counter++] = renumbering_in[cell];
1481  unsigned int this_block_size =
1482  (present_block == n_blocks - 1) ? block_size_last : block_size;
1483 
1484  // Also re-compute the content of cell_partition_data to
1485  // contain the numbers of cells, not blocks
1486  if (cell_partition_data[tick] == block)
1487  cell_partition_data[tick++] = counter_macro;
1488 
1489  for (unsigned int j = 0; j < this_block_size; ++j)
1490  irregular[counter_macro++] =
1491  irregular_cells[present_block * block_size + j];
1492  }
1493  AssertDimension(tick + 1, cell_partition_data.size());
1494  cell_partition_data.back() = counter_macro;
1495 
1496  irregular_cells.swap(irregular);
1497  AssertDimension(counter, n_active_cells);
1498  AssertDimension(counter_macro, n_cell_batches);
1499  // check that the renumbering is one-to-one
1500 #ifdef DEBUG
1501  {
1502  std::vector<unsigned int> sorted_renumbering(renumbering);
1503  std::sort(sorted_renumbering.begin(), sorted_renumbering.end());
1504  for (unsigned int i = 0; i < sorted_renumbering.size(); ++i)
1505  Assert(sorted_renumbering[i] == i, ExcInternalError());
1506  }
1507 #endif
1508  }
1509 
1510  // Update the task_info with the more information for the thread graph.
1511  update_task_info(partition);
1512  }
1513 
1514 
1515 
1516  void
1518  const std::vector<unsigned int> &cell_active_fe_index,
1519  DynamicSparsityPattern & connectivity,
1520  std::vector<unsigned int> & renumbering,
1521  std::vector<unsigned char> & irregular_cells,
1522  const bool hp_bool)
1523  {
1524  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1525  if (n_cell_batches == 0)
1526  return;
1527 
1528  const unsigned int cluster_size = block_size * vectorization_length;
1529 
1530  // Create cell-block partitioning.
1531 
1532  // For each block of cells, this variable saves to which partitions the
1533  // block belongs. Initialize all to n_cell_batches to mark them as not
1534  // yet assigned a partition.
1535  std::vector<unsigned int> cell_partition(n_active_cells,
1537 
1538 
1539  // In element j of this variable, one puts the old number of the block
1540  // that should be the jth block in the new numeration.
1541  std::vector<unsigned int> partition_list(n_active_cells, 0);
1542  std::vector<unsigned int> partition_partition_list(n_active_cells, 0);
1543 
1544  // This vector points to the start of each partition.
1545  std::vector<unsigned int> partition_size(2, 0);
1546 
1547  unsigned int partition = 0;
1548  // Here, we do not block inside the connectivity graph
1549  // blocking_connectivity = false;
1550 
1551  // Make the partitioning of the first layer of the blocks of cells.
1552  make_partitioning(connectivity,
1553  cluster_size,
1554  cell_partition,
1555  partition_list,
1556  partition_size,
1557  partition);
1558 
1559  // Partition within partitions.
1560  make_partitioning_within_partitions_post_blocked(connectivity,
1561  cell_active_fe_index,
1562  partition,
1563  cluster_size,
1564  hp_bool,
1565  cell_partition,
1566  partition_list,
1567  partition_size,
1568  partition_partition_list,
1569  irregular_cells);
1570 
1571  partition_list.swap(renumbering);
1572 
1573  for (unsigned int j = 0; j < renumbering.size(); ++j)
1574  renumbering[j] = partition_list[partition_partition_list[j]];
1575 
1576  for (unsigned int i = 0; i < n_ghost_cells; ++i)
1577  renumbering.push_back(i + n_active_cells);
1578 
1579  update_task_info(partition);
1580  }
1581 
1582 
1583 
1584  void
1586  const std::vector<unsigned char> &irregular_cells,
1587  const DynamicSparsityPattern & connectivity_cells,
1588  DynamicSparsityPattern & connectivity_blocks) const
1589  {
1590  std::vector<std::vector<unsigned int>> cell_blocks(n_blocks);
1591  std::vector<unsigned int> touched_cells(n_active_cells);
1592  unsigned int cell = 0;
1593  for (unsigned int i = 0, mcell = 0; i < n_blocks; ++i)
1594  {
1595  for (unsigned int c = 0;
1596  c < block_size && mcell < *(cell_partition_data.end() - 2);
1597  ++c, ++mcell)
1598  {
1599  unsigned int ncomp = (irregular_cells[mcell] > 0) ?
1600  irregular_cells[mcell] :
1601  vectorization_length;
1602  for (unsigned int c = 0; c < ncomp; ++c, ++cell)
1603  {
1604  cell_blocks[i].push_back(cell);
1605  touched_cells[cell] = i;
1606  }
1607  }
1608  }
1610  for (unsigned int i = 0; i < cell_blocks.size(); ++i)
1611  for (unsigned int col = 0; col < cell_blocks[i].size(); ++col)
1612  {
1614  connectivity_cells.begin(cell_blocks[i][col]);
1615  it != connectivity_cells.end(cell_blocks[i][col]);
1616  ++it)
1617  {
1618  if (touched_cells[it->column()] != i)
1619  connectivity_blocks.add(i, touched_cells[it->column()]);
1620  }
1621  }
1622  }
1623 
1624 
1625 
1626  // Function to create partitioning on the second layer within each
1627  // partition. Version without preblocking.
1628  void
1630  const DynamicSparsityPattern & connectivity,
1631  const std::vector<unsigned int> &cell_active_fe_index,
1632  const unsigned int partition,
1633  const unsigned int cluster_size,
1634  const bool hp_bool,
1635  const std::vector<unsigned int> &cell_partition,
1636  const std::vector<unsigned int> &partition_list,
1637  const std::vector<unsigned int> &partition_size,
1638  std::vector<unsigned int> & partition_partition_list,
1639  std::vector<unsigned char> & irregular_cells)
1640  {
1641  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1642  const unsigned int n_ghost_slots =
1643  *(cell_partition_data.end() - 1) - n_cell_batches;
1644 
1645  // List of cells in previous partition
1646  std::vector<unsigned int> neighbor_list;
1647  // List of cells in current partition for use as neighbors in next
1648  // partition
1649  std::vector<unsigned int> neighbor_neighbor_list;
1650 
1651  std::vector<unsigned int> renumbering(n_active_cells);
1652 
1653  irregular_cells.back() = 0;
1654  irregular_cells.resize(n_active_cells + n_ghost_slots);
1655 
1656  unsigned int max_fe_index = 0;
1657  for (const unsigned int fe_index : cell_active_fe_index)
1658  max_fe_index = std::max(fe_index, max_fe_index);
1659 
1660  Assert(!hp_bool || cell_active_fe_index.size() == n_active_cells,
1661  ExcInternalError());
1662 
1663  {
1664  unsigned int n_cell_batches_before = 0;
1665  // Create partitioning within partitions.
1666 
1667  // For each block of cells, this variable saves to which partitions
1668  // the block belongs. Initialize all to n_cell_batches to mark them as
1669  // not yet assigned a partition.
1670  std::vector<unsigned int> cell_partition_l2(
1672  partition_row_index.clear();
1673  partition_row_index.resize(partition + 1, 0);
1674  cell_partition_data.resize(1, 0);
1675 
1676  unsigned int counter = 0;
1677  unsigned int missing_macros;
1678  for (unsigned int part = 0; part < partition; ++part)
1679  {
1680  neighbor_neighbor_list.resize(0);
1681  neighbor_list.resize(0);
1682  bool work = true;
1683  unsigned int partition_l2 = 0;
1684  unsigned int start_up = partition_size[part];
1685  unsigned int partition_counter = 0;
1686  while (work)
1687  {
1688  if (neighbor_list.size() == 0)
1689  {
1690  work = false;
1691  partition_counter = 0;
1692  for (unsigned int j = start_up;
1693  j < partition_size[part + 1];
1694  ++j)
1695  if (cell_partition[partition_list[j]] == part &&
1696  cell_partition_l2[partition_list[j]] ==
1698  {
1699  start_up = j;
1700  work = true;
1701  partition_counter = 1;
1702  // To start up, set the start_up cell to partition
1703  // and list all its neighbors.
1704  AssertIndexRange(start_up, partition_size[part + 1]);
1705  cell_partition_l2[partition_list[start_up]] =
1706  partition_l2;
1707  neighbor_neighbor_list.push_back(
1708  partition_list[start_up]);
1709  partition_partition_list[counter++] =
1710  partition_list[start_up];
1711  start_up++;
1712  break;
1713  }
1714  }
1715  else
1716  {
1717  partition_counter = 0;
1718  for (const unsigned int neighbor : neighbor_list)
1719  {
1720  Assert(cell_partition[neighbor] == part,
1721  ExcInternalError());
1722  Assert(cell_partition_l2[neighbor] == partition_l2 - 1,
1723  ExcInternalError());
1724  auto neighbor_it = connectivity.begin(neighbor);
1725  const auto end_it = connectivity.end(neighbor);
1726  for (; neighbor_it != end_it; ++neighbor_it)
1727  {
1728  if (cell_partition[neighbor_it->column()] == part &&
1729  cell_partition_l2[neighbor_it->column()] ==
1731  {
1732  cell_partition_l2[neighbor_it->column()] =
1733  partition_l2;
1734  neighbor_neighbor_list.push_back(
1735  neighbor_it->column());
1736  partition_partition_list[counter++] =
1737  neighbor_it->column();
1738  partition_counter++;
1739  }
1740  }
1741  }
1742  }
1743  if (partition_counter > 0)
1744  {
1745  int index_before = neighbor_neighbor_list.size(),
1746  index = index_before;
1747  {
1748  // put the cells into separate lists for each FE index
1749  // within one partition-partition
1750  missing_macros = 0;
1751  std::vector<unsigned int> remaining_per_cell_batch(
1752  max_fe_index + 1);
1753  std::vector<std::vector<unsigned int>>
1754  renumbering_fe_index;
1755  unsigned int cell;
1756  bool filled = true;
1757  if (hp_bool == true)
1758  {
1759  renumbering_fe_index.resize(max_fe_index + 1);
1760  for (cell = counter - partition_counter;
1761  cell < counter;
1762  ++cell)
1763  {
1764  renumbering_fe_index
1765  [cell_active_fe_index.empty() ?
1766  0 :
1767  cell_active_fe_index
1768  [partition_partition_list[cell]]]
1769  .push_back(partition_partition_list[cell]);
1770  }
1771  // check how many more cells are needed in the lists
1772  for (unsigned int j = 0; j < max_fe_index + 1; ++j)
1773  {
1774  remaining_per_cell_batch[j] =
1775  renumbering_fe_index[j].size() %
1776  vectorization_length;
1777  if (remaining_per_cell_batch[j] != 0)
1778  filled = false;
1779  missing_macros +=
1780  ((renumbering_fe_index[j].size() +
1781  vectorization_length - 1) /
1782  vectorization_length);
1783  }
1784  }
1785  else
1786  {
1787  remaining_per_cell_batch.resize(1);
1788  remaining_per_cell_batch[0] =
1789  partition_counter % vectorization_length;
1790  missing_macros =
1791  partition_counter / vectorization_length;
1792  if (remaining_per_cell_batch[0] != 0)
1793  {
1794  filled = false;
1795  missing_macros++;
1796  }
1797  }
1798  missing_macros =
1799  cluster_size - (missing_macros % cluster_size);
1800 
1801  // now we realized that there are some cells missing.
1802  while (missing_macros > 0 || filled == false)
1803  {
1804  if (index == 0)
1805  {
1806  index = neighbor_neighbor_list.size();
1807  if (index == index_before)
1808  {
1809  if (missing_macros != 0)
1810  {
1811  neighbor_neighbor_list.resize(0);
1812  }
1813  start_up--;
1814  break; // not connected - start again
1815  }
1816  index_before = index;
1817  }
1818  index--;
1819  unsigned int additional =
1820  neighbor_neighbor_list[index];
1821 
1822  // go through the neighbors of the last cell in the
1823  // current partition and check if we find some to
1824  // fill up with.
1826  connectivity.begin(
1827  additional),
1828  end =
1829  connectivity.end(
1830  additional);
1831  for (; neighbor != end; ++neighbor)
1832  {
1833  if (cell_partition[neighbor->column()] == part &&
1834  cell_partition_l2[neighbor->column()] ==
1836  {
1837  unsigned int this_index = 0;
1838  if (hp_bool == true)
1839  this_index =
1840  cell_active_fe_index.empty() ?
1841  0 :
1842  cell_active_fe_index[neighbor
1843  ->column()];
1844 
1845  // Only add this cell if we need more macro
1846  // cells in the current block or if there is
1847  // a macro cell with the FE index that is
1848  // not yet fully populated
1849  if (missing_macros > 0 ||
1850  remaining_per_cell_batch[this_index] > 0)
1851  {
1852  cell_partition_l2[neighbor->column()] =
1853  partition_l2;
1854  neighbor_neighbor_list.push_back(
1855  neighbor->column());
1856  if (hp_bool == true)
1857  renumbering_fe_index[this_index]
1858  .push_back(neighbor->column());
1859  partition_partition_list[counter] =
1860  neighbor->column();
1861  counter++;
1862  partition_counter++;
1863  if (remaining_per_cell_batch
1864  [this_index] == 0 &&
1865  missing_macros > 0)
1866  missing_macros--;
1867  remaining_per_cell_batch[this_index]++;
1868  if (remaining_per_cell_batch
1869  [this_index] ==
1870  vectorization_length)
1871  {
1872  remaining_per_cell_batch[this_index] =
1873  0;
1874  }
1875  if (missing_macros == 0)
1876  {
1877  filled = true;
1878  for (unsigned int fe_ind = 0;
1879  fe_ind < max_fe_index + 1;
1880  ++fe_ind)
1881  if (remaining_per_cell_batch
1882  [fe_ind] != 0)
1883  filled = false;
1884  }
1885  if (filled == true)
1886  break;
1887  }
1888  }
1889  }
1890  }
1891  if (hp_bool == true)
1892  {
1893  // set the renumbering according to their active FE
1894  // index within one partition-partition which was
1895  // implicitly assumed above
1896  cell = counter - partition_counter;
1897  for (unsigned int j = 0; j < max_fe_index + 1; ++j)
1898  {
1899  for (const unsigned int jj :
1900  renumbering_fe_index[j])
1901  renumbering[cell++] = jj;
1902  if (renumbering_fe_index[j].size() %
1903  vectorization_length !=
1904  0)
1905  irregular_cells[renumbering_fe_index[j].size() /
1906  vectorization_length +
1907  n_cell_batches_before] =
1908  renumbering_fe_index[j].size() %
1909  vectorization_length;
1910  n_cell_batches_before +=
1911  (renumbering_fe_index[j].size() +
1912  vectorization_length - 1) /
1913  vectorization_length;
1914  renumbering_fe_index[j].resize(0);
1915  }
1916  }
1917  else
1918  {
1919  n_cell_batches_before +=
1920  partition_counter / vectorization_length;
1921  if (partition_counter % vectorization_length != 0)
1922  {
1923  irregular_cells[n_cell_batches_before] =
1924  partition_counter % vectorization_length;
1925  n_cell_batches_before++;
1926  }
1927  }
1928  }
1929  cell_partition_data.push_back(n_cell_batches_before);
1930  partition_l2++;
1931  }
1932  neighbor_list = neighbor_neighbor_list;
1933  neighbor_neighbor_list.resize(0);
1934  }
1935  partition_row_index[part + 1] =
1936  partition_row_index[part] + partition_l2;
1937  }
1938  }
1939  if (hp_bool == true)
1940  {
1941  partition_partition_list.swap(renumbering);
1942  }
1943  }
1944 
1945 
1946 
1947  // Function to create coloring on the second layer within each partition.
1948  // Version assumes preblocking.
1949  void
1951  const DynamicSparsityPattern & connectivity,
1952  const unsigned int partition,
1953  const std::vector<unsigned int> &cell_partition,
1954  const std::vector<unsigned int> &partition_list,
1955  const std::vector<unsigned int> &partition_size,
1956  std::vector<unsigned int> & partition_color_list)
1957  {
1958  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
1959  std::vector<unsigned int> cell_color(n_blocks, n_cell_batches);
1960  std::vector<bool> color_finder;
1961 
1962  partition_row_index.resize(partition + 1);
1963  cell_partition_data.clear();
1964  unsigned int color_counter = 0, index_counter = 0;
1965  for (unsigned int part = 0; part < partition; ++part)
1966  {
1967  partition_row_index[part] = index_counter;
1968  unsigned int max_color = 0;
1969  for (unsigned int k = partition_size[part];
1970  k < partition_size[part + 1];
1971  k++)
1972  {
1973  unsigned int cell = partition_list[k];
1974  unsigned int n_neighbors = connectivity.row_length(cell);
1975 
1976  // In the worst case, each neighbor has a different color. So we
1977  // find at least one available color between 0 and n_neighbors.
1978  color_finder.resize(n_neighbors + 1);
1979  for (unsigned int j = 0; j <= n_neighbors; ++j)
1980  color_finder[j] = true;
1982  connectivity.begin(cell),
1983  end = connectivity.end(cell);
1984  for (; neighbor != end; ++neighbor)
1985  {
1986  // Mark the color that a neighbor within the partition has
1987  // as taken
1988  if (cell_partition[neighbor->column()] == part &&
1989  cell_color[neighbor->column()] <= n_neighbors)
1990  color_finder[cell_color[neighbor->column()]] = false;
1991  }
1992  // Choose the smallest color that is not taken for the block
1993  cell_color[cell] = 0;
1994  while (color_finder[cell_color[cell]] == false)
1995  cell_color[cell]++;
1996  if (cell_color[cell] > max_color)
1997  max_color = cell_color[cell];
1998  }
1999  // Reorder within partition: First, all blocks that belong the 0 and
2000  // then so on until those with color max (Note that the smaller the
2001  // number the larger the partition)
2002  for (unsigned int color = 0; color <= max_color; ++color)
2003  {
2004  cell_partition_data.push_back(color_counter);
2005  index_counter++;
2006  for (unsigned int k = partition_size[part];
2007  k < partition_size[part + 1];
2008  k++)
2009  {
2010  unsigned int cell = partition_list[k];
2011  if (cell_color[cell] == color)
2012  {
2013  partition_color_list[color_counter++] = cell;
2014  }
2015  }
2016  }
2017  }
2018  cell_partition_data.push_back(n_blocks);
2019  partition_row_index[partition] = index_counter;
2020  AssertDimension(color_counter, n_blocks);
2021  }
2022 
2023 
2024  // Function to create partitioning on the first layer.
2025  void
2027  const unsigned int cluster_size,
2028  std::vector<unsigned int> & cell_partition,
2029  std::vector<unsigned int> & partition_list,
2030  std::vector<unsigned int> & partition_size,
2031  unsigned int & partition) const
2032 
2033  {
2034  // For each block of cells, this variable saves to which partitions the
2035  // block belongs. Initialize all to n_cell_batches to mark them as not
2036  // yet assigned a partition.
2037  // std::vector<unsigned int> cell_partition (n_active_cells,
2038  // numbers::invalid_unsigned_int);
2039  // List of cells in previous partition
2040  std::vector<unsigned int> neighbor_list;
2041  // List of cells in current partition for use as neighbors in next
2042  // partition
2043  std::vector<unsigned int> neighbor_neighbor_list;
2044 
2045  // In element j of this variable, one puts the old number of the block
2046  // that should be the jth block in the new numeration.
2047  // std::vector<unsigned int> partition_list(n_active_cells,0);
2048 
2049  // This vector points to the start of each partition.
2050  // std::vector<unsigned int> partition_size(2,0);
2051 
2052  partition = 0;
2053  unsigned int counter = 0;
2054  unsigned int start_nonboundary =
2055  cell_partition_data.size() == 5 ?
2056  vectorization_length *
2057  (cell_partition_data[2] - cell_partition_data[1]) :
2058  0;
2059 
2060  const unsigned int n_cell_batches = *(cell_partition_data.end() - 2);
2061  if (n_cell_batches == 0)
2062  return;
2063  if (scheme == color)
2064  start_nonboundary = n_cell_batches;
2065  if (scheme == partition_color ||
2066  scheme == color) // blocking_connectivity == true
2067  start_nonboundary = ((start_nonboundary + block_size - 1) / block_size);
2068  unsigned int n_blocks;
2069  if (scheme == partition_color ||
2070  scheme == color) // blocking_connectivity == true
2071  n_blocks = this->n_blocks;
2072  else
2073  n_blocks = n_active_cells;
2074 
2075  if (start_nonboundary > n_blocks)
2076  start_nonboundary = n_blocks;
2077 
2078 
2079  unsigned int start_up = 0;
2080  bool work = true;
2081  unsigned int remainder = cluster_size;
2082 
2083  // this performs a classical breath-first search in the connectivity
2084  // graph of the cells under the restriction that the size of the
2085  // partitions should be a multiple of the given block size
2086  while (work)
2087  {
2088  // put the cells with neighbors on remote MPI processes up front
2089  if (start_nonboundary > 0)
2090  {
2091  for (unsigned int cell = 0; cell < start_nonboundary; ++cell)
2092  {
2093  const unsigned int cell_nn = cell;
2094  cell_partition[cell_nn] = partition;
2095  neighbor_list.push_back(cell_nn);
2096  partition_list[counter++] = cell_nn;
2097  partition_size.back()++;
2098  }
2099  start_nonboundary = 0;
2100  remainder -= (start_nonboundary % cluster_size);
2101  if (remainder == cluster_size)
2102  remainder = 0;
2103  }
2104  else
2105  {
2106  // To start up, set the start_up cell to partition and list all
2107  // its neighbors.
2108  cell_partition[start_up] = partition;
2109  neighbor_list.push_back(start_up);
2110  partition_list[counter++] = start_up;
2111  partition_size.back()++;
2112  start_up++;
2113  remainder--;
2114  if (remainder == cluster_size)
2115  remainder = 0;
2116  }
2117  int index_before = neighbor_list.size(), index = index_before,
2118  index_stop = 0;
2119  while (remainder > 0)
2120  {
2121  if (index == index_stop)
2122  {
2123  index = neighbor_list.size();
2124  if (index == index_before)
2125  {
2126  neighbor_list.resize(0);
2127  goto not_connect;
2128  }
2129  index_stop = index_before;
2130  index_before = index;
2131  }
2132  index--;
2133  unsigned int additional = neighbor_list[index];
2135  connectivity.begin(additional),
2136  end =
2137  connectivity.end(additional);
2138  for (; neighbor != end; ++neighbor)
2139  {
2140  if (cell_partition[neighbor->column()] ==
2142  {
2143  partition_size.back()++;
2144  cell_partition[neighbor->column()] = partition;
2145  neighbor_list.push_back(neighbor->column());
2146  partition_list[counter++] = neighbor->column();
2147  remainder--;
2148  if (remainder == 0)
2149  break;
2150  }
2151  }
2152  }
2153 
2154  while (neighbor_list.size() > 0)
2155  {
2156  partition++;
2157 
2158  // counter for number of cells so far in current partition
2159  unsigned int partition_counter = 0;
2160 
2161  // Mark the start of the new partition
2162  partition_size.push_back(partition_size.back());
2163 
2164  // Loop through the list of cells in previous partition and put
2165  // all their neighbors in current partition
2166  for (const unsigned int cell : neighbor_list)
2167  {
2168  Assert(cell_partition[cell] == partition - 1,
2169  ExcInternalError());
2170  auto neighbor = connectivity.begin(cell);
2171  const auto end = connectivity.end(cell);
2172  for (; neighbor != end; ++neighbor)
2173  {
2174  if (cell_partition[neighbor->column()] ==
2176  {
2177  partition_size.back()++;
2178  cell_partition[neighbor->column()] = partition;
2179 
2180  // collect the cells of the current partition for
2181  // use as neighbors in next partition
2182  neighbor_neighbor_list.push_back(neighbor->column());
2183  partition_list[counter++] = neighbor->column();
2184  partition_counter++;
2185  }
2186  }
2187  }
2188  remainder = cluster_size - (partition_counter % cluster_size);
2189  if (remainder == cluster_size)
2190  remainder = 0;
2191  int index_stop = 0;
2192  int index_before = neighbor_neighbor_list.size(),
2193  index = index_before;
2194  while (remainder > 0)
2195  {
2196  if (index == index_stop)
2197  {
2198  index = neighbor_neighbor_list.size();
2199  if (index == index_before)
2200  {
2201  neighbor_neighbor_list.resize(0);
2202  break;
2203  }
2204  index_stop = index_before;
2205  index_before = index;
2206  }
2207  index--;
2208  unsigned int additional = neighbor_neighbor_list[index];
2210  connectivity.begin(
2211  additional),
2212  end = connectivity.end(
2213  additional);
2214  for (; neighbor != end; ++neighbor)
2215  {
2216  if (cell_partition[neighbor->column()] ==
2218  {
2219  partition_size.back()++;
2220  cell_partition[neighbor->column()] = partition;
2221  neighbor_neighbor_list.push_back(neighbor->column());
2222  partition_list[counter++] = neighbor->column();
2223  remainder--;
2224  if (remainder == 0)
2225  break;
2226  }
2227  }
2228  }
2229 
2230  neighbor_list = neighbor_neighbor_list;
2231  neighbor_neighbor_list.resize(0);
2232  }
2233  not_connect:
2234  // One has to check if the graph is not connected so we have to find
2235  // another partition.
2236  work = false;
2237  for (unsigned int j = start_up; j < n_blocks; ++j)
2238  if (cell_partition[j] == numbers::invalid_unsigned_int)
2239  {
2240  start_up = j;
2241  work = true;
2242  if (remainder == 0)
2243  remainder = cluster_size;
2244  break;
2245  }
2246  }
2247  if (remainder != 0)
2248  partition++;
2249 
2250  AssertDimension(partition_size[partition], n_blocks);
2251  }
2252 
2253 
2254  void
2256  {
2257  evens = (partition + 1) / 2;
2258  odds = partition / 2;
2259  n_blocked_workers = odds - (odds + evens + 1) % 2;
2260  n_workers = evens + odds - n_blocked_workers;
2261  // From here only used for partition partition option.
2262  partition_evens.resize(partition);
2263  partition_odds.resize(partition);
2264  partition_n_blocked_workers.resize(partition);
2265  partition_n_workers.resize(partition);
2266  for (unsigned int part = 0; part < partition; ++part)
2267  {
2268  partition_evens[part] =
2269  (partition_row_index[part + 1] - partition_row_index[part] + 1) / 2;
2270  partition_odds[part] =
2271  (partition_row_index[part + 1] - partition_row_index[part]) / 2;
2272  partition_n_blocked_workers[part] =
2273  partition_odds[part] -
2274  (partition_odds[part] + partition_evens[part] + 1) % 2;
2275  partition_n_workers[part] = partition_evens[part] +
2276  partition_odds[part] -
2277  partition_n_blocked_workers[part];
2278  }
2279  }
2280  } // namespace MatrixFreeFunctions
2281 } // namespace internal
2282 
2283 
2284 
2285 // explicit instantiations of template functions
2286 template void
2287 internal::MatrixFreeFunctions::TaskInfo::print_memory_statistics<std::ostream>(
2288  std::ostream &,
2289  const std::size_t) const;
2290 template void
2292  ConditionalOStream>(ConditionalOStream &, const std::size_t) const;
2293 
2294 
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1112
MPICommunication(MFWorkerInterface &worker_in, const bool do_compress)
Definition: task_info.cc:323
void operator()(const tbb::blocked_range< unsigned int > &r) const
Definition: task_info.cc:254
static const unsigned int invalid_unsigned_int
Definition: types.h:196
virtual void cell_loop_post_range(const unsigned int range_index)=0
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
PartitionWork(MFWorkerInterface &worker_in, const unsigned int partition_in, const TaskInfo &task_info_in, const bool is_blocked_in)
Definition: task_info.cc:281
virtual void cell_loop_pre_range(const unsigned int range_index)=0
virtual void cell(const std::pair< unsigned int, unsigned int > &cell_range)=0
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void parallel_for(Iterator x_begin, Iterator x_end, const Functor &functor, const unsigned int grainsize)
Definition: parallel.h:80
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
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:348
void add(const size_type i, const size_type j)
CellWork(MFWorkerInterface &worker_in, const TaskInfo &task_info_in, const unsigned int partition_in)
Definition: task_info.cc:245
void guess_block_size(const unsigned int dofs_per_cell)
Definition: task_info.cc:1147
std::vector< unsigned int > partition_row_index
Definition: task_info.h:466
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
constexpr int block_size
Definition: cuda_size.h:29
std::vector< unsigned int > partition_n_workers
Definition: task_info.h:576
static const char U
virtual void vector_update_ghosts_finish()=0
Finishes the communication for the update ghost values operation.
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:474
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
void make_connectivity_cells_to_blocks(const std::vector< unsigned char > &irregular_cells, const DynamicSparsityPattern &connectivity_cells, DynamicSparsityPattern &connectivity_blocks) const
Definition: task_info.cc:1585
std::vector< unsigned int > partition_n_blocked_workers
Definition: task_info.h:570
void make_partitioning_within_partitions_post_blocked(const DynamicSparsityPattern &connectivity, const std::vector< unsigned int > &cell_active_fe_index, const unsigned int partition, const unsigned int cluster_size, const bool hp_bool, const std::vector< unsigned int > &cell_partition, const std::vector< unsigned int > &partition_list, const std::vector< unsigned int > &partition_size, std::vector< unsigned int > &partition_partition_list, std::vector< unsigned char > &irregular_cells)
Definition: task_info.cc:1629
CellWork(MFWorkerInterface &worker, const unsigned int partition, const TaskInfo &task_info, const bool is_blocked)
Definition: task_info.cc:117
void make_coloring_within_partitions_pre_blocked(const DynamicSparsityPattern &connectivity, const unsigned int partition, const std::vector< unsigned int > &cell_partition, const std::vector< unsigned int > &partition_list, const std::vector< unsigned int > &partition_size, std::vector< unsigned int > &partition_color_list)
Definition: task_info.cc:1950
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
#define Assert(cond, exc)
Definition: exceptions.h:1461
void make_boundary_cells_divisible(std::vector< unsigned int > &boundary_cells)
Definition: task_info.cc:720
void make_thread_graph_partition_color(DynamicSparsityPattern &connectivity, std::vector< unsigned int > &renumbering, std::vector< unsigned char > &irregular_cells, const bool hp_bool)
Definition: task_info.cc:1175
void update_task_info(const unsigned int partition)
Definition: task_info.cc:2255
void initial_setup_blocks_tasks(const std::vector< unsigned int > &boundary_cells, std::vector< unsigned int > &renumbering, std::vector< unsigned char > &incompletely_filled_vectorization)
Definition: task_info.cc:1081
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:407
virtual void face(const unsigned int range_index)=0
virtual void zero_dst_vector_range(const unsigned int range_index)=0
VectorType::value_type * end(VectorType &V)
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1487
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13725
ActualCellWork(MFWorkerInterface &worker, const unsigned int partition, const TaskInfo &task_info)
Definition: task_info.cc:83
virtual void vector_update_ghosts_start()=0
Starts the communication for the update ghost values operation.
Point< 2 > first
Definition: grid_out.cc:4614
std::vector< unsigned int > face_partition_data
Definition: task_info.h:496
virtual void boundary(const unsigned int range_index)=0
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
size_type row_length(const size_type row) const
void make_thread_graph_partition_partition(const std::vector< unsigned int > &cell_active_fe_index, DynamicSparsityPattern &connectivity, std::vector< unsigned int > &renumbering, std::vector< unsigned char > &irregular_cells, const bool hp_bool)
Definition: task_info.cc:1517
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:406
void make_partitioning(const DynamicSparsityPattern &connectivity, const unsigned int cluster_size, std::vector< unsigned int > &cell_partition, std::vector< unsigned int > &partition_list, std::vector< unsigned int > &partition_size, unsigned int &partition) const
Definition: task_info.cc:2026
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
virtual void vector_compress_start()=0
Starts the communication for the vector compress operation.
std::vector< unsigned int > partition_odds
Definition: task_info.h:564
std::size_t memory_consumption() const
Definition: task_info.cc:703
static ::ExceptionBase & ExcNotImplemented()
PartitionWork(MFWorkerInterface &function_in, const unsigned int partition_in, const TaskInfo &task_info_in, const bool is_blocked_in=false)
Definition: task_info.cc:148
static unsigned int n_threads()
ActualCellWork(MFWorkerInterface **worker_pointer, const unsigned int partition, const TaskInfo &task_info)
Definition: task_info.cc:74
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &mpi_communicator)
Definition: mpi.cc:91
std::vector< unsigned int > partition_evens
Definition: task_info.h:558
void print_memory_statistics(StreamType &out, std::size_t data_length) const
Definition: task_info.cc:688
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()
void make_thread_graph(const std::vector< unsigned int > &cell_active_fe_index, DynamicSparsityPattern &connectivity, std::vector< unsigned int > &renumbering, std::vector< unsigned char > &irregular_cells, const bool hp_bool)
Definition: task_info.cc:1322
virtual void vector_compress_finish()=0
Finishes the communication for the vector compress operation.