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 #ifndef dealii_sparsity_tools_h 17 #define dealii_sparsity_tools_h 18 19 20 #include <deal.II/base/config.h> 21 22 #include <deal.II/base/exceptions.h> 23 24 #include <deal.II/lac/block_sparsity_pattern.h> 25 #include <deal.II/lac/dynamic_sparsity_pattern.h> 26 #include <deal.II/lac/sparsity_pattern.h> 27 28 #include <memory> 29 #include <vector> 30 31 #ifdef DEAL_II_WITH_MPI 32 # include <deal.II/base/index_set.h> 33 34 # include <mpi.h> 35 #endif 36 37 DEAL_II_NAMESPACE_OPEN 38 39 40 /*! @addtogroup Sparsity 41 *@{ 42 */ 43 44 /** 45 * A namespace for functions that deal with things that one can do on sparsity 46 * patterns, such as renumbering rows and columns (or degrees of freedom if 47 * you want) according to the connectivity, or partitioning degrees of 48 * freedom. 49 */ 50 namespace SparsityTools 51 { 52 /** 53 * Enumerator with options for partitioner 54 */ 55 enum class Partitioner 56 { 57 /** 58 * Use METIS partitioner. 59 */ 60 metis = 0, 61 /** 62 * Use ZOLTAN partitioner. 63 */ 64 zoltan 65 }; 66 67 68 /** 69 * Use a partitioning algorithm to generate a partitioning of the degrees of 70 * freedom represented by this sparsity pattern. In effect, we view this 71 * sparsity pattern as a graph of connections between various degrees of 72 * freedom, where each nonzero entry in the sparsity pattern corresponds to 73 * an edge between two nodes in the connection graph. The goal is then to 74 * decompose this graph into groups of nodes so that a minimal number of 75 * edges are cut by the boundaries between node groups. This partitioning is 76 * done by METIS or ZOLTAN, depending upon which partitioner is chosen 77 * in the fourth argument. The default is METIS. Note that METIS and 78 * ZOLTAN can only partition symmetric sparsity patterns, and that of 79 * course the sparsity pattern has to be square. We do not check for 80 * symmetry of the sparsity pattern, since this is an expensive operation, 81 * but rather leave this as the responsibility of caller of this function. 82 * 83 * After calling this function, the output array will have values between 84 * zero and @p n_partitions-1 for each node (i.e. row or column of the 85 * matrix). 86 * 87 * If deal.II was not installed with packages ZOLTAN or METIS, this 88 * function will generate an error when corresponding partition method 89 * is chosen, unless @p n_partitions is one. I.e., you can write a program 90 * so that it runs in the single-processor single-partition case without 91 * the packages installed, and only requires them installed when 92 * multiple partitions are required. 93 * 94 * Note that the sparsity pattern itself is not changed by calling this 95 * function. However, you will likely use the information generated by 96 * calling this function to renumber degrees of freedom, after which you 97 * will of course have to regenerate the sparsity pattern. 98 * 99 * This function will rarely be called separately, since in finite element 100 * methods you will want to partition the mesh, not the matrix. This can be 101 * done by calling @p GridTools::partition_triangulation. 102 */ 103 void 104 partition(const SparsityPattern & sparsity_pattern, 105 const unsigned int n_partitions, 106 std::vector<unsigned int> &partition_indices, 107 const Partitioner partitioner = Partitioner::metis); 108 109 110 /** 111 * This function performs the same operation as the one above, except that 112 * it takes into consideration a set of @p cell_weights, which allow the 113 * partitioner to balance the graph while taking into consideration the 114 * computational effort expended on each cell. 115 * 116 * @note If the @p cell_weights vector is empty, then no weighting is taken 117 * into consideration. If not then the size of this vector must equal to the 118 * number of active cells in the triangulation. 119 */ 120 void 121 partition(const SparsityPattern & sparsity_pattern, 122 const std::vector<unsigned int> &cell_weights, 123 const unsigned int n_partitions, 124 std::vector<unsigned int> & partition_indices, 125 const Partitioner partitioner = Partitioner::metis); 126 127 /** 128 * Using a coloring algorithm provided by ZOLTAN to color nodes whose 129 * connections are represented using a SparsityPattern object. In effect, 130 * we view this sparsity pattern as a graph of connections between nodes, 131 * where each nonzero entry in the @p sparsity_pattern corresponds to 132 * an edge between two nodes. The goal is to assign each node a color index 133 * such that no two directly connected nodes have the same color. 134 * The assigned colors are listed in @p color_indices indexed from one and 135 * the function returns total number of colors used. ZOLTAN coloring 136 * algorithm is run in serial by each processor. Hence all processors have 137 * coloring information of all the nodes. A wrapper function to this 138 * function is available in GraphColoring namespace as well. 139 * 140 * Note that this function requires that SparsityPattern be symmetric 141 * (and hence square) i.e an undirected graph representation. We do not 142 * check for symmetry of the sparsity pattern, since this is an expensive 143 * operation, but rather leave this as the responsibility of caller of 144 * this function. 145 * 146 * @image html color_undirected_graph.png 147 * The usage of the function is illustrated by the image above, showing five 148 * nodes numbered from 0 to 4. The connections shown are bidirectional. 149 * After coloring, it is clear that no two directly connected nodes are 150 * assigned the same color. 151 * 152 * If deal.II was not installed with package ZOLTAN, this function will 153 * generate an error. 154 * 155 * @note The current function is an alternative to 156 * GraphColoring::make_graph_coloring() which is tailored to graph 157 * coloring arising in shared-memory parallel assembly of matrices. 158 */ 159 unsigned int 160 color_sparsity_pattern(const SparsityPattern & sparsity_pattern, 161 std::vector<unsigned int> &color_indices); 162 163 /** 164 * For a given sparsity pattern, compute a re-enumeration of row/column 165 * indices based on the algorithm by Cuthill-McKee. 166 * 167 * This algorithm is a graph renumbering algorithm in which we attempt to 168 * find a new numbering of all nodes of a graph based on their connectivity 169 * to other nodes (i.e. the edges that connect nodes). This connectivity is 170 * here represented by the sparsity pattern. In many cases within the 171 * library, the nodes represent degrees of freedom and edges are nonzero 172 * entries in a matrix, i.e. pairs of degrees of freedom that couple through 173 * the action of a bilinear form. 174 * 175 * The algorithms starts at a node, searches the other nodes for those which 176 * are coupled with the one we started with and numbers these in a certain 177 * way. It then finds the second level of nodes, namely those that couple 178 * with those of the previous level (which were those that coupled with the 179 * initial node) and numbers these. And so on. For the details of the 180 * algorithm, especially the numbering within each level, we refer the 181 * reader to the book of Schwarz (H. R. Schwarz: Methode der finiten 182 * Elemente). 183 * 184 * These algorithms have one major drawback: they require a good starting 185 * node, i.e. node that will have number zero in the output array. A 186 * starting node forming the initial level of nodes can thus be given by the 187 * user, e.g. by exploiting knowledge of the actual topology of the domain. 188 * It is also possible to give several starting indices, which may be used 189 * to simulate a simple upstream numbering (by giving the inflow nodes as 190 * starting values) or to make preconditioning faster (by letting the 191 * Dirichlet boundary indices be starting points). 192 * 193 * If no starting index is given, one is chosen automatically, namely one 194 * with the smallest coordination number (the coordination number is the 195 * number of other nodes this node couples with). This node is usually 196 * located on the boundary of the domain. There is, however, large ambiguity 197 * in this when using the hierarchical meshes used in this library, since in 198 * most cases the computational domain is not approximated by tilting and 199 * deforming elements and by plugging together variable numbers of elements 200 * at vertices, but rather by hierarchical refinement. There is therefore a 201 * large number of nodes with equal coordination numbers. The renumbering 202 * algorithms will therefore not give optimal results. 203 * 204 * If the graph has two or more unconnected components and if no starting 205 * indices are given, the algorithm will number each component 206 * consecutively. However, this requires the determination of a starting 207 * index for each component; as a consequence, the algorithm will produce an 208 * exception if starting indices are given, taking the latter as an 209 * indication that the caller of the function would like to override the 210 * part of the algorithm that chooses starting indices. 211 */ 212 void 213 reorder_Cuthill_McKee( 214 const DynamicSparsityPattern & sparsity, 215 std::vector<DynamicSparsityPattern::size_type> & new_indices, 216 const std::vector<DynamicSparsityPattern::size_type> &starting_indices = 217 std::vector<DynamicSparsityPattern::size_type>()); 218 219 /** 220 * For a given sparsity pattern, compute a re-enumeration of row/column 221 * indices in a hierarchical way, similar to what 222 * DoFRenumbering::hierarchical does for degrees of freedom on 223 * hierarchically refined meshes. 224 * 225 * This algorithm first selects a node with the minimum number of neighbors 226 * and puts that node and its direct neighbors into one chunk. Next, it 227 * selects one of the neighbors of the already selected nodes, adds the node 228 * and its direct neighbors that are not part of one of the previous chunks, 229 * into the next. After this sweep, neighboring nodes are grouped together. 230 * To ensure a similar grouping on a more global level, this grouping is 231 * called recursively on the groups so formed. The recursion stops when no 232 * further grouping is possible. Eventually, the ordering obtained by this 233 * method passes through the indices represented in the sparsity pattern in 234 * a z-like way. 235 * 236 * If the graph has two or more unconnected components, the algorithm will 237 * number each component consecutively, starting with the components with 238 * the lowest number of nodes. 239 */ 240 void 241 reorder_hierarchical( 242 const DynamicSparsityPattern & sparsity, 243 std::vector<DynamicSparsityPattern::size_type> &new_indices); 244 245 #ifdef DEAL_II_WITH_MPI 246 /** 247 * Communicate rows in a dynamic sparsity pattern over MPI. 248 * 249 * @param[in,out] dsp A dynamic sparsity pattern that has been built locally 250 * and for which we need to exchange entries with other processors to make 251 * sure that each processor knows all the elements of the rows of a matrix 252 * it stores and that may eventually be written to. This sparsity pattern 253 * will be changed as a result of this function: All entries in rows that 254 * belong to a different processor are sent to them and added there. 255 * 256 * @param locally_owned_rows An IndexSet describing the rows owned by the 257 * calling MPI process. The index set shall be one-to-one among the 258 * processors in the communicator. 259 * 260 * @param mpi_comm The MPI communicator shared between the processors that 261 * participate in this operation. 262 * 263 * @param locally_relevant_rows The range of elements stored on the local 264 * MPI process. This should be the one used in the constructor of the 265 * DynamicSparsityPattern, and should also be the locally relevant set. Only 266 * rows contained in this set are checked in dsp for transfer. This function 267 * needs to be used with PETScWrappers::MPI::SparseMatrix for it to work 268 * correctly in a parallel computation. 269 */ 270 void 271 distribute_sparsity_pattern(DynamicSparsityPattern &dsp, 272 const IndexSet & locally_owned_rows, 273 const MPI_Comm & mpi_comm, 274 const IndexSet & locally_relevant_rows); 275 276 /** 277 * Communicate rows in a dynamic sparsity pattern over MPI, similar to the 278 * one above but using a vector `rows_per_cpu` containing the number of 279 * rows per CPU for determining ownership. This is typically the value 280 * returned by DoFHandler::n_locally_owned_dofs_per_processor -- given that 281 * the construction of the input to this function involves all-to-all 282 * communication, it is typically slower than the function above for more 283 * than a thousand of processes (and quick enough also for small sizes). 284 */ 285 void 286 distribute_sparsity_pattern( 287 DynamicSparsityPattern & dsp, 288 const std::vector<DynamicSparsityPattern::size_type> &rows_per_cpu, 289 const MPI_Comm & mpi_comm, 290 const IndexSet & myrange); 291 292 /** 293 * Similar to the function above, but for BlockDynamicSparsityPattern 294 * instead. 295 * 296 * @param[in,out] dsp The locally built sparsity pattern to be modified. 297 * 298 * @param locally_owned_rows An IndexSet describing the rows owned by the 299 * calling MPI process. The index set shall be one-to-one among the 300 * processors in the communicator. 301 * 302 * @param mpi_comm The MPI communicator to use. 303 * 304 * @param locally_relevant_rows Typically the locally relevant DoFs. 305 */ 306 void 307 distribute_sparsity_pattern(BlockDynamicSparsityPattern &dsp, 308 const IndexSet & locally_owned_rows, 309 const MPI_Comm & mpi_comm, 310 const IndexSet &locally_relevant_rows); 311 312 /** 313 * @deprecated Use the distribute_sparsity_pattern() with a single index set 314 * for the present MPI process only. 315 */ 316 void 317 distribute_sparsity_pattern(BlockDynamicSparsityPattern &dsp, 318 const std::vector<IndexSet> &owned_set_per_cpu, 319 const MPI_Comm & mpi_comm, 320 const IndexSet & myrange); 321 322 /** 323 * Gather rows in a dynamic sparsity pattern over MPI. 324 * The function is similar to SparsityTools::distribute(), 325 * however instead of distributing sparsity stored in non-owned 326 * rows on this MPI process, this function will gather sparsity 327 * from other MPI processes and will add this to the local 328 * DynamicSparsityPattern. 329 * 330 * @param dsp A dynamic sparsity pattern that has been built locally and which 331 * we need to extend according to the sparsity of rows stored on other MPI 332 * processes. 333 * 334 * @param locally_owned_rows An IndexSet describing the rows owned by the 335 * calling MPI process. The index set shall be one-to-one among the 336 * processors in the communicator. 337 * 338 * @param mpi_comm The MPI communicator shared between the processors that 339 * participate in this operation. 340 * 341 * @param locally_relevant_rows The range of rows this MPI process needs to 342 * gather. Only the part which is not included in the locally owned rows will 343 * be used. 344 */ 345 void 346 gather_sparsity_pattern(DynamicSparsityPattern &dsp, 347 const IndexSet & locally_owned_rows, 348 const MPI_Comm & mpi_comm, 349 const IndexSet & locally_relevant_rows); 350 351 /** 352 * @deprecated Use the gather_sparsity_pattern() method with the index set 353 * for the present processor only. 354 */ 355 DEAL_II_DEPRECATED void 356 gather_sparsity_pattern(DynamicSparsityPattern & dsp, 357 const std::vector<IndexSet> &owned_rows_per_processor, 358 const MPI_Comm & mpi_comm, 359 const IndexSet & ghost_range); 360 361 #endif 362 363 364 /** 365 * Exception 366 */ 367 DeclExceptionMsg(ExcMETISNotInstalled, 368 "The function you called requires METIS, but you did not " 369 "configure deal.II with METIS."); 370 371 /** 372 * Exception 373 */ 374 DeclException1(ExcInvalidNumberOfPartitions, 375 int, 376 << "The number of partitions you gave is " << arg1 377 << ", but must be greater than zero."); 378 379 /** 380 * Exception 381 */ 382 DeclException1(ExcMETISError, 383 int, 384 << " An error with error number " << arg1 385 << " occurred while calling a METIS function"); 386 387 /** 388 * Exception 389 */ 390 DeclException2(ExcInvalidArraySize, 391 int, 392 int, 393 << "The array has size " << arg1 << " but should have size " 394 << arg2); 395 /** 396 * Exception 397 */ 398 DeclExceptionMsg( 399 ExcZOLTANNotInstalled, 400 "The function you called requires ZOLTAN, but you did not " 401 "configure deal.II with ZOLTAN or zoltan_cpp.h is not available."); 402 } // namespace SparsityTools 403 404 /** 405 * @} 406 */ 407 408 DEAL_II_NAMESPACE_CLOSE 409 410 #endif 411