17 #include <deal.II/base/exceptions.h> 18 #include <deal.II/base/std_cxx14/memory.h> 20 #include <deal.II/lac/exceptions.h> 21 #include <deal.II/lac/sparsity_pattern.h> 22 #include <deal.II/lac/sparsity_tools.h> 28 #ifdef DEAL_II_WITH_MPI 29 # include <deal.II/base/mpi.h> 30 # include <deal.II/base/utilities.h> 32 # include <deal.II/lac/block_sparsity_pattern.h> 33 # include <deal.II/lac/dynamic_sparsity_pattern.h> 36 #ifdef DEAL_II_WITH_METIS 43 #ifdef DEAL_II_TRILINOS_WITH_ZOLTAN 44 # include <zoltan_cpp.h> 50 DEAL_II_NAMESPACE_OPEN
58 const std::vector<unsigned int> &cell_weights,
59 const unsigned int n_partitions,
60 std::vector<unsigned int> & partition_indices)
64 #ifndef DEAL_II_WITH_METIS 65 (void)sparsity_pattern;
68 (void)partition_indices;
77 idx_t n =
static_cast<signed int>(sparsity_pattern.
n_rows()),
80 static_cast<int>(n_partitions),
89 nparts = std::min(n, nparts);
92 idx_t options[METIS_NOPTIONS];
93 METIS_SetDefaultOptions(options);
97 std::vector<idx_t> int_rowstart(1);
98 int_rowstart.reserve(sparsity_pattern.
n_rows() + 1);
99 std::vector<idx_t> int_colnums;
105 col < sparsity_pattern.
end(row);
107 int_colnums.push_back(col->column());
108 int_rowstart.push_back(int_colnums.size());
111 std::vector<idx_t> int_partition_indices(sparsity_pattern.
n_rows());
114 std::vector<idx_t> int_cell_weights;
115 if (cell_weights.size() > 0)
117 Assert(cell_weights.size() == sparsity_pattern.
n_rows(),
119 sparsity_pattern.
n_rows()));
120 int_cell_weights.resize(cell_weights.size());
121 std::copy(cell_weights.begin(),
123 int_cell_weights.begin());
127 idx_t *
const p_int_cell_weights =
128 (cell_weights.size() > 0 ? int_cell_weights.data() :
nullptr);
138 ierr = METIS_PartGraphRecursive(&n,
150 int_partition_indices.data());
154 ierr = METIS_PartGraphKway(&n,
166 int_partition_indices.data());
173 std::copy(int_partition_indices.begin(),
174 int_partition_indices.end(),
175 partition_indices.begin());
181 #ifdef DEAL_II_TRILINOS_WITH_ZOLTAN 184 get_number_of_objects(
void *data,
int *ierr)
195 get_object_list(
void *data,
198 ZOLTAN_ID_PTR globalID,
199 ZOLTAN_ID_PTR localID,
211 auto n_dofs = graph->
n_rows();
213 for (
unsigned int i = 0; i < n_dofs; i++)
222 get_num_edges_list(
void *data,
226 ZOLTAN_ID_PTR globalID,
237 for (
int i = 0; i < num_obj; ++i)
240 numEdges[i] = graph->
row_length(globalID[i]) - 1;
249 get_edge_list(
void *data,
256 ZOLTAN_ID_PTR nborGID,
265 ZOLTAN_ID_PTR nextNborGID = nborGID;
266 int * nextNborProc = nborProc;
270 i < static_cast<SparsityPattern::size_type>(num_obj);
278 if (i != col->column())
283 *nextNborGID++ = col->column();
293 const std::vector<unsigned int> &cell_weights,
294 const unsigned int n_partitions,
295 std::vector<unsigned int> & partition_indices)
299 #ifndef DEAL_II_TRILINOS_WITH_ZOLTAN 300 (void)sparsity_pattern;
303 (void)partition_indices;
308 cell_weights.size() == 0,
310 "The cell weighting functionality for Zoltan has not yet been implemented."));
314 std::unique_ptr<Zoltan> zz =
315 std_cxx14::make_unique<Zoltan>(MPI_COMM_SELF);
319 zz->Set_Param(
"DEBUG_LEVEL",
"0");
323 zz->Set_Param(
"NUM_LOCAL_PARTS",
324 std::to_string(n_partitions));
339 zz->Set_Param(
"PHG_EDGE_SIZE_THRESHOLD",
"0.5");
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);
353 int num_gid_entries = 1;
354 int num_lid_entries = 1;
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;
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;
367 const int rc = zz->LB_Partition(changes,
388 std::fill(partition_indices.begin(), partition_indices.end(), 0);
391 for (
int i = 0; i < num_export; i++)
392 partition_indices[export_local_ids[i]] = export_to_part[i];
400 const unsigned int n_partitions,
401 std::vector<unsigned int> &partition_indices,
404 std::vector<unsigned int> cell_weights;
417 const std::vector<unsigned int> &cell_weights,
418 const unsigned int n_partitions,
419 std::vector<unsigned int> & partition_indices,
428 Assert(partition_indices.size() == sparsity_pattern.
n_rows(),
430 sparsity_pattern.
n_rows()));
433 if (n_partitions == 1 || (sparsity_pattern.
n_rows() == 1))
435 std::fill_n(partition_indices.begin(), partition_indices.size(), 0U);
440 partition_metis(sparsity_pattern,
445 partition_zoltan(sparsity_pattern,
456 std::vector<unsigned int> &color_indices)
460 #ifndef DEAL_II_TRILINOS_WITH_ZOLTAN 461 (void)sparsity_pattern;
467 std::unique_ptr<Zoltan> zz = std_cxx14::make_unique<Zoltan>(MPI_COMM_SELF);
471 zz->Set_Param(
"DEBUG_LEVEL",
"0");
472 zz->Set_Param(
"COLORING_PROBLEM",
"DISTANCE-1");
473 zz->Set_Param(
"NUM_GID_ENTRIES",
"1");
474 zz->Set_Param(
"NUM_LID_ENTRIES",
"1");
475 zz->Set_Param(
"OBJ_WEIGHT_DIM",
"0");
476 zz->Set_Param(
"RECOLORING_NUM_OF_ITERATIONS",
"0");
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);
489 int num_gid_entries = 1;
490 const int num_objects = graph.
n_rows();
493 std::vector<ZOLTAN_ID_TYPE> global_ids(num_objects);
494 std::vector<int> color_exp(num_objects);
497 for (
int i = 0; i < num_objects; i++)
501 int rc = zz->Color(num_gid_entries,
511 color_indices.resize(num_objects);
512 Assert(color_exp.size() == color_indices.size(),
515 std::copy(color_exp.begin(), color_exp.end(), color_indices.begin());
517 unsigned int n_colors =
518 *(std::max_element(color_indices.begin(), color_indices.end()));
532 find_unnumbered_starting_index(
534 const std::vector<DynamicSparsityPattern::size_type> &new_indices)
544 if (sparsity.
row_length(row) < min_coordination)
547 starting_point = row;
572 return starting_point;
581 std::vector<DynamicSparsityPattern::size_type> & new_indices,
582 const std::vector<DynamicSparsityPattern::size_type> &starting_indices)
590 "You can't specify more starting indices than there are rows"));
594 "Only valid for sparsity patterns which store all rows."));
597 ExcMessage(
"Invalid starting index: All starting indices need " 598 "to be between zero and the number of rows in the " 599 "sparsity pattern."));
603 std::vector<DynamicSparsityPattern::size_type> last_round_dofs(
607 std::fill(new_indices.begin(),
613 if (last_round_dofs.empty())
614 last_round_dofs.push_back(
615 internal::find_unnumbered_starting_index(sparsity, new_indices));
623 new_indices[last_round_dofs[i]] = next_free_number++;
629 std::vector<DynamicSparsityPattern::size_type> next_round_dofs;
633 i < last_round_dofs.size();
636 sparsity.
begin(last_round_dofs[i]);
637 j < sparsity.
end(last_round_dofs[i]);
639 next_round_dofs.push_back(j->column());
642 std::sort(next_round_dofs.begin(), next_round_dofs.end());
645 std::vector<DynamicSparsityPattern::size_type>::iterator end_sorted;
647 std::unique(next_round_dofs.begin(), next_round_dofs.end());
648 next_round_dofs.erase(end_sorted, next_round_dofs.end());
651 for (
int s = next_round_dofs.size() - 1; s >= 0; --s)
653 next_round_dofs.erase(next_round_dofs.begin() + s);
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(
680 internal::find_unnumbered_starting_index(sparsity,
688 std::multimap<DynamicSparsityPattern::size_type, int>
689 dofs_by_coordination;
692 for (std::vector<DynamicSparsityPattern::size_type>::iterator s =
693 next_round_dofs.begin();
694 s != next_round_dofs.end();
701 const std::pair<const DynamicSparsityPattern::size_type, int>
702 new_entry(coordination, *s);
703 dofs_by_coordination.insert(new_entry);
707 std::multimap<DynamicSparsityPattern::size_type, int>::iterator i;
708 for (i = dofs_by_coordination.begin(); i != dofs_by_coordination.end();
710 new_indices[i->second] = next_free_number++;
713 last_round_dofs = next_round_dofs;
719 Assert((std::find(new_indices.begin(),
722 (next_free_number == sparsity.
n_rows()),
733 std::vector<DynamicSparsityPattern::size_type> &renumbering)
740 "Only valid for sparsity patterns which store all rows."));
742 std::vector<types::global_dof_index> touched_nodes(
744 std::vector<unsigned int> row_lengths(connectivity.
n_rows());
745 std::set<types::global_dof_index> current_neighbors;
746 std::vector<std::vector<types::global_dof_index>> groups;
756 row_lengths[row] = connectivity.
row_length(row);
759 std::vector<unsigned int> n_remaining_neighbors(row_lengths);
772 std::pair<types::global_dof_index, types::global_dof_index>
777 if (row_lengths[i] < min_neighbors.second)
779 min_neighbors = std::make_pair(i, n_remaining_neighbors[i]);
780 if (n_remaining_neighbors[i] <= 1)
788 current_neighbors.clear();
789 current_neighbors.insert(min_neighbors.first);
790 while (!current_neighbors.empty())
796 for (std::set<types::global_dof_index>::iterator it =
797 current_neighbors.begin();
798 it != current_neighbors.end();
803 if (n_remaining_neighbors[*it] < min_neighbors.second)
805 std::make_pair(*it, n_remaining_neighbors[*it]);
812 min_neighbors.second;
813 for (std::set<types::global_dof_index>::iterator it =
814 current_neighbors.begin();
815 it != current_neighbors.end();
817 if (n_remaining_neighbors[*it] == best_row_length)
818 if (row_lengths[*it] > min_neighbors.second)
819 min_neighbors = std::make_pair(*it, row_lengths[*it]);
823 groups.emplace_back();
824 std::vector<types::global_dof_index> &next_group = groups.back();
826 next_group.push_back(min_neighbors.first);
827 touched_nodes[min_neighbors.first] = groups.size() - 1;
829 connectivity.
begin(min_neighbors.first);
830 it != connectivity.
end(min_neighbors.first);
834 next_group.push_back(it->column());
835 touched_nodes[it->column()] = groups.size() - 1;
843 for (
unsigned int i = 0; i < next_group.size(); ++i)
846 connectivity.
begin(next_group[i]);
847 it != connectivity.
end(next_group[i]);
850 if (touched_nodes[it->column()] ==
852 current_neighbors.insert(it->column());
853 n_remaining_neighbors[it->column()]--;
855 current_neighbors.erase(next_group[i]);
866 if (groups.size() < connectivity.
n_rows())
874 connectivity.
begin(groups[i][col]);
875 it != connectivity.
end(groups[i][col]);
877 connectivity_next.
add(i, touched_nodes[it->column()]);
880 std::vector<types::global_dof_index> renumbering_next(groups.size());
887 col < groups[renumbering_next[i]].size();
889 renumbering[count] = groups[renumbering_next[i]][col];
898 renumbering[count] = groups[i][col];
906 std::vector<DynamicSparsityPattern::size_type> &renumbering)
911 internal::reorder_hierarchical(connectivity, renumbering);
917 #ifdef DEAL_II_WITH_MPI 921 const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu,
922 const MPI_Comm & mpi_comm,
926 std::vector<DynamicSparsityPattern::size_type> start_index(
927 rows_per_cpu.size() + 1);
930 start_index[i + 1] = start_index[i] + rows_per_cpu[i];
933 std::vector<DynamicSparsityPattern::size_type>>;
938 unsigned int dest_cpu = 0;
940 DynamicSparsityPattern::size_type n_local_rel_rows = myrange.
n_elements();
941 for (DynamicSparsityPattern::size_type row_idx = 0;
942 row_idx < n_local_rel_rows;
945 DynamicSparsityPattern::size_type row =
949 while (row >= start_index[dest_cpu + 1])
953 if (dest_cpu == myid)
955 row_idx += rows_per_cpu[myid] - 1;
959 DynamicSparsityPattern::size_type rlen = dsp.
row_length(row);
966 std::vector<DynamicSparsityPattern::size_type> &dst =
971 for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
974 DynamicSparsityPattern::size_type column =
976 dst.push_back(column);
981 unsigned int num_receive = 0;
983 std::vector<unsigned int> send_to;
984 send_to.reserve(send_data.size());
985 for (map_vec_t::iterator it = send_data.begin(); it != send_data.end();
987 send_to.push_back(it->first);
995 std::vector<MPI_Request> requests(send_data.size());
1000 unsigned int idx = 0;
1001 for (map_vec_t::iterator it = send_data.begin(); it != send_data.end();
1004 const int ierr = MPI_Isend(&(it->second[0]),
1006 DEAL_II_DOF_INDEX_MPI_TYPE,
1017 std::vector<DynamicSparsityPattern::size_type> recv_buf;
1018 for (
unsigned int index = 0; index < num_receive; ++index)
1022 int ierr = MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, mpi_comm, &status);
1026 ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
1028 recv_buf.resize(len);
1029 ierr = MPI_Recv(recv_buf.data(),
1031 DEAL_II_DOF_INDEX_MPI_TYPE,
1038 std::vector<DynamicSparsityPattern::size_type>::const_iterator ptr =
1040 std::vector<DynamicSparsityPattern::size_type>::const_iterator end =
1044 DynamicSparsityPattern::size_type num = *(ptr++);
1046 DynamicSparsityPattern::size_type row = *(ptr++);
1047 for (
unsigned int c = 0; c < num; ++c)
1059 if (requests.size())
1062 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1069 const std::vector<IndexSet> &owned_set_per_cpu,
1070 const MPI_Comm & mpi_comm,
1077 std::vector<BlockDynamicSparsityPattern::size_type>>;
1078 map_vec_t send_data;
1081 unsigned int dest_cpu = 0;
1083 BlockDynamicSparsityPattern::size_type n_local_rel_rows =
1085 for (BlockDynamicSparsityPattern::size_type row_idx = 0;
1086 row_idx < n_local_rel_rows;
1089 BlockDynamicSparsityPattern::size_type row =
1095 while (!owned_set_per_cpu[dest_cpu].is_element(row))
1098 if (dest_cpu == owned_set_per_cpu.size())
1103 if (dest_cpu == myid)
1106 BlockDynamicSparsityPattern::size_type rlen = dsp.
row_length(row);
1113 std::vector<BlockDynamicSparsityPattern::size_type> &dst =
1114 send_data[dest_cpu];
1116 dst.push_back(rlen);
1118 for (BlockDynamicSparsityPattern::size_type c = 0; c < rlen; ++c)
1121 BlockDynamicSparsityPattern::size_type column =
1123 dst.push_back(column);
1128 unsigned int num_receive = 0;
1130 std::vector<unsigned int> send_to;
1131 send_to.reserve(send_data.size());
1132 for (map_vec_t::iterator it = send_data.begin(); it != send_data.end();
1134 send_to.push_back(it->first);
1142 std::vector<MPI_Request> requests(send_data.size());
1147 unsigned int idx = 0;
1148 for (map_vec_t::iterator it = send_data.begin(); it != send_data.end();
1151 const int ierr = MPI_Isend(&(it->second[0]),
1153 DEAL_II_DOF_INDEX_MPI_TYPE,
1164 std::vector<BlockDynamicSparsityPattern::size_type> recv_buf;
1165 for (
unsigned int index = 0; index < num_receive; ++index)
1169 int ierr = MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, mpi_comm, &status);
1173 ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
1175 recv_buf.resize(len);
1176 ierr = MPI_Recv(recv_buf.data(),
1178 DEAL_II_DOF_INDEX_MPI_TYPE,
1185 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1186 ptr = recv_buf.begin();
1187 std::vector<BlockDynamicSparsityPattern::size_type>::const_iterator
1188 end = recv_buf.end();
1191 BlockDynamicSparsityPattern::size_type num = *(ptr++);
1193 BlockDynamicSparsityPattern::size_type row = *(ptr++);
1194 for (
unsigned int c = 0; c < num; ++c)
1206 if (requests.size())
1209 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1216 DEAL_II_NAMESPACE_CLOSE
size_type row_length(const size_type row) const
const types::global_dof_index invalid_size_type
#define AssertDimension(dim1, dim2)
void add(const size_type i, const size_type j)
static::ExceptionBase & ExcMETISNotInstalled()
size_type n_elements() const
#define AssertThrow(cond, exc)
unsigned long long int global_dof_index
static::ExceptionBase & ExcNotCompressed()
static::ExceptionBase & ExcMessage(std::string arg1)
unsigned int row_length(const size_type row) const
static::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
void add(const size_type i, const size_type j)
const IndexSet & row_index_set() const
#define Assert(cond, exc)
static::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int row_length(const size_type row) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
std::vector< unsigned int > invert_permutation(const std::vector< unsigned int > &permutation)
types::global_dof_index size_type
size_type column_number(const size_type row, const size_type index) const
static::ExceptionBase & ExcNotQuadratic()
types::global_dof_index size_type
static::ExceptionBase & ExcZOLTANNotInstalled()
#define AssertThrowMPI(error_code)
size_type column_number(const size_type row, const unsigned int index) const
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
bool exists(const size_type i, const size_type j) const
const types::global_dof_index invalid_dof_index
size_type nth_index_in_set(const unsigned int local_index) const
std::vector< unsigned int > compute_point_to_point_communication_pattern(const MPI_Comm &mpi_comm, const std::vector< unsigned int > &destinations)
bool is_compressed() const
types::global_dof_index size_type
std::size_t n_nonzero_elements() const
static::ExceptionBase & ExcMETISError(int arg1)
static::ExceptionBase & ExcInternalError()