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,
56template <std::
floating_po
int T>
57void dot_blas(std::span<const T>
A, std::array<std::size_t, 2>
Ashape,
58 std::span<const T>
B, std::array<std::size_t, 2>
Bshape,
61 static_assert(std::is_same_v<T, float>
or std::is_same_v<T, double>);
76 if constexpr (std::is_same_v<T, float>)
81 else if constexpr (std::is_same_v<T, double>)
94template <
typename U,
typename V>
95std::pair<std::vector<typename U::value_type>, std::array<std::size_t, 2>>
98 std::vector<typename U::value_type>
result(
u.size() *
v.size());
99 for (std::size_t
i = 0;
i <
u.size(); ++
i)
100 for (std::size_t
j = 0;
j <
v.size(); ++
j)
103 return {std::move(
result), {
u.size(),
v.size()}};
110template <
typename U,
typename V>
111std::array<typename U::value_type, 3>
cross(
const U&
u,
const V&
v)
115 return {
u[1] *
v[2] -
u[2] *
v[1],
u[2] *
v[0] -
u[0] *
v[2],
116 u[0] *
v[1] -
u[1] *
v[0]};
125template <std::
floating_po
int T>
126std::pair<std::vector<T>, std::vector<T>>
eigh(std::span<const T>
A,
130 std::vector<T> M(
A.begin(),
A.end());
133 std::vector<T>
w(
n, 0);
142 std::vector<T>
work(1);
143 std::vector<int>
iwork(1);
146 if constexpr (std::is_same_v<T, float>)
151 else if constexpr (std::is_same_v<T, double>)
158 throw std::runtime_error(
"Could not find workspace size for syevd.");
165 if constexpr (std::is_same_v<T, float>)
170 else if constexpr (std::is_same_v<T, double>)
176 throw std::runtime_error(
"Eigenvalue computation did not converge.");
178 return {std::move(
w), std::move(M)};
185template <std::
floating_po
int T>
187solve(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
188 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
190 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
191 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
195 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
198 stdex::mdarray<T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>,
199 MDSPAN_IMPL_STANDARD_NAMESPACE::layout_left>
200 _A(
A.extents()),
_B(
B.extents());
201 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
202 for (std::size_t
j = 0;
j <
A.extent(1); ++
j)
204 for (std::size_t
i = 0;
i <
B.extent(0); ++
i)
205 for (std::size_t
j = 0;
j <
B.extent(1); ++
j)
208 int N =
_A.extent(0);
210 int lda =
_A.extent(0);
211 int ldb =
_B.extent(0);
213 std::vector<int>
piv(
N);
215 if constexpr (std::is_same_v<T, float>)
217 else if constexpr (std::is_same_v<T, double>)
220 throw std::runtime_error(
"Call to dgesv failed: " + std::to_string(
info));
223 std::vector<T>
rb(
_B.extent(0) *
_B.extent(1));
224 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
225 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
226 r(
rb.data(),
_B.extents());
227 for (std::size_t
i = 0;
i <
_B.extent(0); ++
i)
228 for (std::size_t
j = 0;
j <
_B.extent(1); ++
j)
237template <std::
floating_po
int T>
239 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
240 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
245 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
246 stdex::mdarray<T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>,
247 MDSPAN_IMPL_STANDARD_NAMESPACE::layout_left>
249 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
250 for (std::size_t
j = 0;
j <
A.extent(1); ++
j)
253 std::vector<T>
B(
A.extent(1), 1);
254 int N =
_A.extent(0);
256 int lda =
_A.extent(0);
260 std::vector<int>
piv(
N);
262 if constexpr (std::is_same_v<T, float>)
264 else if constexpr (std::is_same_v<T, double>)
269 throw std::runtime_error(
"dgesv failed due to invalid value: "
270 + std::to_string(
info));
282template <std::
floating_po
int T>
283std::vector<std::size_t>
286 std::size_t dim =
A.second[0];
293 if constexpr (std::is_same_v<T, float>)
295 else if constexpr (std::is_same_v<T, double>)
300 throw std::runtime_error(
"LU decomposition failed: "
301 + std::to_string(
info));
304 std::vector<std::size_t>
perm(dim);
305 for (std::size_t
i = 0;
i < dim; ++
i)
316template <
typename U,
typename V,
typename W>
322 if (
A.extent(0) *
B.extent(1) *
A.extent(1) < 512)
324 std::fill_n(
C.data_handle(),
C.extent(0) *
C.extent(1), 0);
325 for (std::size_t
i = 0;
i <
A.extent(0); ++
i)
326 for (std::size_t
j = 0;
j <
B.extent(1); ++
j)
327 for (std::size_t
k = 0;
k <
A.extent(1); ++
k)
332 using T =
typename std::decay_t<U>::value_type;
334 std::span(
A.data_handle(),
A.size()), {A.extent(0), A.extent(1)},
335 std::span(
B.data_handle(),
B.size()), {B.extent(0), B.extent(1)},
336 std::span(
C.data_handle(),
C.size()));
343template <std::
floating_po
int T>
344std::vector<T>
eye(std::size_t
n)
346 std::vector<T>
I(
n *
n, 0);
348 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
349 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
350 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
352 for (std::size_t
i = 0;
i <
n; ++
i)
361template <std::
floating_po
int T>
363 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
364 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
366 std::size_t
start = 0)
368 for (std::size_t
i =
start;
i < wcoeffs.extent(0); ++
i)
373 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
374 a += wcoeffs(
i,
k) * wcoeffs(
j,
k);
375 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
376 wcoeffs(
i,
k) -=
a * wcoeffs(
j,
k);
380 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
381 norm += wcoeffs(
i,
k) * wcoeffs(
i,
k);
382 if (std::abs(
norm) < 4 * std::numeric_limits<T>::epsilon())
384 throw std::runtime_error(
385 "Cannot orthogonalise the rows of a matrix with incomplete row rank");
388 for (std::size_t
k = 0;
k < wcoeffs.extent(1); ++
k)
389 wcoeffs(
i,
k) /= std::sqrt(
norm);
A finite element.
Definition finite-element.h:139
void dot(const U &A, const V &B, W &&C)
Definition math.h:317
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)
Definition math.h:187
std::array< typename U::value_type, 3 > cross(const U &u, const V &v)
Definition math.h:111
bool is_singular(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > A)
Definition math.h:238
std::vector< std::size_t > transpose_lu(std::pair< std::vector< T >, std::array< std::size_t, 2 > > &A)
Definition math.h:284
void orthogonalise(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > wcoeffs, std::size_t start=0)
Definition math.h:362
std::pair< std::vector< T >, std::vector< T > > eigh(std::span< const T > A, std::size_t n)
Definition math.h:126
std::vector< T > eye(std::size_t n)
Definition math.h:344
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:96