11 #include <dolfinx/common/IndexMap.h>
16 #include <xtl/xspan.hpp>
23 template <
typename T,
class Allocator = std::allocator<T>>
28 Vector(
const std::shared_ptr<const common::IndexMap>&
map,
int bs,
29 const Allocator& alloc = Allocator())
31 _x(
bs * (
map->size_local() +
map->num_ghosts()), alloc)
34 _datatype = MPI::mpi_type<T>();
37 MPI_Type_contiguous(
bs, dolfinx::MPI::mpi_type<T>(), &_datatype);
38 MPI_Type_commit(&_datatype);
44 : _map(x._map), _bs(x._bs), _request(MPI_REQUEST_NULL),
45 _buffer_send_fwd(x._buffer_send_fwd),
46 _buffer_recv_fwd(x._buffer_recv_fwd), _x(x._x)
48 MPI_Type_dup(x._datatype, &_datatype);
53 : _map(std::move(x._map)), _bs(std::move(x._bs)),
54 _datatype(std::exchange(x._datatype, MPI_DATATYPE_NULL)),
55 _request(std::exchange(x._request, MPI_REQUEST_NULL)),
56 _buffer_send_fwd(std::move(x._buffer_send_fwd)),
57 _buffer_recv_fwd(std::move(x._buffer_recv_fwd)), _x(std::move(x._x))
65 MPI_Type_free(&_datatype);
81 const std::vector<std::int32_t>& indices
82 = _map->scatter_fwd_indices().array();
83 _buffer_send_fwd.resize(_bs * indices.size());
84 for (std::size_t i = 0; i < indices.size(); ++i)
86 std::copy_n(std::next(_x.cbegin(), _bs * indices[i]), _bs,
87 std::next(_buffer_send_fwd.begin(), _bs * i));
90 _buffer_recv_fwd.resize(_bs * _map->num_ghosts());
91 _map->scatter_fwd_begin(xtl::span<const T>(_buffer_send_fwd), _datatype,
92 _request, xtl::span<T>(_buffer_recv_fwd));
100 const std::int32_t local_size = _bs * _map->size_local();
101 xtl::span xremote(_x.data() + local_size, _map->num_ghosts() * _bs);
102 _map->scatter_fwd_end(_request);
105 const std::vector<std::int32_t>& scatter_fwd_ghost_pos
106 = _map->scatter_fwd_ghost_positions();
107 for (std::size_t i = 0; i < _map->num_ghosts(); ++i)
109 const int pos = scatter_fwd_ghost_pos[i];
110 std::copy_n(std::next(_buffer_recv_fwd.cbegin(), _bs * pos), _bs,
111 std::next(xremote.begin(), _bs * i));
128 const std::int32_t local_size = _bs * _map->size_local();
129 xtl::span<const T> xremote(_x.data() + local_size,
130 _map->num_ghosts() * _bs);
131 const std::vector<std::int32_t>& scatter_fwd_ghost_pos
132 = _map->scatter_fwd_ghost_positions();
133 _buffer_recv_fwd.resize(_bs * scatter_fwd_ghost_pos.size());
134 for (std::size_t i = 0; i < scatter_fwd_ghost_pos.size(); ++i)
136 const int pos = scatter_fwd_ghost_pos[i];
137 std::copy_n(std::next(xremote.cbegin(), _bs * i), _bs,
138 std::next(_buffer_recv_fwd.begin(), _bs * pos));
142 _buffer_send_fwd.resize(_bs * _map->scatter_fwd_indices().array().size());
143 _map->scatter_rev_begin(xtl::span<const T>(_buffer_recv_fwd), _datatype,
144 _request, xtl::span<T>(_buffer_send_fwd));
156 _map->scatter_rev_end(_request);
159 const std::vector<std::int32_t>& shared_indices
160 = _map->scatter_fwd_indices().array();
163 case common::IndexMap::Mode::insert:
164 for (std::size_t i = 0; i < shared_indices.size(); ++i)
166 std::copy_n(std::next(_buffer_send_fwd.cbegin(), _bs * i), _bs,
167 std::next(_x.begin(), _bs * shared_indices[i]));
170 case common::IndexMap::Mode::add:
171 for (std::size_t i = 0; i < shared_indices.size(); ++i)
172 for (
int j = 0; j < _bs; ++j)
173 _x[shared_indices[i] * _bs + j] += _buffer_send_fwd[i * _bs + j];
200 const std::int32_t size_local = _map->size_local();
201 double local_linf = 0.0;
204 auto local_max_entry = std::max_element(
205 _x.begin(), std::next(_x.begin(), size_local),
206 [](T a, T b) { return std::norm(a) < std::norm(b); });
207 local_linf = std::abs(*local_max_entry);
211 MPI_Allreduce(&local_linf, &linf, 1, MPI_DOUBLE, MPI_MAX,
212 _map->comm(common::IndexMap::Direction::forward));
216 throw std::runtime_error(
"Norm type not supported");
224 const std::int32_t size_local = _map->size_local();
225 double result = std::transform_reduce(
226 _x.begin(), std::next(_x.begin(), size_local), 0.0, std::plus<double>(),
227 [](T val) { return std::norm(val); });
229 MPI_Allreduce(&result, &norm2, 1, MPI_DOUBLE, MPI_SUM,
230 _map->comm(common::IndexMap::Direction::forward));
235 std::shared_ptr<const common::IndexMap>
map()
const {
return _map; }
238 constexpr
int bs()
const {
return _bs; }
241 const std::vector<T, Allocator>&
array()
const {
return _x; }
248 std::shared_ptr<const common::IndexMap> _map;
254 MPI_Datatype _datatype = MPI_DATATYPE_NULL;
255 MPI_Request _request = MPI_REQUEST_NULL;
256 std::vector<T> _buffer_send_fwd, _buffer_recv_fwd;
259 std::vector<T, Allocator> _x;
268 template <
typename T,
class Allocator = std::allocator<T>>
271 const std::int32_t local_size = a.
bs() * a.
map()->size_local();
272 if (local_size != b.
bs() * b.
map()->size_local())
273 throw std::runtime_error(
"Incompatible vector sizes");
274 const std::vector<T>& x_a = a.
array();
275 const std::vector<T>& x_b = b.
array();
277 const T local = std::transform_reduce(
278 x_a.begin(), x_a.begin() + local_size, x_b.begin(),
static_cast<T
>(0),
282 if constexpr (std::is_same<T, std::complex<double>>::value
283 or std::is_same<T, std::complex<float>>::value)
284 return std::conj(a) * b;
290 MPI_Allreduce(&local, &result, 1, dolfinx::MPI::mpi_type<T>(), MPI_SUM,
291 a.map()->comm(common::IndexMap::Direction::forward));