DOLFINx
DOLFINx C++ interface
Vector.h
1// Copyright (C) 2020 Garth N. Wells
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 "utils.h"
10#include <complex>
11#include <dolfinx/common/IndexMap.h>
12#include <dolfinx/common/Scatterer.h>
13#include <limits>
14#include <memory>
15#include <numeric>
16#include <span>
17#include <vector>
18
19namespace dolfinx::la
20{
21
23
24template <typename T, class Allocator = std::allocator<T>>
25class Vector
26{
27public:
29 using value_type = T;
30
32 using allocator_type = Allocator;
33
35 Vector(const std::shared_ptr<const common::IndexMap>& map, int bs,
36 const Allocator& alloc = Allocator())
37 : _map(map), _scatterer(std::make_shared<common::Scatterer>(*_map, bs)),
38 _bs(bs), _buffer_local(_scatterer->local_buffer_size(), alloc),
39 _buffer_remote(_scatterer->remote_buffer_size(), alloc),
40 _x(bs * (map->size_local() + map->num_ghosts()), alloc)
41 {
42 }
43
45 Vector(const Vector& x)
46 : _map(x._map), _scatterer(x._scatterer), _bs(x._bs),
47 _request(MPI_REQUEST_NULL), _buffer_local(x._buffer_local),
48 _buffer_remote(x._buffer_remote), _x(x._x)
49 {
50 }
51
54 : _map(std::move(x._map)), _scatterer(std::move(x._scatterer)),
55 _bs(std::move(x._bs)),
56 _request(std::exchange(x._request, MPI_REQUEST_NULL)),
57 _buffer_local(std::move(x._buffer_local)),
58 _buffer_remote(std::move(x._buffer_remote)), _x(std::move(x._x))
59 {
60 }
61
62 // Assignment operator (disabled)
63 Vector& operator=(const Vector& x) = delete;
64
66 Vector& operator=(Vector&& x) = default;
67
70 void set(T v) { std::fill(_x.begin(), _x.end(), v); }
71
75 {
76 const std::int32_t local_size = _bs * _map->size_local();
77 std::span<const T> x_local(_x.data(), local_size);
78
79 auto pack = [](const auto& in, const auto& idx, auto& out)
80 {
81 for (std::size_t i = 0; i < idx.size(); ++i)
82 out[i] = in[idx[i]];
83 };
84 pack(x_local, _scatterer->local_indices(), _buffer_local);
85
86 _scatterer->scatter_fwd_begin(std::span<const T>(_buffer_local),
87 std::span<T>(_buffer_remote), _request);
88 }
89
93 {
94 const std::int32_t local_size = _bs * _map->size_local();
95 const std::int32_t num_ghosts = _bs * _map->num_ghosts();
96 std::span<T> x_remote(_x.data() + local_size, num_ghosts);
97 _scatterer->scatter_fwd_end(_request);
98
99 auto unpack = [](const auto& in, const auto& idx, auto& out, auto op)
100 {
101 for (std::size_t i = 0; i < idx.size(); ++i)
102 out[idx[i]] = op(out[idx[i]], in[i]);
103 };
104
105 unpack(_buffer_remote, _scatterer->remote_indices(), x_remote,
106 [](auto /*a*/, auto b) { return b; });
107 }
108
112 {
113 this->scatter_fwd_begin();
114 this->scatter_fwd_end();
115 }
116
120 {
121 const std::int32_t local_size = _bs * _map->size_local();
122 const std::int32_t num_ghosts = _bs * _map->num_ghosts();
123 std::span<T> x_remote(_x.data() + local_size, num_ghosts);
124
125 auto pack = [](const auto& in, const auto& idx, auto& out)
126 {
127 for (std::size_t i = 0; i < idx.size(); ++i)
128 out[i] = in[idx[i]];
129 };
130 pack(x_remote, _scatterer->remote_indices(), _buffer_remote);
131
132 _scatterer->scatter_rev_begin(std::span<const T>(_buffer_remote),
133 std::span<T>(_buffer_local), _request);
134 }
135
142 template <class BinaryOperation>
143 void scatter_rev_end(BinaryOperation op)
144 {
145 const std::int32_t local_size = _bs * _map->size_local();
146 std::span<T> x_local(_x.data(), local_size);
147 _scatterer->scatter_rev_end(_request);
148
149 auto unpack = [](const auto& in, const auto& idx, auto& out, auto op)
150 {
151 for (std::size_t i = 0; i < idx.size(); ++i)
152 out[idx[i]] = op(out[idx[i]], in[i]);
153 };
154 unpack(_buffer_local, _scatterer->local_indices(), x_local, op);
155 }
156
162 template <class BinaryOperation>
163 void scatter_rev(BinaryOperation op)
164 {
165 this->scatter_rev_begin();
166 this->scatter_rev_end(op);
167 }
168
170 std::shared_ptr<const common::IndexMap> map() const { return _map; }
171
173 constexpr int bs() const { return _bs; }
174
176 std::span<const T> array() const { return std::span<const T>(_x); }
177
179 std::span<T> mutable_array() { return std::span(_x); }
180
182 constexpr allocator_type allocator() const { return _x.get_allocator(); }
183
184private:
185 // Map describing the data layout
186 std::shared_ptr<const common::IndexMap> _map;
187
188 // Scatter for managing MPI communication
189 std::shared_ptr<const common::Scatterer> _scatterer;
190
191 // Block size
192 int _bs;
193
194 // MPI request handle
195 MPI_Request _request = MPI_REQUEST_NULL;
196
197 // Buffers for ghost scatters
198 std::vector<T, Allocator> _buffer_local, _buffer_remote;
199
200 // Vector data
201 std::vector<T, Allocator> _x;
202};
203
210template <typename T, class Allocator>
212{
213 const std::int32_t local_size = a.bs() * a.map()->size_local();
214 if (local_size != b.bs() * b.map()->size_local())
215 throw std::runtime_error("Incompatible vector sizes");
216 std::span<const T> x_a = a.array().subspan(0, local_size);
217 std::span<const T> x_b = b.array().subspan(0, local_size);
218
219 const T local = std::transform_reduce(
220 x_a.begin(), x_a.end(), x_b.begin(), static_cast<T>(0), std::plus{},
221 [](T a, T b) -> T
222 {
223 if constexpr (std::is_same<T, std::complex<double>>::value
224 or std::is_same<T, std::complex<float>>::value)
225 {
226 return std::conj(a) * b;
227 }
228 else
229 return a * b;
230 });
231
232 T result;
233 MPI_Allreduce(&local, &result, 1, dolfinx::MPI::mpi_type<T>(), MPI_SUM,
234 a.map()->comm());
235 return result;
236}
237
240template <typename T, class Allocator>
242{
243 T result = inner_product(a, a);
244 return std::real(result);
245}
246
251template <typename T, class Allocator>
252auto norm(const Vector<T, Allocator>& a, Norm type = Norm::l2)
253{
254 switch (type)
255 {
256 case Norm::l2:
257 return std::sqrt(squared_norm(a));
258 case Norm::linf:
259 {
260 const std::int32_t size_local = a.bs() * a.map()->size_local();
261 std::span<const T> x_a = a.array().subspan(0, size_local);
262 auto max_pos = std::max_element(x_a.begin(), x_a.end(),
263 [](T a, T b)
264 { return std::norm(a) < std::norm(b); });
265 auto local_linf = std::abs(*max_pos);
266 decltype(local_linf) linf = 0;
267 MPI_Allreduce(&local_linf, &linf, 1, MPI::mpi_type<decltype(linf)>(),
268 MPI_MAX, a.map()->comm());
269 return linf;
270 }
271 default:
272 throw std::runtime_error("Norm type not supported");
273 }
274}
275
281template <typename T, typename U>
282void orthonormalize(const std::span<Vector<T, U>>& basis, double tol = 1.0e-10)
283{
284 // Loop over each vector in basis
285 for (std::size_t i = 0; i < basis.size(); ++i)
286 {
287 // Orthogonalize vector i with respect to previously orthonormalized
288 // vectors
289 for (std::size_t j = 0; j < i; ++j)
290 {
291 // basis_i <- basis_i - dot_ij basis_j
292 T dot_ij = inner_product(basis[i], basis[j]);
293 std::transform(basis[j].array().begin(), basis[j].array().end(),
294 basis[i].array().begin(), basis[i].mutable_array().begin(),
295 [dot_ij](auto xj, auto xi) { return xi - dot_ij * xj; });
296 }
297
298 // Normalise basis function
299 double norm = la::norm(basis[i], la::Norm::l2);
300 if (norm < tol)
301 {
302 throw std::runtime_error(
303 "Linear dependency detected. Cannot orthogonalize.");
304 }
305 std::transform(basis[i].array().begin(), basis[i].array().end(),
306 basis[i].mutable_array().begin(),
307 [norm](auto x) { return x / norm; });
308 }
309}
310
315template <typename T, typename U>
316bool is_orthonormal(const std::span<const Vector<T, U>>& basis,
317 double tol = 1.0e-10)
318{
319 for (std::size_t i = 0; i < basis.size(); i++)
320 {
321 for (std::size_t j = i; j < basis.size(); j++)
322 {
323 const double delta_ij = (i == j) ? 1.0 : 0.0;
324 T dot_ij = inner_product(basis[i], basis[j]);
325 if (std::abs(delta_ij - dot_ij) > tol)
326 return false;
327 }
328 }
329
330 return true;
331}
332
333} // namespace dolfinx::la
Distributed vector.
Definition: Vector.h:26
std::span< const T > array() const
Get local part of the vector (const version)
Definition: Vector.h:176
void scatter_fwd_begin()
Begin scatter of local data from owner to ghosts on other ranks.
Definition: Vector.h:74
Allocator allocator_type
The allocator type.
Definition: Vector.h:32
Vector(const std::shared_ptr< const common::IndexMap > &map, int bs, const Allocator &alloc=Allocator())
Create a distributed vector.
Definition: Vector.h:35
constexpr int bs() const
Get block size.
Definition: Vector.h:173
std::shared_ptr< const common::IndexMap > map() const
Get IndexMap.
Definition: Vector.h:170
void set(T v)
Set all entries (including ghosts)
Definition: Vector.h:70
std::span< T > mutable_array()
Get local part of the vector.
Definition: Vector.h:179
Vector & operator=(Vector &&x)=default
Move Assignment operator.
void scatter_fwd()
Scatter local data to ghost positions on other ranks.
Definition: Vector.h:111
Vector(Vector &&x)
Move constructor.
Definition: Vector.h:53
Vector(const Vector &x)
Copy constructor.
Definition: Vector.h:45
void scatter_fwd_end()
End scatter of local data from owner to ghosts on other ranks.
Definition: Vector.h:92
void scatter_rev_begin()
Start scatter of ghost data to owner.
Definition: Vector.h:119
constexpr allocator_type allocator() const
Get the allocator associated with the container.
Definition: Vector.h:182
void scatter_rev_end(BinaryOperation op)
End scatter of ghost data to owner. This process may receive data from more than one process,...
Definition: Vector.h:143
void scatter_rev(BinaryOperation op)
Scatter ghost data to owner. This process may receive data from more than one process,...
Definition: Vector.h:163
T value_type
The value type.
Definition: Vector.h:29
Linear algebra interface.
Definition: sparsitybuild.h:15
Norm
Norm types.
Definition: utils.h:13
bool is_orthonormal(const std::span< const Vector< T, U > > &basis, double tol=1.0e-10)
Test if basis is orthonormal.
Definition: Vector.h:316
T inner_product(const Vector< T, Allocator > &a, const Vector< T, Allocator > &b)
Compute the inner product of two vectors. The two vectors must have the same parallel layout.
Definition: Vector.h:211
void orthonormalize(const std::span< Vector< T, U > > &basis, double tol=1.0e-10)
Orthonormalize a set of vectors.
Definition: Vector.h:282
auto squared_norm(const Vector< T, Allocator > &a)
Compute the squared L2 norm of vector.
Definition: Vector.h:241
auto norm(const Vector< T, Allocator > &a, Norm type=Norm::l2)
Compute the norm of the vector.
Definition: Vector.h:252