27#ifdef DEAL_II_WITH_MPI
35#ifdef DEAL_II_WITH_METIS
42#ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
43# include <zoltan_cpp.h>
57 const std::vector<unsigned int> &cell_weights,
58 const unsigned int n_partitions,
59 std::vector<unsigned int> &partition_indices)
63#ifndef DEAL_II_WITH_METIS
64 (void)sparsity_pattern;
67 (void)partition_indices;
75 idx_t n =
static_cast<signed int>(sparsity_pattern.
n_rows());
90 idx_t options[METIS_NOPTIONS];
91 METIS_SetDefaultOptions(options);
95 std::vector<idx_t> int_rowstart(1);
96 int_rowstart.reserve(sparsity_pattern.
n_rows() + 1);
97 std::vector<idx_t> int_colnums;
103 col < sparsity_pattern.
end(row);
105 int_colnums.push_back(col->column());
106 int_rowstart.push_back(int_colnums.size());
109 std::vector<idx_t> int_partition_indices(sparsity_pattern.
n_rows());
112 std::vector<idx_t> int_cell_weights;
113 if (cell_weights.size() > 0)
115 Assert(cell_weights.size() == sparsity_pattern.
n_rows(),
117 sparsity_pattern.
n_rows()));
118 int_cell_weights.resize(cell_weights.size());
119 std::copy(cell_weights.begin(),
121 int_cell_weights.begin());
125 idx_t *
const p_int_cell_weights =
126 (cell_weights.size() > 0 ? int_cell_weights.data() :
nullptr);
137 ierr = METIS_PartGraphRecursive(&n,
149 int_partition_indices.data());
153 ierr = METIS_PartGraphKway(&n,
165 int_partition_indices.data());
172 std::copy(int_partition_indices.begin(),
173 int_partition_indices.end(),
174 partition_indices.begin());
180#ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
183 get_number_of_objects(
void *
data,
int *ierr)
194 get_object_list(
void *
data,
197 ZOLTAN_ID_PTR globalID,
198 ZOLTAN_ID_PTR localID,
210 auto n_dofs = graph->
n_rows();
212 for (
unsigned int i = 0; i < n_dofs; ++i)
221 get_num_edges_list(
void *
data,
225 ZOLTAN_ID_PTR globalID,
236 for (
int i = 0; i < num_obj; ++i)
239 numEdges[i] = graph->
row_length(globalID[i]) - 1;
248 get_edge_list(
void *
data,
255 ZOLTAN_ID_PTR nborGID,
264 ZOLTAN_ID_PTR nextNborGID = nborGID;
265 int *nextNborProc = nborProc;
269 i < static_cast<SparsityPattern::size_type>(num_obj);
277 if (i != col->column())
282 *nextNborGID++ = col->column();
292 const std::vector<unsigned int> &cell_weights,
293 const unsigned int n_partitions,
294 std::vector<unsigned int> &partition_indices)
298#ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
299 (void)sparsity_pattern;
302 (void)partition_indices;
307 cell_weights.empty(),
309 "The cell weighting functionality for Zoltan has not yet been implemented."));
312 std::unique_ptr<Zoltan> zz = std::make_unique<Zoltan>(MPI_COMM_SELF);
316 zz->Set_Param(
"DEBUG_LEVEL",
"0");
320 zz->Set_Param(
"NUM_LOCAL_PARTS",
321 std::to_string(n_partitions));
336 zz->Set_Param(
"PHG_EDGE_SIZE_THRESHOLD",
"0.5");
343 zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
344 zz->Set_Obj_List_Fn(get_object_list, &graph);
345 zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
346 zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
350 int num_gid_entries = 1;
351 int num_lid_entries = 1;
353 ZOLTAN_ID_PTR import_global_ids =
nullptr;
354 ZOLTAN_ID_PTR import_local_ids =
nullptr;
355 int *import_procs =
nullptr;
356 int *import_to_part =
nullptr;
358 ZOLTAN_ID_PTR export_global_ids =
nullptr;
359 ZOLTAN_ID_PTR export_local_ids =
nullptr;
360 int *export_procs =
nullptr;
361 int *export_to_part =
nullptr;
364 const int rc = zz->LB_Partition(changes,
384 std::fill(partition_indices.begin(), partition_indices.end(), 0);
388 for (
int i = 0; i < num_export; ++i)
389 partition_indices[export_local_ids[i]] = export_to_part[i];
397 const unsigned int n_partitions,
398 std::vector<unsigned int> &partition_indices,
401 std::vector<unsigned int> cell_weights;
414 const std::vector<unsigned int> &cell_weights,
415 const unsigned int n_partitions,
416 std::vector<unsigned int> &partition_indices,
425 Assert(partition_indices.size() == sparsity_pattern.
n_rows(),
427 sparsity_pattern.
n_rows()));
430 if (n_partitions == 1 || (sparsity_pattern.
n_rows() == 1))
432 std::fill_n(partition_indices.begin(), partition_indices.size(), 0U);
437 partition_metis(sparsity_pattern,
442 partition_zoltan(sparsity_pattern,
453 std::vector<unsigned int> &color_indices)
457#ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
458 (void)sparsity_pattern;
464 std::unique_ptr<Zoltan> zz = std::make_unique<Zoltan>(MPI_COMM_SELF);
468 zz->Set_Param(
"DEBUG_LEVEL",
"0");
469 zz->Set_Param(
"COLORING_PROBLEM",
"DISTANCE-1");
470 zz->Set_Param(
"NUM_GID_ENTRIES",
"1");
471 zz->Set_Param(
"NUM_LID_ENTRIES",
"1");
472 zz->Set_Param(
"OBJ_WEIGHT_DIM",
"0");
473 zz->Set_Param(
"RECOLORING_NUM_OF_ITERATIONS",
"0");
480 zz->Set_Num_Obj_Fn(get_number_of_objects, &graph);
481 zz->Set_Obj_List_Fn(get_object_list, &graph);
482 zz->Set_Num_Edges_Multi_Fn(get_num_edges_list, &graph);
483 zz->Set_Edge_List_Multi_Fn(get_edge_list, &graph);
486 int num_gid_entries = 1;
487 const int num_objects = graph.
n_rows();
490 std::vector<ZOLTAN_ID_TYPE> global_ids(num_objects);
491 std::vector<int> color_exp(num_objects);
494 for (
int i = 0; i < num_objects; ++i)
498 int rc = zz->Color(num_gid_entries,
506 color_indices.resize(num_objects);
507 Assert(color_exp.size() == color_indices.size(),
510 std::copy(color_exp.begin(), color_exp.end(), color_indices.begin());
512 unsigned int n_colors =
513 *(std::max_element(color_indices.begin(), color_indices.end()));
529 const std::vector<DynamicSparsityPattern::size_type> &new_indices)
539 if (sparsity.
row_length(row) < min_coordination)
542 starting_point = row;
567 return starting_point;
576 std::vector<DynamicSparsityPattern::size_type> &new_indices,
577 const std::vector<DynamicSparsityPattern::size_type> &starting_indices)
585 "You can't specify more starting indices than there are rows"));
589 "Only valid for sparsity patterns which store all rows."));
590 for (
const auto starting_index : starting_indices)
592 (void)starting_index;
594 ExcMessage(
"Invalid starting index: All starting indices need "
595 "to be between zero and the number of rows in the "
596 "sparsity pattern."));
601 std::vector<DynamicSparsityPattern::size_type> last_round_dofs(
605 std::fill(new_indices.begin(),
611 if (last_round_dofs.empty())
612 last_round_dofs.push_back(
619 for (
const auto &last_round_dof : last_round_dofs)
620 new_indices[last_round_dof] = next_free_number++;
623 std::vector<DynamicSparsityPattern::size_type> next_round_dofs;
627 std::vector<std::pair<unsigned int, DynamicSparsityPattern::size_type>>
628 dofs_by_coordination;
633 next_round_dofs.clear();
636 for (
const auto dof : last_round_dofs)
638 const unsigned int row_length = sparsity.
row_length(dof);
639 for (
unsigned int i = 0; i < row_length; ++i)
645 next_round_dofs.push_back(column);
650 new_indices[column] = 0;
659 if (next_round_dofs.empty())
661 if (std::find(new_indices.begin(),
672 Assert(starting_indices.empty(),
673 ExcMessage(
"The input graph appears to have more than one "
674 "component, but as stated in the documentation "
675 "we only want to reorder such graphs if no "
676 "starting indices are given. The function was "
677 "called with starting indices, however."));
679 next_round_dofs.push_back(
685 dofs_by_coordination.clear();
687 dofs_by_coordination.emplace_back(sparsity.
row_length(next_round_dof),
689 std::sort(dofs_by_coordination.begin(), dofs_by_coordination.end());
692 for (
const auto &i : dofs_by_coordination)
693 new_indices[i.second] = next_free_number++;
696 last_round_dofs.swap(next_round_dofs);
702 Assert((std::find(new_indices.begin(),
705 (next_free_number == sparsity.
n_rows()),
716 std::vector<DynamicSparsityPattern::size_type> &renumbering)
723 "Only valid for sparsity patterns which store all rows."));
740 std::vector<types::global_dof_index> touched_nodes(n_nodes, unseen_node);
742 std::vector<unsigned int> row_lengths(n_nodes);
743 std::vector<types::global_dof_index> current_neighbors;
744 std::vector<types::global_dof_index> group_starts(1);
745 std::vector<types::global_dof_index> group_indices;
746 group_indices.reserve(n_nodes);
756 row_lengths[row] = connectivity.
row_length(row);
759 std::vector<unsigned int> n_remaining_neighbors(row_lengths);
777 if (touched_nodes[i] == unseen_node)
778 if (row_lengths[i] < candidate_valence)
781 candidate_valence = n_remaining_neighbors[i];
782 if (candidate_valence <= 1)
790 current_neighbors = {candidate_index};
791 touched_nodes[candidate_index] = available_node;
809 unsigned int candidate_row_length = 0;
810 const unsigned int loop_length = current_neighbors.size();
811 unsigned int write_index = 0;
812 for (
unsigned int i = 0; i < loop_length; ++i)
815 Assert(touched_nodes[node] != unseen_node,
817 if (touched_nodes[node] == available_node)
819 current_neighbors[write_index] = node;
821 if (n_remaining_neighbors[node] < candidate_valence ||
822 (n_remaining_neighbors[node] == candidate_valence &&
823 (row_lengths[node] > candidate_row_length ||
824 (row_lengths[node] == candidate_row_length &&
825 node < candidate_index))))
827 candidate_index = node;
828 candidate_valence = n_remaining_neighbors[node];
829 candidate_row_length = row_lengths[node];
833 current_neighbors.resize(write_index);
838 Assert(touched_nodes[node] == available_node,
843 if (current_neighbors.empty())
848 group_indices.push_back(candidate_index);
849 touched_nodes[candidate_index] = group_starts.size() - 1;
850 const auto end_it = connectivity.
end(candidate_index);
851 for (
auto it = connectivity.
begin(candidate_index); it != end_it;
853 if (touched_nodes[it->column()] >= available_node)
855 group_indices.push_back(it->column());
856 touched_nodes[it->column()] = group_starts.size() - 1;
858 group_starts.push_back(group_indices.size());
867 group_starts[group_starts.size() - 2];
868 index < group_starts.back();
871 auto it = connectivity.
begin(group_indices[index]);
872 const auto end_row = connectivity.
end(group_indices[index]);
873 for (; it != end_row; ++it)
875 if (touched_nodes[it->column()] == unseen_node)
877 current_neighbors.push_back(it->column());
878 touched_nodes[it->column()] = available_node;
880 n_remaining_neighbors[it->column()]--;
892 const unsigned int n_groups = group_starts.size() - 1;
893 if (n_groups < n_nodes)
899 index < group_starts[row + 1];
902 auto it = connectivity.
begin(group_indices[index]);
903 const auto end_it = connectivity.
end(group_indices[index]);
904 for (; it != end_it; ++it)
905 connectivity_next.
add(row, touched_nodes[it->column()]);
909 std::vector<types::global_dof_index> renumbering_next(n_groups);
916 group_starts[renumbering_next[row]];
917 index < group_starts[renumbering_next[row] + 1];
919 renumbering[c] = group_indices[index];
927 renumbering[c++] = i;
935 std::vector<DynamicSparsityPattern::size_type> &renumbering)
946#ifdef DEAL_II_WITH_MPI
952 const IndexSet &locally_relevant_rows)
955 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
958 IndexSet requested_rows(locally_relevant_rows);
961 std::vector<unsigned int> index_owner =
976 rows_data[index_owner[i]].push_back(row);
980 const auto rows_data_received =
987 for (
const auto &
data : rows_data_received)
989 for (
const auto &row :
data.second)
998 send_data[
data.first].push_back(row);
999 send_data[
data.first].push_back(rlen);
1001 send_data[
data.first].push_back(
1007 const auto received_data =
1011 for (
const auto &
data : received_data)
1013 const auto &recv_buf =
data.second;
1014 auto ptr = recv_buf.begin();
1015 const auto end = recv_buf.end();
1018 const auto row = *(ptr++);
1021 const auto n_entries = *(ptr++);
1043 const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu,
1048 std::vector<DynamicSparsityPattern::size_type> start_index(
1049 rows_per_cpu.size() + 1);
1052 start_index[i + 1] = start_index[i] + rows_per_cpu[i];
1054 IndexSet owned(start_index.back());
1055 owned.
add_range(start_index[myid], start_index[myid] + rows_per_cpu[myid]);
1064 const IndexSet &locally_owned_rows,
1066 const IndexSet &locally_relevant_rows)
1071 "The DynamicSparsityPattern must be initialized with an IndexSet that contains locally relevant indices."));
1073 IndexSet requested_rows(locally_relevant_rows);
1076 std::vector<unsigned int> index_owner =
1082 std::map<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
1084 map_vec_t send_data;
1100 send_data[index_owner[i]].push_back(row);
1101 send_data[index_owner[i]].push_back(rlen);
1106 send_data[index_owner[i]].push_back(column);
1113 for (
const auto &
data : receive_data)
1115 const auto &recv_buf =
data.second;
1116 auto ptr = recv_buf.begin();
1117 const auto end = recv_buf.end();
1120 const auto row = *(ptr++);
1122 const auto n_entries = *(ptr++);
1136 const std::vector<IndexSet> &owned_set_per_cpu,
1142 owned_set_per_cpu[myid],
1151 const IndexSet &locally_owned_rows,
1153 const IndexSet &locally_relevant_rows)
1157 std::vector<BlockDynamicSparsityPattern::size_type>>;
1158 map_vec_t send_data;
1160 IndexSet requested_rows(locally_relevant_rows);
1163 std::vector<unsigned int> index_owner =
1182 std::vector<BlockDynamicSparsityPattern::size_type> &dst =
1183 send_data[index_owner[i]];
1185 dst.push_back(rlen);
1192 dst.push_back(column);
1196 unsigned int num_receive = 0;
1198 std::vector<unsigned int> send_to;
1199 send_to.reserve(send_data.size());
1200 for (
const auto &sparsity_line : send_data)
1201 send_to.push_back(sparsity_line.first);
1208 std::vector<MPI_Request> requests(send_data.size());
1220 unsigned int idx = 0;
1221 for (
const auto &sparsity_line : send_data)
1223 const int ierr = MPI_Isend(
1224 sparsity_line.second.data(),
1225 sparsity_line.second.size(),
1226 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1227 sparsity_line.first,
1237 std::vector<BlockDynamicSparsityPattern::size_type> recv_buf;
1238 for (
unsigned int index = 0; index < num_receive; ++index)
1241 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, mpi_comm, &status);
1245 ierr = MPI_Get_count(
1247 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1251 recv_buf.resize(len);
1255 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1262 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1263 ptr = recv_buf.begin();
1264 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1265 end = recv_buf.end();
1271 for (
unsigned int c = 0; c < num; ++c)
1283 if (requests.size() > 0)
1286 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
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)
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 n_elements() const
void subtract_set(const IndexSet &other)
void add_range(const size_type begin, const size_type end)
size_type nth_index_in_set(const size_type local_index) const
types::global_dof_index size_type
bool is_compressed() const
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
unsigned int row_length(const size_type row) const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcMETISNotInstalled()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMETISError(int arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
static ::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcZOLTANNotInstalled()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
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)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
unsigned int compute_n_point_to_point_communications(const MPI_Comm mpi_comm, const std::vector< unsigned int > &destinations)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
constexpr types::global_dof_index invalid_dof_index
constexpr types::global_dof_index invalid_size_type
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)