Loading...
Searching...
No Matches
Vector.h
// Copyright (C) 2020-2025 Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#pragma once
#include "utils.h"
#include <algorithm>
#include <cmath>
#include <complex>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/common/Scatterer.h>
#include <dolfinx/common/types.h>
#include <limits>
#include <memory>
#include <numeric>
#include <span>
#include <type_traits>
#include <vector>
namespace dolfinx::la
{
template <class F, class Container, class ScatterContainer>
concept VectorPackKernel = requires(F f, ScatterContainer idx, Container x) {
f(idx.cbegin(), idx.cend(), x.cbegin(), x.begin());
};
template <class GetPtr, class U, class T>
concept GetPtrConcept = requires(GetPtr f, U x) {
{ f(x) } -> std::same_as<T*>;
} and requires(GetPtr f, const U x) {
{ f(x) } -> std::same_as<const T*>;
};
template <typename T, typename Container = std::vector<T>,
typename ScatterContainer = std::vector<std::int32_t>>
class Vector
{
static_assert(std::is_same_v<typename Container::value_type, T>);
template <typename, typename, typename>
friend class Vector;
private:
auto get_pack()
{
return [](typename ScatterContainer::const_iterator idx_first,
typename ScatterContainer::const_iterator idx_last,
const auto in_first, auto out_first)
{
// out[i] = in[idx[i]]
std::transform(idx_first, idx_last, out_first,
[in_first](auto p) { return *std::next(in_first, p); });
};
}
auto get_unpack()
{
return [](typename ScatterContainer::const_iterator idx_first,
typename ScatterContainer::const_iterator idx_last,
const auto in_first, auto out_first)
{
// out[idx[i]] = in[i]
for (typename ScatterContainer::const_iterator idx = idx_first;
idx != idx_last; ++idx)
{
std::size_t d = std::distance(idx_first, idx);
*std::next(out_first, *idx) = *std::next(in_first, d);
}
};
}
template <typename BinaryOp>
auto get_unpack_op(BinaryOp op)
{
return [op](typename ScatterContainer::const_iterator idx_first,
typename ScatterContainer::const_iterator idx_last,
const auto in_first, auto out_first)
{
// out[idx[i]] = op(out[idx[i]], in[i])
for (typename ScatterContainer::const_iterator idx = idx_first;
idx != idx_last; ++idx)
{
std::size_t d = std::distance(idx_first, idx);
auto& out = *std::next(out_first, *idx);
out = op(out, *std::next(in_first, d));
}
};
}
public:
using container_type = Container;
using value_type = container_type::value_type;
static_assert(std::is_same_v<value_type, typename container_type::value_type>,
"Scalar type and container value type must be the same.");
Vector(std::shared_ptr<const common::IndexMap> map, int bs)
: _map(map), _bs(bs), _x(bs * (map->size_local() + map->num_ghosts())),
_scatterer(
std::make_shared<common::Scatterer<ScatterContainer>>(*_map, bs)),
_buffer_local(_scatterer->local_indices().size()),
_buffer_remote(_scatterer->remote_indices().size())
{
}
Vector(const Vector& x) = default;
Vector(Vector&& x) = default;
private:
std::shared_ptr<const common::Scatterer<ScatterContainer>>
scatter_ptr(auto sc) const
{
using SC = typename std::remove_cv<
typename decltype(sc)::element_type>::type::container_type;
if constexpr (std::is_same_v<ScatterContainer, SC>)
return sc; // Scatters use same container
else // Scatters use different containers, so copy
return std::make_shared<common::Scatterer<ScatterContainer>>(*sc);
}
public:
template <typename T0, typename Container0, typename ScatterContainer0>
explicit Vector(const Vector<T0, Container0, ScatterContainer0>& x)
: _map(x.index_map()), _bs(x.bs()), _x(x._x.begin(), x._x.end()),
_scatterer(scatter_ptr(x._scatterer)), _request(MPI_REQUEST_NULL),
_buffer_local(_scatterer->local_indices().size()),
_buffer_remote(_scatterer->remote_indices().size())
{
}
// Assignment operator (disabled)
Vector& operator=(const Vector& x) = delete;
Vector& operator=(Vector&& x) = default;
[[deprecated("Use std::ranges::fill(u.array(), v) instead.")]] void
set(value_type v)
{
std::ranges::fill(_x, v);
}
template <typename U, typename GetPtr>
requires VectorPackKernel<U, container_type, ScatterContainer>
and GetPtrConcept<GetPtr, Container, T>
void scatter_fwd_begin(U pack, GetPtr get_ptr)
{
pack(_scatterer->local_indices().begin(), _scatterer->local_indices().end(),
_x.begin(), _buffer_local.begin());
_scatterer->scatter_fwd_begin(get_ptr(_buffer_local),
get_ptr(_buffer_remote), _request);
}
void scatter_fwd_begin()
requires requires(Container c) {
{ c.data() } -> std::same_as<T*>;
}
{
scatter_fwd_begin(get_pack(), [](auto&& x) { return x.data(); });
}
template <typename U>
requires VectorPackKernel<U, container_type, ScatterContainer>
void scatter_fwd_end(U unpack)
{
_scatterer->scatter_end(_request);
unpack(_scatterer->remote_indices().begin(),
_scatterer->remote_indices().end(), _buffer_remote.begin(),
std::next(_x.begin(), _bs * _map->size_local()));
}
void scatter_fwd_end()
requires requires(Container c) {
{ c.data() } -> std::same_as<T*>;
}
{
this->scatter_fwd_end(get_unpack());
}
void scatter_fwd()
requires requires(Container c) {
{ c.data() } -> std::same_as<T*>;
}
{
this->scatter_fwd_begin(get_pack(), [](auto&& x) { return x.data(); });
this->scatter_fwd_end(get_unpack());
}
template <typename U, typename GetPtr>
requires VectorPackKernel<U, container_type, ScatterContainer>
and GetPtrConcept<GetPtr, Container, T>
void scatter_rev_begin(U pack, GetPtr get_ptr)
{
std::int32_t local_size = _bs * _map->size_local();
pack(_scatterer->remote_indices().begin(),
_scatterer->remote_indices().end(), std::next(_x.begin(), local_size),
_buffer_remote.begin());
_scatterer->scatter_rev_begin(get_ptr(_buffer_remote),
get_ptr(_buffer_local), _request);
}
void scatter_rev_begin()
requires requires(Container c) {
{ c.data() } -> std::same_as<T*>;
}
{
scatter_rev_begin(get_pack(), [](auto&& x) { return x.data(); });
}
template <typename U>
requires VectorPackKernel<U, container_type, ScatterContainer>
void scatter_rev_end(U unpack)
{
_scatterer->scatter_end(_request);
unpack(_scatterer->local_indices().begin(),
_scatterer->local_indices().end(), _buffer_local.begin(),
_x.begin());
}
template <class BinaryOperation>
requires requires(Container c) {
{ c.data() } -> std::same_as<T*>;
}
void scatter_rev(BinaryOperation op)
{
this->scatter_rev_begin();
this->scatter_rev_end(get_unpack_op(op));
}
std::shared_ptr<const common::IndexMap> index_map() const { return _map; }
constexpr int bs() const { return _bs; }
container_type& array() { return _x; }
const container_type& array() const { return _x; }
[[deprecated("Use array() instead.")]] container_type& mutable_array()
{
return _x;
}
private:
// Map describing the data layout
std::shared_ptr<const common::IndexMap> _map;
// Block size
int _bs;
// Vector data
container_type _x;
// Scatter for managing MPI communication
std::shared_ptr<const common::Scatterer<ScatterContainer>> _scatterer;
// MPI request handle
MPI_Request _request = MPI_REQUEST_NULL;
// Buffers for ghost scatters
container_type _buffer_local, _buffer_remote;
}; // namespace dolfinx::la
template <class V>
auto inner_product(const V& a, const V& b)
{
using T = typename V::value_type;
const std::int32_t local_size = a.bs() * a.index_map()->size_local();
if (local_size != b.bs() * b.index_map()->size_local())
throw std::runtime_error("Incompatible vector sizes");
const T local = std::transform_reduce(
a.array().begin(), std::next(a.array().begin(), local_size),
b.array().begin(), static_cast<T>(0), std::plus{},
[](T a, T b) -> T
{
if constexpr (std::is_same<T, std::complex<double>>::value
or std::is_same<T, std::complex<float>>::value)
{
return std::conj(a) * b;
}
else
return a * b;
});
T result;
MPI_Allreduce(&local, &result, 1, dolfinx::MPI::mpi_t<T>, MPI_SUM,
a.index_map()->comm());
return result;
}
template <class V>
auto squared_norm(const V& a)
{
using T = typename V::value_type;
T result = inner_product(a, a);
return std::real(result);
}
template <class V>
auto norm(const V& x, Norm type = Norm::l2)
{
using T = typename V::value_type;
switch (type)
{
case Norm::l1:
{
std::int32_t size_local = x.bs() * x.index_map()->size_local();
using U = typename dolfinx::scalar_value_t<T>;
U local_l1 = std::accumulate(
x.array().begin(), std::next(x.array().begin(), size_local), U(0),
[](auto norm, auto x) { return norm + std::abs(x); });
U l1(0);
MPI_Allreduce(&local_l1, &l1, 1, MPI::mpi_t<U>, MPI_SUM,
x.index_map()->comm());
return l1;
}
case Norm::l2:
return std::sqrt(squared_norm(x));
case Norm::linf:
{
std::int32_t size_local = x.bs() * x.index_map()->size_local();
auto max_pos = std::max_element(
x.array().begin(), std::next(x.array().begin(), size_local),
[](T a, T b) { return std::norm(a) < std::norm(b); });
auto local_linf = std::abs(*max_pos);
decltype(local_linf) linf = 0;
MPI_Allreduce(&local_linf, &linf, 1, MPI::mpi_t<decltype(linf)>, MPI_MAX,
x.index_map()->comm());
return linf;
}
default:
throw std::runtime_error("Norm type not supported");
}
}
template <class V>
void orthonormalize(std::vector<std::reference_wrapper<V>> basis)
{
using T = typename V::value_type;
using U = typename dolfinx::scalar_value_t<T>;
// Loop over each vector in basis
for (std::size_t i = 0; i < basis.size(); ++i)
{
// Orthogonalize vector i with respect to previously orthonormalized
// vectors
V& bi = basis[i].get();
for (std::size_t j = 0; j < i; ++j)
{
const V& bj = basis[j].get();
// basis_i <- basis_i - dot_ij basis_j
auto dot_ij = inner_product(bi, bj);
std::ranges::transform(bj.array(), bi.array(), bi.array().begin(),
[dot_ij](auto xj, auto xi)
{ return xi - dot_ij * xj; });
}
// Normalise basis function
auto norm = la::norm(bi, la::Norm::l2);
if (norm * norm < std::numeric_limits<U>::epsilon())
{
throw std::runtime_error(
"Linear dependency detected. Cannot orthogonalize.");
}
std::ranges::transform(bi.array(), bi.array().begin(),
[norm](auto x) { return x / norm; });
}
}
template <class V>
bool is_orthonormal(
std::vector<std::reference_wrapper<const V>> basis,
dolfinx::scalar_value_t<typename V::value_type> eps = std::numeric_limits<
dolfinx::scalar_value_t<typename V::value_type>>::epsilon())
{
using T = typename V::value_type;
for (std::size_t i = 0; i < basis.size(); i++)
{
for (std::size_t j = i; j < basis.size(); ++j)
{
T delta_ij = (i == j) ? T(1) : T(0);
auto dot_ij = inner_product(basis[i].get(), basis[j].get());
if (std::norm(delta_ij - dot_ij) > eps)
return false;
}
}
return true;
}
} // namespace dolfinx::la
dolfinx::la::Vector
A vector that can be distributed across processes.
Definition Vector.h:50
dolfinx::la::Vector::mutable_array
container_type & mutable_array()
Get local part of the vector.
Definition Vector.h:396
dolfinx::la::Vector::array
container_type & array()
Get the process-local part of the vector.
Definition Vector.h:387
dolfinx::la::Vector::value_type
container_type::value_type value_type
Scalar type.
Definition Vector.h:121
dolfinx::la::Vector::array
const container_type & array() const
Get the process-local part of the vector (const version).
Definition Vector.h:392
dolfinx::la::Vector::scatter_rev_begin
void scatter_rev_begin(U pack, GetPtr get_ptr)
Start scatter (send) of ghost entry data to the owning process of an index.
Definition Vector.h:314
dolfinx::la::Vector::Vector
Vector(const Vector< T0, Container0, ScatterContainer0 > &x)
Copy-convert vector, possibly using to different container types.
Definition Vector.h:178
dolfinx::la::Vector::bs
constexpr int bs() const
Get block size.
Definition Vector.h:382
dolfinx::la::Vector::Vector
Vector(std::shared_ptr< const common::IndexMap > map, int bs)
Create a distributed vector.
Definition Vector.h:131
dolfinx::la::Vector::scatter_rev
void scatter_rev(BinaryOperation op)
Scatter (send) of ghost data values to the owning process and assign/accumulate into the owned data e...
Definition Vector.h:372
dolfinx::la::Vector::scatter_fwd_end
void scatter_fwd_end(U unpack)
End scatter (send) of local data values that are ghosted on other processes.
Definition Vector.h:256
dolfinx::la::Vector::scatter_rev_end
void scatter_rev_end(U unpack)
End scatter of ghost data to owner and update owned entries.
Definition Vector.h:349
dolfinx::la::Vector::operator=
Vector & operator=(Vector &&x)=default
Move assignment operator.
dolfinx::la::Vector::index_map
std::shared_ptr< const common::IndexMap > index_map() const
Get IndexMap.
Definition Vector.h:379
dolfinx::la::Vector::scatter_fwd_end
void scatter_fwd_end()
End scatter (send) of local data values that are ghosted on other processes (simplified CPU version).
Definition Vector.h:272
dolfinx::la::Vector::scatter_fwd
void scatter_fwd()
Scatter (send) of local data values that are ghosted on other processes and update ghost entry values...
Definition Vector.h:289
dolfinx::la::Vector::container_type
Container container_type
Container type.
Definition Vector.h:118
dolfinx::la::Vector::scatter_fwd_begin
void scatter_fwd_begin(U pack, GetPtr get_ptr)
Begin scatter (send) of local data that is ghosted on other processes.
Definition Vector.h:219
dolfinx::la::Vector::set
void set(value_type v)
Set all entries (including ghosts).
Definition Vector.h:197
dolfinx::la::Vector::Vector
Vector(Vector &&x)=default
Move constructor.
dolfinx::la::Vector::scatter_fwd_begin
void scatter_fwd_begin()
Begin scatter (send) of local data that is ghosted on other processes (simplified CPU version).
Definition Vector.h:235
dolfinx::la::Vector::Vector
Vector(const Vector &x)=default
Copy constructor.
dolfinx::la::Vector::scatter_rev_begin
void scatter_rev_begin()
Start scatter (send) of ghost entry data to the owning process of an index (simplified CPU version).
Definition Vector.h:332
dolfinx::la::GetPtrConcept
Access to pointer function concept.
Definition Vector.h:33
dolfinx::la::VectorPackKernel
la::Vector scatter pack/unpack function concept.
Definition Vector.h:27
dolfinx::MPI::mpi_t
MPI_Datatype mpi_t
Retrieves the MPI data type associated to the provided type.
Definition MPI.h:280
dolfinx::common
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
dolfinx::la
Linear algebra interface.
Definition sparsitybuild.h:15
dolfinx::la::orthonormalize
void orthonormalize(std::vector< std::reference_wrapper< V > > basis)
Orthonormalize a set of vectors.
Definition Vector.h:520
dolfinx::la::squared_norm
auto squared_norm(const V &a)
Compute the squared L2 norm of vector.
Definition Vector.h:463
dolfinx::la::norm
auto norm(const V &x, Norm type=Norm::l2)
Compute the norm of the vector.
Definition Vector.h:477
dolfinx::la::inner_product
auto inner_product(const V &a, const V &b)
Compute the inner product of two vectors.
Definition Vector.h:432
dolfinx::la::is_orthonormal
bool is_orthonormal(std::vector< std::reference_wrapper< const V > > basis, dolfinx::scalar_value_t< typename V::value_type > eps=std::numeric_limits< dolfinx::scalar_value_t< typename V::value_type > >::epsilon())
Test if basis is orthonormal.
Definition Vector.h:563
dolfinx::la::Norm
Norm
Norm types.
Definition utils.h:17