Note: this is documentation for an old release. View the latest documentation at docs.fenicsproject.org/dolfinx/v0.9.0/cpp/doxygen/d6/dbe/Vector_8h_source.html
DOLFINx  0.5.1
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 
19 namespace dolfinx::la
20 {
21 
23 
24 template <typename T, class Allocator = std::allocator<T>>
25 class Vector
26 {
27 public:
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 
111  void scatter_fwd()
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 
184 private:
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 
210 template <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 
240 template <typename T, class Allocator>
242 {
243  T result = inner_product(a, a);
244  return std::real(result);
245 }
246 
251 template <typename T, class Allocator>
252 auto 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 
281 template <typename T, typename U>
282 void 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 
315 template <typename T, typename U>
316 bool 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::shared_ptr< const common::IndexMap > map() const
Get IndexMap.
Definition: Vector.h:170
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
std::span< const T > array() const
Get local part of the vector (const version)
Definition: Vector.h:176
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
void set(T v)
Set all entries (including ghosts)
Definition: Vector.h:70
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
std::span< T > mutable_array()
Get local part of the vector.
Definition: Vector.h:179
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
auto norm(const Vector< T, Allocator > &a, Norm type=Norm::l2)
Compute the norm of the vector.
Definition: Vector.h:252
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
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