Basix
Loading...
Searching...
No Matches
precompute.h
1// Copyright (c) 2020 Matthew Scroggs
2// FEniCS Project
3// SPDX-License-Identifier: MIT
4
5#pragma once
6
7#include "math.h"
8#include "mdspan.hpp"
9#include <concepts>
10#include <span>
11#include <tuple>
12#include <type_traits>
13#include <vector>
14
17{
18
19namespace impl
20{
23template <typename T, typename = void>
24struct scalar_value_type
25{
27 typedef T value_type;
28};
30template <typename T>
31struct scalar_value_type<T, std::void_t<typename T::value_type>>
32{
33 typedef typename T::value_type value_type;
34};
36template <typename T>
37using scalar_value_type_t = typename scalar_value_type<T>::value_type;
38} // namespace impl
39
87void prepare_permutation(std::span<std::size_t> perm);
88
140template <typename E>
141void apply_permutation(std::span<const std::size_t> perm, std::span<E> data,
142 std::size_t offset = 0, std::size_t block_size = 1)
143{
144 for (std::size_t i = 0; i < perm.size(); ++i)
145 {
146 for (std::size_t b = 0; b < block_size; ++b)
147 {
148 std::swap(data[block_size * (offset + i) + b],
149 data[block_size * (offset + perm[i]) + b]);
150 }
151 }
152}
153
155template <typename E>
156void apply_permutation_mapped(std::span<const std::size_t> perm,
157 std::span<E> data, std::span<const int> emap,
158 std::size_t block_size = 1)
159{
160 for (std::size_t i = 0; i < perm.size(); ++i)
161 {
162 for (std::size_t b = 0; b < block_size; ++b)
163 {
164 std::swap(data[block_size * emap[i] + b],
165 data[block_size * emap[perm[i]] + b]);
166 }
167 }
168}
169
176template <typename E>
177void apply_permutation_to_transpose(std::span<const std::size_t> perm,
178 std::span<E> data, std::size_t offset = 0,
179 std::size_t block_size = 1)
180{
181 const std::size_t dim = perm.size();
182 const std::size_t data_size
183 = (data.size() + (dim < block_size ? block_size - dim : 0)) / block_size;
184 for (std::size_t b = 0; b < block_size; ++b)
185 {
186 for (std::size_t i = 0; i < dim; ++i)
187 {
188 std::swap(data[data_size * b + offset + i],
189 data[data_size * b + offset + perm[i]]);
190 }
191 }
192}
193
210template <std::floating_point T>
211std::vector<std::size_t>
212prepare_matrix(std::pair<std::vector<T>, std::array<std::size_t, 2>>& A)
213{
214 return math::transpose_lu<T>(A);
215}
216
250template <typename T, typename E>
252 std::span<const std::size_t> v_size_t,
253 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
254 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
255 M,
256 std::span<E> data, std::size_t offset = 0, std::size_t block_size = 1)
257{
258 using U = typename impl::scalar_value_type_t<E>;
259
260 const std::size_t dim = v_size_t.size();
262 for (std::size_t b = 0; b < block_size; ++b)
263 {
264 for (std::size_t i = 0; i < dim; ++i)
265 {
266 for (std::size_t j = i + 1; j < dim; ++j)
267 {
268 data[block_size * (offset + i) + b]
269 += static_cast<U>(M(i, j)) * data[block_size * (offset + j) + b];
270 }
271 }
272 for (std::size_t i = 1; i <= dim; ++i)
273 {
274 data[block_size * (offset + dim - i) + b]
275 *= static_cast<U>(M(dim - i, dim - i));
276 for (std::size_t j = 0; j < dim - i; ++j)
277 {
278 data[block_size * (offset + dim - i) + b]
279 += static_cast<U>(M(dim - i, j))
280 * data[block_size * (offset + j) + b];
281 }
282 }
283 }
284}
285
292template <typename T, typename E>
294 std::span<const std::size_t> v_size_t,
295 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
296 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
297 M,
298 std::span<E> data, std::size_t offset = 0, std::size_t block_size = 1)
299{
300 using U = typename impl::scalar_value_type_t<E>;
301
302 const std::size_t dim = v_size_t.size();
303 const std::size_t data_size
304 = (data.size() + (dim < block_size ? block_size - dim : 0)) / block_size;
306 for (std::size_t b = 0; b < block_size; ++b)
307 {
308 for (std::size_t i = 0; i < dim; ++i)
309 {
310 for (std::size_t j = i + 1; j < dim; ++j)
311 {
312 data[data_size * b + offset + i]
313 += static_cast<U>(M(i, j)) * data[data_size * b + offset + j];
314 }
315 }
316 for (std::size_t i = 1; i <= dim; ++i)
317 {
318 data[data_size * b + offset + dim - i]
319 *= static_cast<U>(M(dim - i, dim - i));
320 for (std::size_t j = 0; j < dim - i; ++j)
321 {
322 data[data_size * b + offset + dim - i]
323 += static_cast<U>(M(dim - i, j)) * data[data_size * b + offset + j];
324 }
325 }
326 }
327}
328
329} // namespace basix::precompute
A finite element.
Definition finite-element.h:139
Matrix and permutation precomputation.
Definition precompute.h:17
void apply_permutation(std::span< const std::size_t > perm, std::span< E > data, std::size_t offset=0, std::size_t block_size=1)
Definition precompute.h:141
void apply_permutation_to_transpose(std::span< const std::size_t > perm, std::span< E > data, std::size_t offset=0, std::size_t block_size=1)
Definition precompute.h:177
void apply_matrix_to_transpose(std::span< const std::size_t > v_size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > M, std::span< E > data, std::size_t offset=0, std::size_t block_size=1)
Apply a (precomputed) matrix to some transposed data.
Definition precompute.h:293
std::vector< std::size_t > prepare_matrix(std::pair< std::vector< T >, std::array< std::size_t, 2 > > &A)
Definition precompute.h:212
void prepare_permutation(std::span< std::size_t > perm)
Definition precompute.cpp:10
void apply_permutation_mapped(std::span< const std::size_t > perm, std::span< E > data, std::span< const int > emap, std::size_t block_size=1)
Permutation of mapped data.
Definition precompute.h:156
void apply_matrix(std::span< const std::size_t > v_size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > M, std::span< E > data, std::size_t offset=0, std::size_t block_size=1)
Apply a (precomputed) matrix.
Definition precompute.h:251