Basix
precompute.h
1// Copyright (c) 2020 Matthew Scroggs
2// FEniCS Project
3// SPDX-License-Identifier: MIT
4
5#pragma once
6
7#include "mdspan.hpp"
8#include <span>
9#include <tuple>
10#include <vector>
11
14{
15
63void prepare_permutation(const std::span<std::size_t>& perm);
64
116template <typename E>
117void apply_permutation(const std::span<const std::size_t>& perm,
118 const std::span<E>& data, std::size_t offset = 0,
119 std::size_t block_size = 1)
120{
121 for (std::size_t b = 0; b < block_size; ++b)
122 {
123 for (std::size_t i = 0; i < perm.size(); ++i)
124 {
125 std::swap(data[block_size * (offset + i) + b],
126 data[block_size * (offset + perm[i]) + b]);
127 }
128 }
129}
130
137template <typename E>
138void apply_permutation_to_transpose(const std::span<const std::size_t>& perm,
139 const std::span<E>& data,
140 std::size_t offset = 0,
141 std::size_t block_size = 1)
142{
143 const std::size_t dim = perm.size();
144 const std::size_t data_size
145 = (data.size() + (dim < block_size ? block_size - dim : 0)) / block_size;
146 for (std::size_t b = 0; b < block_size; ++b)
147 {
148 for (std::size_t i = 0; i < dim; ++i)
149 {
150 std::swap(data[data_size * b + offset + i],
151 data[data_size * b + offset + perm[i]]);
152 }
153 }
154}
155
170std::vector<std::size_t>
171prepare_matrix(std::pair<std::vector<double>, std::array<std::size_t, 2>>& A);
172
206template <typename T, typename E>
207void apply_matrix(const std::span<const std::size_t>& v_size_t,
208 const std::experimental::mdspan<
209 const T, std::experimental::dextents<std::size_t, 2>>& M,
210 const std::span<E>& data, std::size_t offset = 0,
211 std::size_t block_size = 1)
212{
213 const std::size_t dim = v_size_t.size();
214 apply_permutation(v_size_t, data, offset, block_size);
215 for (std::size_t b = 0; b < block_size; ++b)
216 {
217 for (std::size_t i = 0; i < dim; ++i)
218 {
219 for (std::size_t j = i + 1; j < dim; ++j)
220 {
221 data[block_size * (offset + i) + b]
222 += M(i, j) * data[block_size * (offset + j) + b];
223 }
224 }
225 for (std::size_t i = 1; i <= dim; ++i)
226 {
227 data[block_size * (offset + dim - i) + b] *= M(dim - i, dim - i);
228 for (std::size_t j = 0; j < dim - i; ++j)
229 {
230 data[block_size * (offset + dim - i) + b]
231 += M(dim - i, j) * data[block_size * (offset + j) + b];
232 }
233 }
234 }
235}
236
243template <typename T, typename E>
245 const std::span<const std::size_t>& v_size_t,
246 const std::experimental::mdspan<
247 const T, std::experimental::dextents<std::size_t, 2>>& M,
248 const std::span<E>& data, std::size_t offset = 0,
249 std::size_t block_size = 1)
250{
251 const std::size_t dim = v_size_t.size();
252 const std::size_t data_size
253 = (data.size() + (dim < block_size ? block_size - dim : 0)) / block_size;
254 apply_permutation_to_transpose(v_size_t, data, offset, block_size);
255 for (std::size_t b = 0; b < block_size; ++b)
256 {
257 for (std::size_t i = 0; i < dim; ++i)
258 {
259 for (std::size_t j = i + 1; j < dim; ++j)
260 {
261 data[data_size * b + offset + i]
262 += M(i, j) * data[data_size * b + offset + j];
263 }
264 }
265 for (std::size_t i = 1; i <= dim; ++i)
266 {
267 data[data_size * b + offset + dim - i] *= M(dim - i, dim - i);
268 for (std::size_t j = 0; j < dim - i; ++j)
269 {
270 data[data_size * b + offset + dim - i]
271 += M(dim - i, j) * data[data_size * b + offset + j];
272 }
273 }
274 }
275}
276
277} // namespace basix::precompute
Matrix and permutation precomputation.
Definition: precompute.h:14
void apply_permutation(const std::span< const std::size_t > &perm, const std::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:117
void apply_matrix_to_transpose(const std::span< const std::size_t > &v_size_t, const std::experimental::mdspan< const T, std::experimental::dextents< std::size_t, 2 > > &M, const 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:244
void apply_matrix(const std::span< const std::size_t > &v_size_t, const std::experimental::mdspan< const T, std::experimental::dextents< std::size_t, 2 > > &M, const std::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Apply a (precomputed) matrix.
Definition: precompute.h:207
void apply_permutation_to_transpose(const std::span< const std::size_t > &perm, const std::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:138
std::vector< std::size_t > prepare_matrix(std::pair< std::vector< double >, std::array< std::size_t, 2 > > &A)
Definition: precompute.cpp:26
void prepare_permutation(const std::span< std::size_t > &perm)
Definition: precompute.cpp:17