DOLFINx
DOLFINx C++ interface
IndexMap.h
1// Copyright (C) 2015-2019 Chris Richardson, Garth N. Wells and Igor Baratta
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include <array>
10#include <cstdint>
11#include <dolfinx/common/MPI.h>
12#include <dolfinx/graph/AdjacencyList.h>
13#include <map>
14#include <memory>
15#include <tuple>
16#include <utility>
17#include <vector>
18#include <xtl/xspan.hpp>
19
20namespace dolfinx::common
21{
22// Forward declaration
23class IndexMap;
24
32std::vector<int32_t>
33compute_owned_indices(const xtl::span<const std::int32_t>& indices,
34 const IndexMap& map);
35
45std::tuple<std::int64_t, std::vector<std::int32_t>,
46 std::vector<std::vector<std::int64_t>>,
47 std::vector<std::vector<int>>>
49 const std::vector<
50 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
51
58
60{
61public:
63 enum class Mode
64 {
65 insert,
66 add
67 };
68
70 enum class Direction
71 {
72 reverse, // Ghost to owner
73 forward, // Owner to ghost
74 };
75
83 IndexMap(MPI_Comm comm, std::int32_t local_size);
84
97 IndexMap(MPI_Comm comm, std::int32_t local_size,
98 const xtl::span<const int>& dest_ranks,
99 const xtl::span<const std::int64_t>& ghosts,
100 const xtl::span<const int>& src_ranks);
101
102private:
103 template <typename U, typename V, typename W, typename X>
104 IndexMap(std::array<std::int64_t, 2> local_range, std::size_t size_global,
105 MPI_Comm comm, U&& comm_owner_to_ghost, U&& comm_ghost_to_owner,
106 V&& displs_recv_fwd, V&& ghost_pos_recv_fwd, W&& ghosts,
107 X&& shared_indices)
108 : _local_range(local_range), _size_global(size_global), _comm(comm),
109 _comm_owner_to_ghost(std::forward<U>(comm_owner_to_ghost)),
110 _comm_ghost_to_owner(std::forward<U>(comm_ghost_to_owner)),
111 _displs_recv_fwd(std::forward<V>(displs_recv_fwd)),
112 _ghost_pos_recv_fwd(std::forward<V>(ghost_pos_recv_fwd)),
113 _ghosts(std::forward<W>(ghosts)),
114 _shared_indices(std::forward<X>(shared_indices))
115 {
116 _sizes_recv_fwd.resize(_displs_recv_fwd.size() - 1, 0);
117 std::adjacent_difference(_displs_recv_fwd.cbegin() + 1,
118 _displs_recv_fwd.cend(), _sizes_recv_fwd.begin());
119
120 const std::vector<int32_t>& displs_send = _shared_indices->offsets();
121 _sizes_send_fwd.resize(_shared_indices->num_nodes(), 0);
122 std::adjacent_difference(displs_send.cbegin() + 1, displs_send.cend(),
123 _sizes_send_fwd.begin());
124 }
125
126public:
127 // Copy constructor
128 IndexMap(const IndexMap& map) = delete;
129
131 IndexMap(IndexMap&& map) = default;
132
134 ~IndexMap() = default;
135
137 IndexMap& operator=(IndexMap&& map) = default;
138
139 // Copy assignment
140 IndexMap& operator=(const IndexMap& map) = delete;
141
143 std::array<std::int64_t, 2> local_range() const noexcept;
144
146 std::int32_t num_ghosts() const noexcept;
147
149 std::int32_t size_local() const noexcept;
150
152 std::int64_t size_global() const noexcept;
153
156 const std::vector<std::int64_t>& ghosts() const noexcept;
157
160 MPI_Comm comm() const;
161
166 MPI_Comm comm(Direction dir) const;
167
171 void local_to_global(const xtl::span<const std::int32_t>& local,
172 const xtl::span<std::int64_t>& global) const;
173
178 void global_to_local(const xtl::span<const std::int64_t>& global,
179 const xtl::span<std::int32_t>& local) const;
180
184 std::vector<std::int64_t> global_indices() const;
185
199 const graph::AdjacencyList<std::int32_t>&
200 scatter_fwd_indices() const noexcept;
201
207 const std::vector<std::int32_t>& scatter_fwd_ghost_positions() const noexcept;
208
222 std::vector<int> ghost_owners() const;
223
230 std::map<std::int32_t, std::set<int>> compute_shared_indices() const;
231
242 std::pair<IndexMap, std::vector<std::int32_t>>
243 create_submap(const xtl::span<const std::int32_t>& indices) const;
244
264 template <typename T>
265 void scatter_fwd_begin(const xtl::span<const T>& send_buffer,
266 MPI_Datatype& data_type, MPI_Request& request,
267 const xtl::span<T>& recv_buffer) const
268 {
269 // Send displacement
270 const std::vector<int32_t>& displs_send_fwd = _shared_indices->offsets();
271
272 // Return early if there are no incoming or outgoing edges
273 if (_displs_recv_fwd.size() == 1 and displs_send_fwd.size() == 1)
274 return;
275
276 // Get block size
277 int n;
278 MPI_Type_size(data_type, &n);
279 n /= sizeof(T);
280 if (static_cast<int>(send_buffer.size()) != n * displs_send_fwd.back())
281 throw std::runtime_error("Incompatible send buffer size.");
282 if (static_cast<int>(recv_buffer.size()) != n * _displs_recv_fwd.back())
283 throw std::runtime_error("Incompatible receive buffer size..");
284
285 // Start send/receive
286 MPI_Ineighbor_alltoallv(send_buffer.data(), _sizes_send_fwd.data(),
287 displs_send_fwd.data(), data_type,
288 recv_buffer.data(), _sizes_recv_fwd.data(),
289 _displs_recv_fwd.data(), data_type,
290 _comm_owner_to_ghost.comm(), &request);
291 }
292
299 void scatter_fwd_end(MPI_Request& request) const
300 {
301 // Return early if there are no incoming or outgoing edges
302 const std::vector<int32_t>& displs_send_fwd = _shared_indices->offsets();
303 if (_displs_recv_fwd.size() == 1 and displs_send_fwd.size() == 1)
304 return;
305
306 // Wait for communication to complete
307 MPI_Wait(&request, MPI_STATUS_IGNORE);
308 }
309
320 template <typename T>
321 void scatter_fwd(const xtl::span<const T>& local_data,
322 xtl::span<T> remote_data, int n) const
323 {
324 MPI_Datatype data_type;
325 if (n == 1)
326 data_type = dolfinx::MPI::mpi_type<T>();
327 else
328 {
329 MPI_Type_contiguous(n, dolfinx::MPI::mpi_type<T>(), &data_type);
330 MPI_Type_commit(&data_type);
331 }
332
333 const std::vector<std::int32_t>& indices = _shared_indices->array();
334 std::vector<T> send_buffer(n * indices.size());
335 for (std::size_t i = 0; i < indices.size(); ++i)
336 {
337 std::copy_n(std::next(local_data.cbegin(), n * indices[i]), n,
338 std::next(send_buffer.begin(), n * i));
339 }
340
341 MPI_Request request;
342 std::vector<T> buffer_recv(n * _displs_recv_fwd.back());
343 scatter_fwd_begin(xtl::span<const T>(send_buffer), data_type, request,
344 xtl::span<T>(buffer_recv));
345 scatter_fwd_end(request);
346
347 // Copy into ghost area ("remote_data")
348 assert(remote_data.size() == n * _ghost_pos_recv_fwd.size());
349 for (std::size_t i = 0; i < _ghost_pos_recv_fwd.size(); ++i)
350 {
351 std::copy_n(std::next(buffer_recv.cbegin(), n * _ghost_pos_recv_fwd[i]),
352 n, std::next(remote_data.begin(), n * i));
353 }
354
355 if (n != 1)
356 MPI_Type_free(&data_type);
357 }
358
379 template <typename T>
380 void scatter_rev_begin(const xtl::span<const T>& send_buffer,
381 MPI_Datatype& data_type, MPI_Request& request,
382 const xtl::span<T>& recv_buffer) const
383 {
384 // Get displacement vector
385 const std::vector<int32_t>& displs_send_fwd = _shared_indices->offsets();
386
387 // Return early if there are no incoming or outgoing edges
388 if (_displs_recv_fwd.size() == 1 and displs_send_fwd.size() == 1)
389 return;
390
391 // Get block size
392 int n;
393 MPI_Type_size(data_type, &n);
394 n /= sizeof(T);
395 if (static_cast<int>(send_buffer.size()) != n * _ghosts.size())
396 throw std::runtime_error("Inconsistent send buffer size.");
397 if (static_cast<int>(recv_buffer.size()) != n * displs_send_fwd.back())
398 throw std::runtime_error("Inconsistent receive buffer size.");
399
400 // Send and receive data
401 MPI_Ineighbor_alltoallv(send_buffer.data(), _sizes_recv_fwd.data(),
402 _displs_recv_fwd.data(), data_type,
403 recv_buffer.data(), _sizes_send_fwd.data(),
404 displs_send_fwd.data(), data_type,
405 _comm_ghost_to_owner.comm(), &request);
406 }
407
414 void scatter_rev_end(MPI_Request& request) const
415 {
416 // Return early if there are no incoming or outgoing edges
417 const std::vector<int32_t>& displs_send_fwd = _shared_indices->offsets();
418 if (_displs_recv_fwd.size() == 1 and displs_send_fwd.size() == 1)
419 return;
420
421 // Wait for communication to complete
422 MPI_Wait(&request, MPI_STATUS_IGNORE);
423 }
424
434 template <typename T>
435 void scatter_rev(xtl::span<T> local_data,
436 const xtl::span<const T>& remote_data, int n,
437 IndexMap::Mode op) const
438 {
439 MPI_Datatype data_type;
440 if (n == 1)
441 data_type = dolfinx::MPI::mpi_type<T>();
442 else
443 {
444 MPI_Type_contiguous(n, dolfinx::MPI::mpi_type<T>(), &data_type);
445 MPI_Type_commit(&data_type);
446 }
447
448 // Pack send buffer
449 std::vector<T> buffer_send;
450 buffer_send.resize(n * _displs_recv_fwd.back());
451 for (std::size_t i = 0; i < _ghost_pos_recv_fwd.size(); ++i)
452 {
453 std::copy_n(std::next(remote_data.cbegin(), n * i), n,
454 std::next(buffer_send.begin(), n * _ghost_pos_recv_fwd[i]));
455 }
456
457 // Exchange data
458 MPI_Request request;
459 std::vector<T> buffer_recv(n * _shared_indices->array().size());
460 scatter_rev_begin(xtl::span<const T>(buffer_send), data_type, request,
461 xtl::span<T>(buffer_recv));
462 scatter_rev_end(request);
463
464 // Copy or accumulate into "local_data"
465 assert(local_data.size() == n * this->size_local());
466 const std::vector<std::int32_t>& shared_indices = _shared_indices->array();
467 switch (op)
468 {
469 case Mode::insert:
470 for (std::size_t i = 0; i < shared_indices.size(); ++i)
471 {
472 std::copy_n(std::next(buffer_recv.cbegin(), n * i), n,
473 std::next(local_data.begin(), n * shared_indices[i]));
474 }
475 break;
476 case Mode::add:
477 for (std::size_t i = 0; i < shared_indices.size(); ++i)
478 {
479 for (int j = 0; j < n; ++j)
480 local_data[shared_indices[i] * n + j] += buffer_recv[i * n + j];
481 }
482 break;
483 }
484
485 if (n != 1)
486 MPI_Type_free(&data_type);
487 }
488
489private:
490 // Range of indices (global) owned by this process
491 std::array<std::int64_t, 2> _local_range;
492
493 // Number indices across communicator
494 std::int64_t _size_global;
495
496 // MPI communicator (duplicated of 'input' communicator)
497 dolfinx::MPI::Comm _comm;
498
499 // Communicator where the source ranks own the indices in the callers
500 // halo, and the destination ranks 'ghost' indices owned by the
501 // caller. I.e.,
502 // - in-edges (src) are from ranks that own my ghosts
503 // - out-edges (dest) go to ranks that 'ghost' my owned indices
504 dolfinx::MPI::Comm _comm_owner_to_ghost;
505
506 // Communicator where the source ranks have ghost indices that are
507 // owned by the caller, and the destination ranks are the owners of
508 // indices in the callers halo region. I.e.,
509 // - in-edges (src) are from ranks that 'ghost' my owned indicies
510 // - out-edges (dest) are to the owning ranks of my ghost indices
511 dolfinx::MPI::Comm _comm_ghost_to_owner;
512
513 // MPI sizes and displacements for forward (owner -> ghost) scatter
514 // Note: '_displs_send_fwd' can be got from _shared_indices->offsets()
515 std::vector<std::int32_t> _sizes_send_fwd, _sizes_recv_fwd, _displs_recv_fwd;
516
517 // Position in the recv buffer for a forward scatter for the ith ghost
518 // index (_ghost[i]) entry
519 std::vector<std::int32_t> _ghost_pos_recv_fwd;
520
521 // Local-to-global map for ghost indices
522 std::vector<std::int64_t> _ghosts;
523
524 // List of owned local indices that are in the ghost (halo) region on
525 // other ranks, grouped by rank in the neighbor communicator
526 // (destination ranks in forward communicator and source ranks in the
527 // reverse communicator), i.e. `_shared_indices.num_nodes() ==
528 // size(_comm_owner_to_ghost)`. The array _shared_indices.offsets() is
529 // equivalent to 'displs_send_fwd'.
530 std::unique_ptr<graph::AdjacencyList<std::int32_t>> _shared_indices;
531};
532
533} // namespace dolfinx::common
A duplicate MPI communicator and manage lifetime of the communicator.
Definition: MPI.h:40
This class represents the distribution index arrays across processes. An index array is a contiguous ...
Definition: IndexMap.h:60
~IndexMap()=default
Destructor.
Direction
Edge directions of neighborhood communicator.
Definition: IndexMap.h:71
std::vector< int > ghost_owners() const
Compute the owner on the neighborhood communicator of each ghost index.
Definition: IndexMap.cpp:627
const std::vector< std::int32_t > & scatter_fwd_ghost_positions() const noexcept
Position of ghost entries in the receive buffer after a forward scatter, e.g. for a receive buffer b ...
Definition: IndexMap.cpp:622
std::int32_t num_ghosts() const noexcept
Number of ghost indices on this process.
Definition: IndexMap.cpp:542
Mode
Mode for reverse scatter operation.
Definition: IndexMap.h:64
void scatter_rev_end(MPI_Request &request) const
Complete a non-blocking send of ghost values to the owning rank. This function complete the communica...
Definition: IndexMap.h:414
std::array< std::int64_t, 2 > local_range() const noexcept
Range of indices (global) owned by this process.
Definition: IndexMap.cpp:537
std::int32_t size_local() const noexcept
Number of indices owned by on this process.
Definition: IndexMap.cpp:544
std::int64_t size_global() const noexcept
Number indices across communicator.
Definition: IndexMap.cpp:549
IndexMap(IndexMap &&map)=default
Move constructor.
const graph::AdjacencyList< std::int32_t > & scatter_fwd_indices() const noexcept
Local (owned) indices shared with neighbor processes, i.e. are ghosts on other processes,...
Definition: IndexMap.cpp:615
void scatter_fwd(const xtl::span< const T > &local_data, xtl::span< T > remote_data, int n) const
Send n values for each index that is owned to processes that have the index as a ghost....
Definition: IndexMap.h:321
IndexMap & operator=(IndexMap &&map)=default
Move assignment.
const std::vector< std::int64_t > & ghosts() const noexcept
Local-to-global map for ghosts (local indexing beyond end of local range)
Definition: IndexMap.cpp:551
void global_to_local(const xtl::span< const std::int64_t > &global, const xtl::span< std::int32_t > &local) const
Compute local indices for array of global indices.
Definition: IndexMap.cpp:575
std::pair< IndexMap, std::vector< std::int32_t > > create_submap(const xtl::span< const std::int32_t > &indices) const
Create new index map from a subset of indices in this index map. The order of the indices is preserve...
Definition: IndexMap.cpp:785
void scatter_fwd_end(MPI_Request &request) const
Complete a non-blocking send from the local owner of to process ranks that have the index as a ghost....
Definition: IndexMap.h:299
std::map< std::int32_t, std::set< int > > compute_shared_indices() const
Definition: IndexMap.cpp:656
void local_to_global(const xtl::span< const std::int32_t > &local, const xtl::span< std::int64_t > &global) const
Compute global indices for array of local indices.
Definition: IndexMap.cpp:556
IndexMap(MPI_Comm comm, std::int32_t local_size)
Create an non-overlapping index map with local_size owned on this process.
Definition: IndexMap.cpp:371
void scatter_fwd_begin(const xtl::span< const T > &send_buffer, MPI_Datatype &data_type, MPI_Request &request, const xtl::span< T > &recv_buffer) const
Start a non-blocking send of owned data to ranks that ghost the data. The communication is completed ...
Definition: IndexMap.h:265
MPI_Comm comm() const
Return the MPI communicator used to create the index map.
Definition: IndexMap.cpp:641
std::vector< std::int64_t > global_indices() const
Global indices.
Definition: IndexMap.cpp:601
void scatter_rev_begin(const xtl::span< const T > &send_buffer, MPI_Datatype &data_type, MPI_Request &request, const xtl::span< T > &recv_buffer) const
Start a non-blocking send of ghost values to the owning rank. The non-blocking communication is compl...
Definition: IndexMap.h:380
void scatter_rev(xtl::span< T > local_data, const xtl::span< const T > &remote_data, int n, IndexMap::Mode op) const
Send n values for each ghost index to owning to the process.
Definition: IndexMap.h:435
Miscellaneous classes, functions and types.
std::vector< int32_t > compute_owned_indices(const xtl::span< const std::int32_t > &indices, const IndexMap &map)
Given a vector of indices (local numbering, owned or ghost) and an index map, this function returns t...
Definition: IndexMap.cpp:132
std::tuple< std::int64_t, std::vector< std::int32_t >, std::vector< std::vector< std::int64_t > >, std::vector< std::vector< int > > > stack_index_maps(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Compute layout data and ghost indices for a stacked (concatenated) index map, i.e....
Definition: IndexMap.cpp:243