Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sparsity_tools.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 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 
20 #include <deal.II/lac/exceptions.h>
23 
24 #include <algorithm>
25 #include <functional>
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
73  // METIS. Note that this is particularly
74  // simple, since METIS wants exactly our
75  // compressed row storage format. we only
76  // have to set up a few auxiliary arrays
77  idx_t n = static_cast<signed int>(sparsity_pattern.n_rows()),
78  ncon = 1, // number of balancing constraints (should be >0)
79  nparts =
80  static_cast<int>(n_partitions), // number of subdomains to create
81  dummy; // the numbers of edges cut by the
82  // resulting partition
83 
84  // We can not partition n items into more than n parts. METIS will
85  // generate non-sensical output (everything is owned by a single process)
86  // and complain with a message (but won't return an error code!):
87  // ***Cannot bisect a graph with 0 vertices!
88  // ***You are trying to partition a graph into too many parts!
89  nparts = std::min(n, nparts);
90 
91  // use default options for METIS
92  idx_t options[METIS_NOPTIONS];
93  METIS_SetDefaultOptions(options);
94 
95  // one more nuisance: we have to copy our own data to arrays that store
96  // signed integers :-(
97  std::vector<idx_t> int_rowstart(1);
98  int_rowstart.reserve(sparsity_pattern.n_rows() + 1);
99  std::vector<idx_t> int_colnums;
100  int_colnums.reserve(sparsity_pattern.n_nonzero_elements());
101  for (SparsityPattern::size_type row = 0; row < sparsity_pattern.n_rows();
102  ++row)
103  {
104  for (SparsityPattern::iterator col = sparsity_pattern.begin(row);
105  col < sparsity_pattern.end(row);
106  ++col)
107  int_colnums.push_back(col->column());
108  int_rowstart.push_back(int_colnums.size());
109  }
110 
111  std::vector<idx_t> int_partition_indices(sparsity_pattern.n_rows());
112 
113  // Setup cell weighting option
114  std::vector<idx_t> int_cell_weights;
115  if (cell_weights.size() > 0)
116  {
117  Assert(cell_weights.size() == sparsity_pattern.n_rows(),
118  ExcDimensionMismatch(cell_weights.size(),
119  sparsity_pattern.n_rows()));
120  int_cell_weights.resize(cell_weights.size());
121  std::copy(cell_weights.begin(),
122  cell_weights.end(),
123  int_cell_weights.begin());
124  }
125  // Set a pointer to the optional cell weighting information.
126  // METIS expects a null pointer if there are no weights to be considered.
127  idx_t *const p_int_cell_weights =
128  (cell_weights.size() > 0 ? int_cell_weights.data() : nullptr);
129 
130 
131  // Make use of METIS' error code.
132  int ierr;
133 
134  // Select which type of partitioning to create
135 
136  // Use recursive if the number of partitions is less than or equal to 8
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.size() == 0,
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 =
315  std_cxx14::make_unique<Zoltan>(MPI_COMM_SELF);
316 
317  // General parameters
318  // DEBUG_LEVEL call must precede the call to LB_METHOD
319  zz->Set_Param("DEBUG_LEVEL", "0"); // set level of debug info
320  zz->Set_Param(
321  "LB_METHOD",
322  "GRAPH"); // graph based partition method (LB-load balancing)
323  zz->Set_Param("NUM_LOCAL_PARTS",
324  std::to_string(n_partitions)); // set number of partitions
325 
326  // The PHG partitioner is a hypergraph partitioner that Zoltan could use
327  // for graph partitioning.
328  // If number of vertices in hyperedge divided by total vertices in
329  // hypergraph exceeds PHG_EDGE_SIZE_THRESHOLD,
330  // then the hyperedge will be omitted as such (dense) edges will likely
331  // incur high communication costs regardless of the partition.
332  // PHG_EDGE_SIZE_THRESHOLD value is raised to 0.5 from the default
333  // value of 0.25 so that the PHG partitioner doesn't throw warning saying
334  // "PHG_EDGE_SIZE_THRESHOLD is low ..." after removing all dense edges.
335  // For instance, in two dimensions if the triangulation being partitioned
336  // is two quadrilaterals sharing an edge and if PHG_EDGE_SIZE_THRESHOLD
337  // value is set to 0.25, PHG will remove all the edges throwing the
338  // above warning.
339  zz->Set_Param("PHG_EDGE_SIZE_THRESHOLD", "0.5");
340 
341  // Need a non-const object equal to sparsity_pattern
342  SparsityPattern graph;
343  graph.copy_from(sparsity_pattern);
344 
345  // Set query functions
346  zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
347  zz->Set_Obj_List_Fn(get_object_list, &graph);
348  zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
349  zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
350 
351  // Variables needed by partition function
352  int changes = 0;
353  int num_gid_entries = 1;
354  int num_lid_entries = 1;
355  int num_import = 0;
356  ZOLTAN_ID_PTR import_global_ids = nullptr;
357  ZOLTAN_ID_PTR import_local_ids = nullptr;
358  int * import_procs = nullptr;
359  int * import_to_part = nullptr;
360  int num_export = 0;
361  ZOLTAN_ID_PTR export_global_ids = nullptr;
362  ZOLTAN_ID_PTR export_local_ids = nullptr;
363  int * export_procs = nullptr;
364  int * export_to_part = nullptr;
365 
366  // call partitioner
367  const int rc = zz->LB_Partition(changes,
368  num_gid_entries,
369  num_lid_entries,
370  num_import,
371  import_global_ids,
372  import_local_ids,
373  import_procs,
374  import_to_part,
375  num_export,
376  export_global_ids,
377  export_local_ids,
378  export_procs,
379  export_to_part);
380  (void)rc;
381 
382  // check for error code in partitioner
383  Assert(rc == ZOLTAN_OK, ExcInternalError());
384 
385  // By default, all indices belong to part 0. After zoltan partition
386  // some are migrated to different part ID, which is stored in
387  // export_to_part array.
388  std::fill(partition_indices.begin(), partition_indices.end(), 0);
389 
390  // copy from export_to_part to partition_indices, whose part_ids != 0.
391  Assert(export_to_part != nullptr, ExcInternalError());
392  for (int i = 0; i < num_export; i++)
393  partition_indices[export_local_ids[i]] = export_to_part[i];
394 #endif
395  }
396  } // namespace
397 
398 
399  void
400  partition(const SparsityPattern & sparsity_pattern,
401  const unsigned int n_partitions,
402  std::vector<unsigned int> &partition_indices,
403  const Partitioner partitioner)
404  {
405  std::vector<unsigned int> cell_weights;
406 
407  // Call the other more general function
408  partition(sparsity_pattern,
409  cell_weights,
410  n_partitions,
411  partition_indices,
412  partitioner);
413  }
414 
415 
416  void
417  partition(const SparsityPattern & sparsity_pattern,
418  const std::vector<unsigned int> &cell_weights,
419  const unsigned int n_partitions,
420  std::vector<unsigned int> & partition_indices,
421  const Partitioner partitioner)
422  {
423  Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
424  ExcNotQuadratic());
425  Assert(sparsity_pattern.is_compressed(),
427 
428  Assert(n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions));
429  Assert(partition_indices.size() == sparsity_pattern.n_rows(),
430  ExcInvalidArraySize(partition_indices.size(),
431  sparsity_pattern.n_rows()));
432 
433  // check for an easy return
434  if (n_partitions == 1 || (sparsity_pattern.n_rows() == 1))
435  {
436  std::fill_n(partition_indices.begin(), partition_indices.size(), 0U);
437  return;
438  }
439 
440  if (partitioner == Partitioner::metis)
441  partition_metis(sparsity_pattern,
442  cell_weights,
443  n_partitions,
444  partition_indices);
445  else if (partitioner == Partitioner::zoltan)
446  partition_zoltan(sparsity_pattern,
447  cell_weights,
448  n_partitions,
449  partition_indices);
450  else
451  AssertThrow(false, ExcInternalError());
452  }
453 
454 
455  unsigned int
456  color_sparsity_pattern(const SparsityPattern & sparsity_pattern,
457  std::vector<unsigned int> &color_indices)
458  {
459  // Make sure that ZOLTAN is actually
460  // installed and detected
461 #ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
462  (void)sparsity_pattern;
463  (void)color_indices;
465  return 0;
466 #else
467  // coloring algorithm is run in serial by each processor.
468  std::unique_ptr<Zoltan> zz = std_cxx14::make_unique<Zoltan>(MPI_COMM_SELF);
469 
470  // Coloring parameters
471  // DEBUG_LEVEL must precede all other calls
472  zz->Set_Param("DEBUG_LEVEL", "0"); // level of debug info
473  zz->Set_Param("COLORING_PROBLEM", "DISTANCE-1"); // Standard coloring
474  zz->Set_Param("NUM_GID_ENTRIES", "1"); // 1 entry represents global ID
475  zz->Set_Param("NUM_LID_ENTRIES", "1"); // 1 entry represents local ID
476  zz->Set_Param("OBJ_WEIGHT_DIM", "0"); // object weights not used
477  zz->Set_Param("RECOLORING_NUM_OF_ITERATIONS", "0");
478 
479  // Zoltan::Color function requires a non-const SparsityPattern object
480  SparsityPattern graph;
481  graph.copy_from(sparsity_pattern);
482 
483  // Set query functions required by coloring function
484  zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
485  zz->Set_Obj_List_Fn(get_object_list, &graph);
486  zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
487  zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
488 
489  // Variables needed by coloring function
490  int num_gid_entries = 1;
491  const int num_objects = graph.n_rows();
492 
493  // Preallocate input variables. Element type fixed by ZOLTAN.
494  std::vector<ZOLTAN_ID_TYPE> global_ids(num_objects);
495  std::vector<int> color_exp(num_objects);
496 
497  // Set ids for which coloring needs to be done
498  for (int i = 0; i < num_objects; i++)
499  global_ids[i] = i;
500 
501  // Call ZOLTAN coloring algorithm
502  int rc = zz->Color(num_gid_entries,
503  num_objects,
504  global_ids.data(),
505  color_exp.data());
506 
507  (void)rc;
508  // Check for error code
509  Assert(rc == ZOLTAN_OK, ExcInternalError());
510 
511  // Allocate and assign color indices
512  color_indices.resize(num_objects);
513  Assert(color_exp.size() == color_indices.size(),
514  ExcDimensionMismatch(color_exp.size(), color_indices.size()));
515 
516  std::copy(color_exp.begin(), color_exp.end(), color_indices.begin());
517 
518  unsigned int n_colors =
519  *(std::max_element(color_indices.begin(), color_indices.end()));
520  return n_colors;
521 #endif
522  }
523 
524 
525  namespace internal
526  {
534  const DynamicSparsityPattern & sparsity,
535  const std::vector<DynamicSparsityPattern::size_type> &new_indices)
536  {
537  DynamicSparsityPattern::size_type starting_point =
539  DynamicSparsityPattern::size_type min_coordination = sparsity.n_rows();
540  for (DynamicSparsityPattern::size_type row = 0; row < sparsity.n_rows();
541  ++row)
542  // look over all as-yet unnumbered indices
543  if (new_indices[row] == numbers::invalid_size_type)
544  {
545  if (sparsity.row_length(row) < min_coordination)
546  {
547  min_coordination = sparsity.row_length(row);
548  starting_point = row;
549  }
550  }
551 
552  // now we still have to care for the case that no unnumbered dof has a
553  // coordination number less than sparsity.n_rows(). this rather exotic
554  // case only happens if we only have one cell, as far as I can see,
555  // but there may be others as well.
556  //
557  // if that should be the case, we can chose an arbitrary dof as
558  // starting point, e.g. the first unnumbered one
559  if (starting_point == numbers::invalid_size_type)
560  {
561  for (DynamicSparsityPattern::size_type i = 0; i < new_indices.size();
562  ++i)
563  if (new_indices[i] == numbers::invalid_size_type)
564  {
565  starting_point = i;
566  break;
567  }
568 
569  Assert(starting_point != numbers::invalid_size_type,
570  ExcInternalError());
571  }
572 
573  return starting_point;
574  }
575  } // namespace internal
576 
577 
578 
579  void
581  const DynamicSparsityPattern & sparsity,
582  std::vector<DynamicSparsityPattern::size_type> & new_indices,
583  const std::vector<DynamicSparsityPattern::size_type> &starting_indices)
584  {
585  Assert(sparsity.n_rows() == sparsity.n_cols(),
586  ExcDimensionMismatch(sparsity.n_rows(), sparsity.n_cols()));
587  Assert(sparsity.n_rows() == new_indices.size(),
588  ExcDimensionMismatch(sparsity.n_rows(), new_indices.size()));
589  Assert(starting_indices.size() <= sparsity.n_rows(),
590  ExcMessage(
591  "You can't specify more starting indices than there are rows"));
592  Assert(sparsity.row_index_set().size() == 0 ||
593  sparsity.row_index_set().size() == sparsity.n_rows(),
594  ExcMessage(
595  "Only valid for sparsity patterns which store all rows."));
596  for (const auto starting_index : starting_indices)
597  {
598  (void)starting_index;
599  Assert(starting_index < sparsity.n_rows(),
600  ExcMessage("Invalid starting index: All starting indices need "
601  "to be between zero and the number of rows in the "
602  "sparsity pattern."));
603  }
604 
605  // store the indices of the dofs renumbered in the last round. Default to
606  // starting points
607  std::vector<DynamicSparsityPattern::size_type> last_round_dofs(
608  starting_indices);
609 
610  // initialize the new_indices array with invalid values
611  std::fill(new_indices.begin(),
612  new_indices.end(),
614 
615  // if no starting indices were given: find dof with lowest coordination
616  // number
617  if (last_round_dofs.empty())
618  last_round_dofs.push_back(
619  internal::find_unnumbered_starting_index(sparsity, new_indices));
620 
621  // store next free dof index
622  DynamicSparsityPattern::size_type next_free_number = 0;
623 
624  // enumerate the first round dofs
625  for (DynamicSparsityPattern::size_type i = 0; i != last_round_dofs.size();
626  ++i)
627  new_indices[last_round_dofs[i]] = next_free_number++;
628 
629  // now do as many steps as needed to renumber all dofs
630  while (true)
631  {
632  // store the indices of the dofs to be renumbered in the next round
633  std::vector<DynamicSparsityPattern::size_type> next_round_dofs;
634 
635  // find all neighbors of the dofs numbered in the last round
636  for (const auto dof : last_round_dofs)
637  for (DynamicSparsityPattern::iterator j = sparsity.begin(dof);
638  j < sparsity.end(dof);
639  ++j)
640  next_round_dofs.push_back(j->column());
641 
642  // sort dof numbers
643  std::sort(next_round_dofs.begin(), next_round_dofs.end());
644 
645  // delete multiple entries
646  std::vector<DynamicSparsityPattern::size_type>::iterator end_sorted;
647  end_sorted =
648  std::unique(next_round_dofs.begin(), next_round_dofs.end());
649  next_round_dofs.erase(end_sorted, next_round_dofs.end());
650 
651  // eliminate dofs which are already numbered
652  for (int s = next_round_dofs.size() - 1; s >= 0; --s)
653  if (new_indices[next_round_dofs[s]] != numbers::invalid_size_type)
654  next_round_dofs.erase(next_round_dofs.begin() + s);
655 
656  // check whether there are any new dofs in the list. if there are
657  // none, then we have completely numbered the current component of the
658  // graph. check if there are as yet unnumbered components of the graph
659  // that we would then have to do next
660  if (next_round_dofs.empty())
661  {
662  if (std::find(new_indices.begin(),
663  new_indices.end(),
664  numbers::invalid_size_type) == new_indices.end())
665  // no unnumbered indices, so we can leave now
666  break;
667 
668  // otherwise find a valid starting point for the next component of
669  // the graph and continue with numbering that one. we only do so
670  // if no starting indices were provided by the user (see the
671  // documentation of this function) so produce an error if we got
672  // here and starting indices were given
673  Assert(starting_indices.empty(),
674  ExcMessage("The input graph appears to have more than one "
675  "component, but as stated in the documentation "
676  "we only want to reorder such graphs if no "
677  "starting indices are given. The function was "
678  "called with starting indices, however."))
679 
680  next_round_dofs.push_back(
682  new_indices));
683  }
684 
685 
686 
687  // store for each coordination number the dofs with these coordination
688  // number
689  std::multimap<DynamicSparsityPattern::size_type, int>
690  dofs_by_coordination;
691 
692  // find coordination number for each of these dofs
693  for (const types::global_dof_index next_round_dof : next_round_dofs)
694  {
695  const DynamicSparsityPattern::size_type coordination =
696  sparsity.row_length(next_round_dof);
697 
698  // insert this dof at its coordination number
699  const std::pair<const DynamicSparsityPattern::size_type, int>
700  new_entry(coordination, next_round_dof);
701  dofs_by_coordination.insert(new_entry);
702  }
703 
704  // assign new DoF numbers to the elements of the present front:
705  std::multimap<DynamicSparsityPattern::size_type, int>::iterator i;
706  for (i = dofs_by_coordination.begin(); i != dofs_by_coordination.end();
707  ++i)
708  new_indices[i->second] = next_free_number++;
709 
710  // after that: copy this round's dofs for the next round
711  last_round_dofs = next_round_dofs;
712  }
713 
714  // test for all indices numbered. this mostly tests whether the
715  // front-marching-algorithm (which Cuthill-McKee actually is) has reached
716  // all points.
717  Assert((std::find(new_indices.begin(),
718  new_indices.end(),
719  numbers::invalid_size_type) == new_indices.end()) &&
720  (next_free_number == sparsity.n_rows()),
721  ExcInternalError());
722  }
723 
724 
725 
726  namespace internal
727  {
728  void
730  const DynamicSparsityPattern & connectivity,
731  std::vector<DynamicSparsityPattern::size_type> &renumbering)
732  {
733  AssertDimension(connectivity.n_rows(), connectivity.n_cols());
734  AssertDimension(connectivity.n_rows(), renumbering.size());
735  Assert(connectivity.row_index_set().size() == 0 ||
736  connectivity.row_index_set().size() == connectivity.n_rows(),
737  ExcMessage(
738  "Only valid for sparsity patterns which store all rows."));
739 
740  std::vector<types::global_dof_index> touched_nodes(
741  connectivity.n_rows(), numbers::invalid_dof_index);
742  std::vector<unsigned int> row_lengths(connectivity.n_rows());
743  std::set<types::global_dof_index> current_neighbors;
744  std::vector<std::vector<types::global_dof_index>> groups;
745 
746  // First collect the number of neighbors for each node. We use this
747  // field to find next nodes with the minimum number of non-touched
748  // neighbors in the field n_remaining_neighbors, so we will count down
749  // on this field. We also cache the row lengths because we need this
750  // data frequently and getting it from the sparsity pattern is more
751  // expensive.
752  for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row)
753  {
754  row_lengths[row] = connectivity.row_length(row);
755  Assert(row_lengths[row] > 0, ExcInternalError());
756  }
757  std::vector<unsigned int> n_remaining_neighbors(row_lengths);
758 
759  // This outer loop is typically traversed only once, unless the global
760  // graph is not connected
761  while (true)
762  {
763  // Find cell with the minimal number of neighbors (typically a
764  // corner node when based on FEM meshes). If no cell is left, we are
765  // done. Together with the outer while loop, this loop can possibly
766  // be of quadratic complexity in the number of disconnected
767  // partitions, i.e. up to connectivity.n_rows() in the worst case,
768  // but that is not the usual use case of this loop and thus not
769  // optimized for.
770  std::pair<types::global_dof_index, types::global_dof_index>
771  min_neighbors(numbers::invalid_dof_index,
773  for (types::global_dof_index i = 0; i < touched_nodes.size(); ++i)
774  if (touched_nodes[i] == numbers::invalid_dof_index)
775  if (row_lengths[i] < min_neighbors.second)
776  {
777  min_neighbors = std::make_pair(i, n_remaining_neighbors[i]);
778  if (n_remaining_neighbors[i] <= 1)
779  break;
780  }
781  if (min_neighbors.first == numbers::invalid_dof_index)
782  break;
783 
784  Assert(min_neighbors.second > 0, ExcInternalError());
785 
786  current_neighbors.clear();
787  current_neighbors.insert(min_neighbors.first);
788  while (!current_neighbors.empty())
789  {
790  // Find node with minimum number of untouched neighbors among the
791  // next set of possible neighbors
792  min_neighbors = std::make_pair(numbers::invalid_dof_index,
794  for (const auto current_neighbor : current_neighbors)
795  {
796  Assert(touched_nodes[current_neighbor] ==
798  ExcInternalError());
799  if (n_remaining_neighbors[current_neighbor] <
800  min_neighbors.second)
801  min_neighbors =
802  std::make_pair(current_neighbor,
803  n_remaining_neighbors[current_neighbor]);
804  }
805 
806  // Among the set of nodes with the minimal number of neighbors,
807  // choose the one with the largest number of touched neighbors,
808  // i.e., the one with the largest row length
809  const types::global_dof_index best_row_length =
810  min_neighbors.second;
811  for (const auto current_neighbor : current_neighbors)
812  if (n_remaining_neighbors[current_neighbor] == best_row_length)
813  if (row_lengths[current_neighbor] > min_neighbors.second)
814  min_neighbors =
815  std::make_pair(current_neighbor,
816  row_lengths[current_neighbor]);
817 
818  // Add the pivot and all direct neighbors of the pivot node not
819  // yet touched to the list of new entries.
820  groups.emplace_back();
821  std::vector<types::global_dof_index> &next_group = groups.back();
822 
823  next_group.push_back(min_neighbors.first);
824  touched_nodes[min_neighbors.first] = groups.size() - 1;
826  connectivity.begin(min_neighbors.first);
827  it != connectivity.end(min_neighbors.first);
828  ++it)
829  if (touched_nodes[it->column()] == numbers::invalid_dof_index)
830  {
831  next_group.push_back(it->column());
832  touched_nodes[it->column()] = groups.size() - 1;
833  }
834 
835  // Add all neighbors of the current list not yet touched to the
836  // set of possible next pivots. The added node is no longer a
837  // valid neighbor (here we assume symmetry of the
838  // connectivity). Delete the entries of the current list from
839  // the set of possible next pivots.
840  for (const auto index : next_group)
841  {
843  connectivity.begin(index);
844  it != connectivity.end(index);
845  ++it)
846  {
847  if (touched_nodes[it->column()] ==
849  current_neighbors.insert(it->column());
850  n_remaining_neighbors[it->column()]--;
851  }
852  current_neighbors.erase(index);
853  }
854  }
855  }
856 
857  // Sanity check: for all nodes, there should not be any neighbors left
858  for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row)
859  Assert(n_remaining_neighbors[row] == 0, ExcInternalError());
860 
861  // If the number of groups is smaller than the number of nodes, we
862  // continue by recursively calling this method
863  if (groups.size() < connectivity.n_rows())
864  {
865  // Form the connectivity of the groups
866  DynamicSparsityPattern connectivity_next(groups.size(),
867  groups.size());
868  for (types::global_dof_index i = 0; i < groups.size(); ++i)
869  for (types::global_dof_index col = 0; col < groups[i].size(); ++col)
871  connectivity.begin(groups[i][col]);
872  it != connectivity.end(groups[i][col]);
873  ++it)
874  connectivity_next.add(i, touched_nodes[it->column()]);
875 
876  // Recursively call the reordering
877  std::vector<types::global_dof_index> renumbering_next(groups.size());
878  reorder_hierarchical(connectivity_next, renumbering_next);
879 
880  // Renumber the indices group by group according to the incoming
881  // ordering for the groups
882  for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i)
883  for (types::global_dof_index col = 0;
884  col < groups[renumbering_next[i]].size();
885  ++col, ++count)
886  renumbering[count] = groups[renumbering_next[i]][col];
887  }
888  else
889  {
890  // All groups should have size one and no more recursion is possible,
891  // so use the numbering of the groups
892  for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i)
893  for (types::global_dof_index col = 0; col < groups[i].size();
894  ++col, ++count)
895  renumbering[count] = groups[i][col];
896  }
897  }
898  } // namespace internal
899 
900  void
902  const DynamicSparsityPattern & connectivity,
903  std::vector<DynamicSparsityPattern::size_type> &renumbering)
904  {
905  // the internal renumbering keeps the numbering the wrong way around (but
906  // we cannot invert the numbering inside that method because it is used
907  // recursively), so invert it here
908  internal::reorder_hierarchical(connectivity, renumbering);
909  renumbering = Utilities::invert_permutation(renumbering);
910  }
911 
912 
913 
914 #ifdef DEAL_II_WITH_MPI
915 
916  void
918  const std::vector<IndexSet> &owned_rows_per_processor,
919  const MPI_Comm & mpi_comm,
920  const IndexSet & ghost_range)
921  {
922  const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
924  owned_rows_per_processor[myid],
925  mpi_comm,
926  ghost_range);
927  }
928 
929 
930 
931  void
933  const IndexSet & locally_owned_rows,
934  const MPI_Comm & mpi_comm,
935  const IndexSet & locally_relevant_rows)
936  {
937  using map_vec_t =
938  std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
939 
940  // 1. limit rows to non owned:
941  IndexSet requested_rows(locally_relevant_rows);
942  requested_rows.subtract_set(locally_owned_rows);
943 
944  std::vector<unsigned int> index_owner =
945  Utilities::MPI::compute_index_owner(locally_owned_rows,
946  requested_rows,
947  mpi_comm);
948 
949  // 2. go through requested_rows, figure out the owner and add the row to
950  // request
951  map_vec_t rows_data;
953  i < requested_rows.n_elements();
954  ++i)
955  {
957  requested_rows.nth_index_in_set(i);
958 
959  rows_data[index_owner[i]].push_back(row);
960  }
961 
962  // 3. get what others ask us to send
963  const auto rows_data_received =
964  Utilities::MPI::some_to_some(mpi_comm, rows_data);
965 
966  // 4. now prepare data to be sent in the same format as in
967  // distribute_sparsity_pattern() below, i.e.
968  // rX,num_rX,cols_rX
969  map_vec_t send_data;
970  for (const auto &data : rows_data_received)
971  {
972  for (const auto &row : data.second)
973  {
974  const auto rlen = dsp.row_length(row);
975 
976  // skip empty lines
977  if (rlen == 0)
978  continue;
979 
980  // save entries
981  send_data[data.first].push_back(row); // row index
982  send_data[data.first].push_back(rlen); // number of entries
983  for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
984  send_data[data.first].push_back(
985  dsp.column_number(row, c)); // columns
986  } // loop over rows
987  } // loop over received data
988 
989  // 5. communicate rows
990  const auto received_data =
991  Utilities::MPI::some_to_some(mpi_comm, send_data);
992 
993  // 6. add result to our sparsity
994  for (const auto &data : received_data)
995  {
996  const auto &recv_buf = data.second;
997  auto ptr = recv_buf.begin();
998  const auto end = recv_buf.end();
999  while (ptr != end)
1000  {
1001  const auto row = *(ptr++);
1002  Assert(ptr != end, ExcInternalError());
1003 
1004  const auto n_entries = *(ptr++);
1005  Assert(n_entries > 0, ExcInternalError());
1006  Assert(ptr != end, ExcInternalError());
1007 
1008  // make sure we clear whatever was previously stored
1009  // in these rows. Otherwise we can't guarantee that the
1010  // data is consistent across MPI communicator.
1011  dsp.clear_row(row);
1012 
1013  Assert(ptr + (n_entries - 1) != end, ExcInternalError());
1014  dsp.add_entries(row, ptr, ptr + n_entries, true);
1015  ptr += n_entries;
1016  }
1017  Assert(ptr == end, ExcInternalError());
1018  }
1019  }
1020 
1021 
1022 
1023  void
1025  DynamicSparsityPattern & dsp,
1026  const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu,
1027  const MPI_Comm & mpi_comm,
1028  const IndexSet & myrange)
1029  {
1030  const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1031  std::vector<DynamicSparsityPattern::size_type> start_index(
1032  rows_per_cpu.size() + 1);
1033  start_index[0] = 0;
1034  for (DynamicSparsityPattern::size_type i = 0; i < rows_per_cpu.size(); ++i)
1035  start_index[i + 1] = start_index[i] + rows_per_cpu[i];
1036 
1037  IndexSet owned(start_index.back());
1038  owned.add_range(start_index[myid], start_index[myid] + rows_per_cpu[myid]);
1039 
1040  distribute_sparsity_pattern(dsp, owned, mpi_comm, myrange);
1041  }
1042 
1043 
1044 
1045  void
1047  const IndexSet & locally_owned_rows,
1048  const MPI_Comm & mpi_comm,
1049  const IndexSet & locally_relevant_rows)
1050  {
1051  IndexSet requested_rows(locally_relevant_rows);
1052  requested_rows.subtract_set(locally_owned_rows);
1053 
1054  std::vector<unsigned int> index_owner =
1055  Utilities::MPI::compute_index_owner(locally_owned_rows,
1056  requested_rows,
1057  mpi_comm);
1058 
1059  using map_vec_t =
1060  std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
1061 
1062  map_vec_t send_data;
1063 
1065  i < requested_rows.n_elements();
1066  ++i)
1067  {
1069  requested_rows.nth_index_in_set(i);
1070 
1071  const auto rlen = dsp.row_length(row);
1072 
1073  // skip empty lines
1074  if (!rlen)
1075  continue;
1076 
1077  // save entries
1078  send_data[index_owner[i]].push_back(row); // row index
1079  send_data[index_owner[i]].push_back(rlen); // number of entries
1080  for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1081  {
1082  // columns
1083  const auto column = dsp.column_number(row, c);
1084  send_data[index_owner[i]].push_back(column);
1085  }
1086  }
1087 
1088  const auto receive_data = Utilities::MPI::some_to_some(mpi_comm, send_data);
1089 
1090  // add what we received
1091  for (const auto &data : receive_data)
1092  {
1093  const auto &recv_buf = data.second;
1094  auto ptr = recv_buf.begin();
1095  const auto end = recv_buf.end();
1096  while (ptr != end)
1097  {
1098  const auto row = *(ptr++);
1099  Assert(ptr != end, ExcInternalError());
1100  const auto n_entries = *(ptr++);
1101 
1102  Assert(ptr + (n_entries - 1) != end, ExcInternalError());
1103  dsp.add_entries(row, ptr, ptr + n_entries, true);
1104  ptr += n_entries;
1105  }
1106  Assert(ptr == end, ExcInternalError());
1107  }
1108  }
1109 
1110 
1111 
1112  void
1114  const std::vector<IndexSet> &owned_set_per_cpu,
1115  const MPI_Comm & mpi_comm,
1116  const IndexSet & myrange)
1117  {
1118  const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
1120  owned_set_per_cpu[myid],
1121  mpi_comm,
1122  myrange);
1123  }
1124 
1125 
1126 
1127  void
1129  const IndexSet & locally_owned_rows,
1130  const MPI_Comm & mpi_comm,
1131  const IndexSet &locally_relevant_rows)
1132  {
1133  using map_vec_t =
1135  std::vector<BlockDynamicSparsityPattern::size_type>>;
1136  map_vec_t send_data;
1137 
1138  IndexSet requested_rows(locally_relevant_rows);
1139  requested_rows.subtract_set(locally_owned_rows);
1140 
1141  std::vector<unsigned int> index_owner =
1142  Utilities::MPI::compute_index_owner(locally_owned_rows,
1143  requested_rows,
1144  mpi_comm);
1145 
1147  i < requested_rows.n_elements();
1148  ++i)
1149  {
1151  requested_rows.nth_index_in_set(i);
1152 
1153  BlockDynamicSparsityPattern::size_type rlen = dsp.row_length(row);
1154 
1155  // skip empty lines
1156  if (!rlen)
1157  continue;
1158 
1159  // save entries
1160  std::vector<BlockDynamicSparsityPattern::size_type> &dst =
1161  send_data[index_owner[i]];
1162 
1163  dst.push_back(rlen); // number of entries
1164  dst.push_back(row); // row index
1165  for (BlockDynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1166  {
1167  // columns
1168  BlockDynamicSparsityPattern::size_type column =
1169  dsp.column_number(row, c);
1170  dst.push_back(column);
1171  }
1172  }
1173 
1174  unsigned int num_receive = 0;
1175  {
1176  std::vector<unsigned int> send_to;
1177  send_to.reserve(send_data.size());
1178  for (const auto &sparsity_line : send_data)
1179  send_to.push_back(sparsity_line.first);
1180 
1181  num_receive =
1183  send_to);
1184  }
1185 
1186  std::vector<MPI_Request> requests(send_data.size());
1187 
1188 
1189  // send data
1190 
1191  static Utilities::MPI::CollectiveMutex mutex;
1192  Utilities::MPI::CollectiveMutex::ScopedLock lock(mutex, mpi_comm);
1193 
1194  const int mpi_tag = Utilities::MPI::internal::Tags::
1196 
1197  {
1198  unsigned int idx = 0;
1199  for (const auto &sparsity_line : send_data)
1200  {
1201  const int ierr =
1202  MPI_Isend(DEAL_II_MPI_CONST_CAST(sparsity_line.second.data()),
1203  sparsity_line.second.size(),
1205  sparsity_line.first,
1206  mpi_tag,
1207  mpi_comm,
1208  &requests[idx++]);
1209  AssertThrowMPI(ierr);
1210  }
1211  }
1212 
1213  {
1214  // receive
1215  std::vector<BlockDynamicSparsityPattern::size_type> recv_buf;
1216  for (unsigned int index = 0; index < num_receive; ++index)
1217  {
1218  MPI_Status status;
1219  int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, mpi_comm, &status);
1220  AssertThrowMPI(ierr);
1221 
1222  int len;
1223  ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
1224  AssertThrowMPI(ierr);
1225 
1226  recv_buf.resize(len);
1227  ierr = MPI_Recv(recv_buf.data(),
1228  len,
1230  status.MPI_SOURCE,
1231  status.MPI_TAG,
1232  mpi_comm,
1233  &status);
1234  AssertThrowMPI(ierr);
1235 
1236  std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1237  ptr = recv_buf.begin();
1238  std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1239  end = recv_buf.end();
1240  while (ptr != end)
1241  {
1242  BlockDynamicSparsityPattern::size_type num = *(ptr++);
1243  Assert(ptr != end, ExcInternalError());
1244  BlockDynamicSparsityPattern::size_type row = *(ptr++);
1245  for (unsigned int c = 0; c < num; ++c)
1246  {
1247  Assert(ptr != end, ExcInternalError());
1248  dsp.add(row, *ptr);
1249  ++ptr;
1250  }
1251  }
1252  Assert(ptr == end, ExcInternalError());
1253  }
1254  }
1255 
1256  // complete all sends, so that we can safely destroy the buffers.
1257  if (requests.size())
1258  {
1259  const int ierr =
1260  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1261  AssertThrowMPI(ierr);
1262  }
1263  }
1264 #endif
1265 } // namespace SparsityTools
1266 
const types::global_dof_index invalid_size_type
Definition: types.h:200
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
iterator end() const
size_type column_number(const size_type row, const unsigned int index) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
iterator begin() const
void clear_row(const size_type row)
void gather_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
size_type column_number(const size_type row, const size_type index) const
void add(const size_type i, const size_type j)
static ::ExceptionBase & ExcMETISNotInstalled()
unsigned int compute_n_point_to_point_communications(const MPI_Comm &mpi_comm, const std::vector< unsigned int > &destinations)
Definition: mpi.cc:479
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 >())
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86
static const char U
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
size_type size() const
Definition: index_set.h:1635
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
void add(const size_type i, const size_type j)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const IndexSet & row_index_set() const
SparsityTools::sparsity_tools_distribute_sparsity_pattern()
Definition: mpi_tags.h:75
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
std::string to_string(const T &t)
Definition: patterns.h:2360
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1446
std::vector< unsigned int > compute_index_owner(const IndexSet &owned_indices, const IndexSet &indices_to_look_up, const MPI_Comm &comm)
Definition: mpi.cc:1046
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
const types::global_dof_index * dummy()
Definition: dof_handler.cc:59
DynamicSparsityPattern::size_type find_unnumbered_starting_index(const DynamicSparsityPattern &sparsity, const std::vector< DynamicSparsityPattern::size_type > &new_indices)
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
void reorder_hierarchical(const DynamicSparsityPattern &connectivity, std::vector< DynamicSparsityPattern::size_type > &renumbering)
#define DEAL_II_MPI_CONST_CAST(expr)
Definition: mpi.h:76
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1674
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcZOLTANNotInstalled()
unsigned int color_sparsity_pattern(const SparsityPattern &sparsity_pattern, std::vector< unsigned int > &color_indices)
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1707
bool is_compressed() const
size_type row_length(const size_type row) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
T min(const T &t, const MPI_Comm &mpi_communicator)
unsigned int row_length(const size_type row) const
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
const types::global_dof_index invalid_dof_index
Definition: types.h:206
std::map< unsigned int, T > some_to_some(const MPI_Comm &comm, const std::map< unsigned int, T > &objects_to_send)
unsigned int row_length(const size_type row) const
size_type n_elements() const
Definition: index_set.h:1833
void copy(const T *begin, const T *end, U *dest)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1881
static ::ExceptionBase & ExcMETISError(int arg1)
size_type n_rows() const
static ::ExceptionBase & ExcInternalError()
size_type n_cols() const