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