DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
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 <algorithm>
11#include <cmath>
12#include <complex>
13#include <dolfinx/common/IndexMap.h>
14#include <dolfinx/common/Scatterer.h>
15#include <dolfinx/common/types.h>
16#include <limits>
17#include <memory>
18#include <numeric>
19#include <span>
20#include <type_traits>
21#include <vector>
22
23namespace dolfinx::la
24{
25
30template <typename T, typename Container = std::vector<T>>
31class Vector
32{
33 static_assert(std::is_same_v<typename Container::value_type, T>);
34
35public:
37 using value_type = T;
38
40 using container_type = Container;
41
42 static_assert(std::is_same_v<value_type, typename container_type::value_type>,
43 "Scalar type and container value type must be the same.");
44
48 Vector(std::shared_ptr<const common::IndexMap> map, int bs)
49 : _map(map), _scatterer(std::make_shared<common::Scatterer<>>(*_map, bs)),
50 _bs(bs), _buffer_local(_scatterer->local_buffer_size()),
51 _buffer_remote(_scatterer->remote_buffer_size()),
52 _x(bs * (map->size_local() + map->num_ghosts()))
53 {
54 }
55
57 Vector(const Vector& x)
58 : _map(x._map), _scatterer(x._scatterer), _bs(x._bs),
59 _request(1, MPI_REQUEST_NULL), _buffer_local(x._buffer_local),
60 _buffer_remote(x._buffer_remote), _x(x._x)
61 {
62 }
63
66 : _map(std::move(x._map)), _scatterer(std::move(x._scatterer)),
67 _bs(std::move(x._bs)),
68 _request(std::exchange(x._request, {MPI_REQUEST_NULL})),
69 _buffer_local(std::move(x._buffer_local)),
70 _buffer_remote(std::move(x._buffer_remote)), _x(std::move(x._x))
71 {
72 }
73
74 // Assignment operator (disabled)
75 Vector& operator=(const Vector& x) = delete;
76
78 Vector& operator=(Vector&& x) = default;
79
82 void set(value_type v) { std::ranges::fill(_x, v); }
83
87 {
88 const std::int32_t local_size = _bs * _map->size_local();
89 std::span<const value_type> x_local(_x.data(), local_size);
90
91 auto pack = [](auto&& in, auto&& idx, auto&& out)
92 {
93 for (std::size_t i = 0; i < idx.size(); ++i)
94 out[i] = in[idx[i]];
95 };
96 pack(x_local, _scatterer->local_indices(), _buffer_local);
97
98 _scatterer->scatter_fwd_begin(std::span<const value_type>(_buffer_local),
99 std::span<value_type>(_buffer_remote),
100 std::span<MPI_Request>(_request));
101 }
102
106 {
107 const std::int32_t local_size = _bs * _map->size_local();
108 const std::int32_t num_ghosts = _bs * _map->num_ghosts();
109 std::span<value_type> x_remote(_x.data() + local_size, num_ghosts);
110 _scatterer->scatter_fwd_end(std::span<MPI_Request>(_request));
111
112 auto unpack = [](auto&& in, auto&& idx, auto&& out, auto op)
113 {
114 for (std::size_t i = 0; i < idx.size(); ++i)
115 out[idx[i]] = op(out[idx[i]], in[i]);
116 };
117
118 unpack(_buffer_remote, _scatterer->remote_indices(), x_remote,
119 [](auto /*a*/, auto b) { return b; });
120 }
121
125 {
126 this->scatter_fwd_begin();
127 this->scatter_fwd_end();
128 }
129
133 {
134 const std::int32_t local_size = _bs * _map->size_local();
135 const std::int32_t num_ghosts = _bs * _map->num_ghosts();
136 std::span<value_type> x_remote(_x.data() + local_size, num_ghosts);
137
138 auto pack = [](auto&& in, auto&& idx, auto&& out)
139 {
140 for (std::size_t i = 0; i < idx.size(); ++i)
141 out[i] = in[idx[i]];
142 };
143 pack(x_remote, _scatterer->remote_indices(), _buffer_remote);
144
145 _scatterer->scatter_rev_begin(std::span<const value_type>(_buffer_remote),
146 std::span<value_type>(_buffer_local),
147 _request);
148 }
149
156 template <class BinaryOperation>
157 void scatter_rev_end(BinaryOperation op)
158 {
159 const std::int32_t local_size = _bs * _map->size_local();
160 std::span<value_type> x_local(_x.data(), local_size);
161 _scatterer->scatter_rev_end(_request);
162
163 auto unpack = [](auto&& in, auto&& idx, auto&& out, auto op)
164 {
165 for (std::size_t i = 0; i < idx.size(); ++i)
166 out[idx[i]] = op(out[idx[i]], in[i]);
167 };
168 unpack(_buffer_local, _scatterer->local_indices(), x_local, op);
169 }
170
176 template <class BinaryOperation>
177 void scatter_rev(BinaryOperation op)
178 {
179 this->scatter_rev_begin();
180 this->scatter_rev_end(op);
181 }
182
184 std::shared_ptr<const common::IndexMap> index_map() const { return _map; }
185
187 constexpr int bs() const { return _bs; }
188
190 std::span<const value_type> array() const
191 {
192 return std::span<const value_type>(_x);
193 }
194
196 std::span<value_type> mutable_array() { return std::span(_x); }
197
198private:
199 // Map describing the data layout
200 std::shared_ptr<const common::IndexMap> _map;
201
202 // Scatter for managing MPI communication
203 std::shared_ptr<const common::Scatterer<>> _scatterer;
204
205 // Block size
206 int _bs;
207
208 // MPI request handle
209 std::vector<MPI_Request> _request = {MPI_REQUEST_NULL};
210
211 // Buffers for ghost scatters
212 container_type _buffer_local, _buffer_remote;
213
214 // Vector data
216};
217
224template <class V>
225auto inner_product(const V& a, const V& b)
226{
227 using T = typename V::value_type;
228 const std::int32_t local_size = a.bs() * a.index_map()->size_local();
229 if (local_size != b.bs() * b.index_map()->size_local())
230 throw std::runtime_error("Incompatible vector sizes");
231 std::span<const T> x_a = a.array().subspan(0, local_size);
232 std::span<const T> x_b = b.array().subspan(0, local_size);
233
234 const T local = std::transform_reduce(
235 x_a.begin(), x_a.end(), x_b.begin(), static_cast<T>(0), std::plus{},
236 [](T a, T b) -> T
237 {
238 if constexpr (std::is_same<T, std::complex<double>>::value
239 or std::is_same<T, std::complex<float>>::value)
240 {
241 return std::conj(a) * b;
242 }
243 else
244 return a * b;
245 });
246
247 T result;
248 MPI_Allreduce(&local, &result, 1, dolfinx::MPI::mpi_type<T>(), MPI_SUM,
249 a.index_map()->comm());
250 return result;
251}
252
255template <class V>
256auto squared_norm(const V& a)
257{
258 using T = typename V::value_type;
259 T result = inner_product(a, a);
260 return std::real(result);
261}
262
267template <class V>
268auto norm(const V& x, Norm type = Norm::l2)
269{
270 using T = typename V::value_type;
271 switch (type)
272 {
273 case Norm::l1:
274 {
275 std::int32_t size_local = x.bs() * x.index_map()->size_local();
276 std::span<const T> data = x.array().subspan(0, size_local);
277 using U = typename dolfinx::scalar_value_type_t<T>;
278 U local_l1
279 = std::accumulate(data.begin(), data.end(), U(0),
280 [](auto norm, auto x) { return norm + std::abs(x); });
281 U l1(0);
282 MPI_Allreduce(&local_l1, &l1, 1, MPI::mpi_type<U>(), MPI_SUM,
283 x.index_map()->comm());
284 return l1;
285 }
286 case Norm::l2:
287 return std::sqrt(squared_norm(x));
288 case Norm::linf:
289 {
290 std::int32_t size_local = x.bs() * x.index_map()->size_local();
291 std::span<const T> data = x.array().subspan(0, size_local);
292 auto max_pos = std::ranges::max_element(
293 data, [](T a, T b) { return std::norm(a) < std::norm(b); });
294 auto local_linf = std::abs(*max_pos);
295 decltype(local_linf) linf = 0;
296 MPI_Allreduce(&local_linf, &linf, 1, MPI::mpi_type<decltype(linf)>(),
297 MPI_MAX, x.index_map()->comm());
298 return linf;
299 }
300 default:
301 throw std::runtime_error("Norm type not supported");
302 }
303}
304
310template <class V>
311void orthonormalize(std::vector<std::reference_wrapper<V>> basis)
312{
313 using T = typename V::value_type;
314 using U = typename dolfinx::scalar_value_type_t<T>;
315
316 // Loop over each vector in basis
317 for (std::size_t i = 0; i < basis.size(); ++i)
318 {
319 // Orthogonalize vector i with respect to previously orthonormalized
320 // vectors
321 V& bi = basis[i].get();
322 for (std::size_t j = 0; j < i; ++j)
323 {
324 const V& bj = basis[j].get();
325
326 // basis_i <- basis_i - dot_ij basis_j
327 auto dot_ij = inner_product(bi, bj);
328 std::ranges::transform(bj.array(), bi.array(), bi.mutable_array().begin(),
329 [dot_ij](auto xj, auto xi)
330 { return xi - dot_ij * xj; });
331 }
332
333 // Normalise basis function
334 auto norm = la::norm(bi, la::Norm::l2);
335 if (norm * norm < std::numeric_limits<U>::epsilon())
336 {
337 throw std::runtime_error(
338 "Linear dependency detected. Cannot orthogonalize.");
339 }
340 std::ranges::transform(bi.array(), bi.mutable_array().begin(),
341 [norm](auto x) { return x / norm; });
342 }
343}
344
353template <class V>
355 std::vector<std::reference_wrapper<const V>> basis,
356 dolfinx::scalar_value_type_t<typename V::value_type> eps
357 = std::numeric_limits<
358 dolfinx::scalar_value_type_t<typename V::value_type>>::epsilon())
359{
360 using T = typename V::value_type;
361 for (std::size_t i = 0; i < basis.size(); i++)
362 {
363 for (std::size_t j = i; j < basis.size(); j++)
364 {
365 T delta_ij = (i == j) ? T(1) : T(0);
366 auto dot_ij = inner_product(basis[i].get(), basis[j].get());
367 if (std::norm(delta_ij - dot_ij) > eps)
368 return false;
369 }
370 }
371
372 return true;
373}
374
375} // namespace dolfinx::la
Definition Vector.h:32
void scatter_fwd_begin()
Definition Vector.h:86
constexpr int bs() const
Get block size.
Definition Vector.h:187
Vector(std::shared_ptr< const common::IndexMap > map, int bs)
Definition Vector.h:48
Vector & operator=(Vector &&x)=default
Move Assignment operator.
void scatter_fwd()
Definition Vector.h:124
std::span< const value_type > array() const
Get local part of the vector (const version)
Definition Vector.h:190
std::shared_ptr< const common::IndexMap > index_map() const
Get IndexMap.
Definition Vector.h:184
Container container_type
Container type.
Definition Vector.h:40
Vector(Vector &&x)
Move constructor.
Definition Vector.h:65
Vector(const Vector &x)
Copy constructor.
Definition Vector.h:57
void scatter_fwd_end()
Definition Vector.h:105
void scatter_rev_begin()
Definition Vector.h:132
void set(value_type v)
Definition Vector.h:82
std::span< value_type > mutable_array()
Get local part of the vector.
Definition Vector.h:196
void scatter_rev_end(BinaryOperation op)
Definition Vector.h:157
void scatter_rev(BinaryOperation op)
Definition Vector.h:177
T value_type
Scalar type.
Definition Vector.h:37
constexpr MPI_Datatype mpi_type()
MPI Type.
Definition MPI.h:273
Linear algebra interface.
Definition sparsitybuild.h:15
void orthonormalize(std::vector< std::reference_wrapper< V > > basis)
Definition Vector.h:311
Norm
Norm types.
Definition utils.h:17
auto squared_norm(const V &a)
Definition Vector.h:256
auto norm(const V &x, Norm type=Norm::l2)
Definition Vector.h:268
auto inner_product(const V &a, const V &b)
Definition Vector.h:225
bool is_orthonormal(std::vector< std::reference_wrapper< const V > > basis, dolfinx::scalar_value_type_t< typename V::value_type > eps=std::numeric_limits< dolfinx::scalar_value_type_t< typename V::value_type > >::epsilon())
Test if basis is orthonormal.
Definition Vector.h:354