20 void ssyevd_(
char* jobz,
char* uplo,
int* n,
float* a,
int* lda,
float* w,
21 float* work,
int* lwork,
int* iwork,
int* liwork,
int* info);
22 void dsyevd_(
char* jobz,
char* uplo,
int* n,
double* a,
int* lda,
double* w,
23 double* work,
int* lwork,
int* iwork,
int* liwork,
int* info);
25 void sgesv_(
int* N,
int* NRHS,
float* A,
int* LDA,
int* IPIV,
float* B,
27 void dgesv_(
int* N,
int* NRHS,
double* A,
int* LDA,
int* IPIV,
double* B,
30 void sgemm_(
char* transa,
char* transb,
int* m,
int* n,
int* k,
float* alpha,
31 float* a,
int* lda,
float* b,
int* ldb,
float* beta,
float* c,
33 void dgemm_(
char* transa,
char* transb,
int* m,
int* n,
int* k,
double* alpha,
34 double* a,
int* lda,
double* b,
int* ldb,
double* beta,
double* c,
37 int sgetrf_(
const int* m,
const int* n,
float* a,
const int* lda,
int* lpiv,
39 int dgetrf_(
const int* m,
const int* n,
double* a,
const int* lda,
int* lpiv,
55template <std::
floating_po
int T>
56void dot_blas(std::span<const T>
A, std::array<std::size_t, 2>
Ashape,
57 std::span<const T>
B, std::array<std::size_t, 2>
Bshape,
60 static_assert(std::is_same_v<T, float>
or std::is_same_v<T, double>);
75 if constexpr (std::is_same_v<T, float>)
80 else if constexpr (std::is_same_v<T, double>)
93template <
typename U,
typename V>
94std::pair<std::vector<typename U::value_type>, std::array<std::size_t, 2>>
97 std::vector<typename U::value_type>
result(
u.size() *
v.size());
98 for (std::size_t
i = 0;
i <
u.size(); ++
i)
99 for (std::size_t
j = 0;
j <
v.size(); ++
j)
101 return {std::move(
result), {
u.size(),
v.size()}};
108template <
typename U,
typename V>
109std::array<typename U::value_type, 3>
cross(
const U&
u,
const V&
v)
113 return {
u[1] *
v[2] -
u[2] *
v[1],
u[2] *
v[0] -
u[0] *
v[2],
114 u[0] *
v[1] -
u[1] *
v[0]};
123template <std::
floating_po
int T>
124std::pair<std::vector<T>, std::vector<T>>
eigh(std::span<const T>
A,
128 std::vector<T> M(
A.begin(),
A.end());
131 std::vector<T>
w(
n, 0);
140 std::vector<T>
work(1);
141 std::vector<int>
iwork(1);
144 if constexpr (std::is_same_v<T, float>)
149 else if constexpr (std::is_same_v<T, double>)
156 throw std::runtime_error(
"Could not find workspace size for syevd.");
163 if constexpr (std::is_same_v<T, float>)
168 else if constexpr (std::is_same_v<T, double>)
174 throw std::runtime_error(
"Eigenvalue computation did not converge.");
176 return {std::move(
w), std::move(M)};
183template <std::
floating_po
int T>
185solve(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
186 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
188 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
189 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
193 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
196 stdex::mdarray<T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>,
197 MDSPAN_IMPL_STANDARD_NAMESPACE::layout_left>
198 _A(
A.extents()),
_B(
B.extents());
199 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
200 for (std::size_t
j = 0;
j <
A.extent(1); ++
j)
202 for (std::size_t
i = 0;
i <
B.extent(0); ++
i)
203 for (std::size_t
j = 0;
j <
B.extent(1); ++
j)
206 int N =
_A.extent(0);
208 int lda =
_A.extent(0);
209 int ldb =
_B.extent(0);
211 std::vector<int>
piv(
N);
213 if constexpr (std::is_same_v<T, float>)
215 else if constexpr (std::is_same_v<T, double>)
218 throw std::runtime_error(
"Call to dgesv failed: " + std::to_string(
info));
221 std::vector<T>
rb(
_B.extent(0) *
_B.extent(1));
222 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
223 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
224 r(
rb.data(),
_B.extents());
225 for (std::size_t
i = 0;
i <
_B.extent(0); ++
i)
226 for (std::size_t
j = 0;
j <
_B.extent(1); ++
j)
235template <std::
floating_po
int T>
237 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
238 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
243 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
244 stdex::mdarray<T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>,
245 MDSPAN_IMPL_STANDARD_NAMESPACE::layout_left>
247 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
248 for (std::size_t
j = 0;
j <
A.extent(1); ++
j)
251 std::vector<T>
B(
A.extent(1), 1);
252 int N =
_A.extent(0);
254 int lda =
_A.extent(0);
258 std::vector<int>
piv(
N);
260 if constexpr (std::is_same_v<T, float>)
262 else if constexpr (std::is_same_v<T, double>)
267 throw std::runtime_error(
"dgesv failed due to invalid value: "
268 + std::to_string(
info));
281template <std::
floating_po
int T>
282std::vector<std::size_t>
285 std::size_t dim =
A.second[0];
292 if constexpr (std::is_same_v<T, float>)
294 else if constexpr (std::is_same_v<T, double>)
299 throw std::runtime_error(
"LU decomposition failed: "
300 + std::to_string(
info));
303 std::vector<std::size_t>
perm(dim);
304 for (std::size_t
i = 0;
i < dim; ++
i)
315template <
typename U,
typename V,
typename W>
321 if (
A.extent(0) *
B.extent(1) *
A.extent(1) < 512)
323 std::fill_n(
C.data_handle(),
C.extent(0) *
C.extent(1), 0);
324 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
325 for (std::size_t
j = 0;
j <
B.extent(1); ++
j)
326 for (std::size_t
k = 0;
k <
A.extent(1); ++
k)
331 using T =
typename std::decay_t<U>::value_type;
333 std::span(
A.data_handle(),
A.size()), {A.extent(0), A.extent(1)},
334 std::span(
B.data_handle(),
B.size()), {B.extent(0), B.extent(1)},
335 std::span(
C.data_handle(),
C.size()));
342template <std::
floating_po
int T>
343std::vector<T>
eye(std::size_t
n)
345 std::vector<T>
I(
n *
n, 0);
347 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
348 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
349 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
351 for (std::size_t
i = 0;
i <
n; ++
i)
360template <std::
floating_po
int T>
362 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
363 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
365 std::size_t
start = 0)
367 for (std::size_t
i =
start;
i < wcoeffs.extent(0); ++
i)
370 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
371 norm += wcoeffs(
i,
k) * wcoeffs(
i,
k);
374 if (
norm < 2 * std::numeric_limits<T>::epsilon())
376 throw std::runtime_error(
377 "Cannot orthogonalise the rows of a matrix with incomplete row rank");
380 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
383 for (std::size_t
j =
i + 1;
j < wcoeffs.extent(0); ++
j)
386 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
387 a += wcoeffs(
i,
k) * wcoeffs(
j,
k);
388 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
389 wcoeffs(
j,
k) -=
a * wcoeffs(
i,
k);
A finite element.
Definition finite-element.h:139
Mathematical functions.
Definition math.h:48
void dot(const U &A, const V &B, W &&C)
Compute C = A * B.
Definition math.h:316
std::vector< T > solve(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > A, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > B)
Solve A X = B.
Definition math.h:185
std::array< typename U::value_type, 3 > cross(const U &u, const V &v)
Definition math.h:109
bool is_singular(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > A)
Check if A is a singular matrix,.
Definition math.h:236
std::vector< std::size_t > transpose_lu(std::pair< std::vector< T >, std::array< std::size_t, 2 > > &A)
Compute the LU decomposition of the transpose of a square matrix A.
Definition math.h:283
void orthogonalise(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > wcoeffs, std::size_t start=0)
Orthogonalise the rows of a matrix (in place).
Definition math.h:361
std::pair< std::vector< T >, std::vector< T > > eigh(std::span< const T > A, std::size_t n)
Definition math.h:124
std::vector< T > eye(std::size_t n)
Build an identity matrix.
Definition math.h:343
std::pair< std::vector< typename U::value_type >, std::array< std::size_t, 2 > > outer(const U &u, const V &v)
Compute the outer product of vectors u and v.
Definition math.h:95