6#include <dune/functions/functionspacebases/flatvectorview.hh>
8#include <amdis/common/Math.hpp>
9#include <amdis/operations/Arithmetic.hpp>
10#include <amdis/operations/MaxMin.hpp>
19 auto negate(A
const& a)
21 return multiplies(a, -1);
25 template <
class A,
class B>
26 auto plus(A
const& a, B
const& b)
28 if constexpr(IsNumber<A>::value && IsNumber<B>::value)
31 using T = std::common_type_t<typename FieldTraits<A>::field_type,
typename FieldTraits<B>::field_type>;
32 typename MakeMatVec<A,T>::type c(a);
34 auto b_ = Dune::Functions::flatVectorView(b);
35 auto c_ = Dune::Functions::flatVectorView(c);
36 assert(
int(b_.size()) ==
int(c_.size()));
37 for(
int i = 0; i <
int(b_.size()); ++i)
45 template <
class A,
class B>
46 auto minus(A
const& a, B
const& b)
48 if constexpr(IsNumber<A>::value && IsNumber<B>::value)
51 using T = std::common_type_t<typename FieldTraits<A>::field_type,
typename FieldTraits<B>::field_type>;
52 typename MakeMatVec<A,T>::type c(a);
54 auto b_ = Dune::Functions::flatVectorView(b);
55 auto c_ = Dune::Functions::flatVectorView(c);
56 assert(
int(b_.size()) ==
int(c_.size()));
57 for(
int i = 0; i <
int(b_.size()); ++i)
64 template <
class A,
class B,
65 std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value,
int>>
66 auto multiplies(A
const& a, B
const& b)
71 template <
class A,
class B,
72 std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value,
int>>
73 auto multiplies(A
const& a, B
const& b)
75 using T = std::common_type_t<typename FieldTraits<A>::field_type,
typename FieldTraits<B>::field_type>;
76 if constexpr(IsNumber<A>::value) {
77 typename MakeMatVec<B,T>::type b_(b);
80 typename MakeMatVec<A,T>::type a_(a);
85 template <
class T,
int N,
class S>
86 auto multiplies(FieldVector<T,N>
const& a, FieldVector<S,N>
const& b)
92 template <
class Mat,
class Vec,
93 std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value,
int>>
94 auto multiplies(Mat
const& mat, Vec
const& vec)
96 static_assert(
int(Mat::cols) ==
int(Vec::dimension),
"");
97 using T = std::common_type_t<typename FieldTraits<Vec>::field_type,
typename FieldTraits<Mat>::field_type>;
98 FieldVector<T,Mat::rows> y;
104 template <
class Vec,
class Mat,
105 std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value,
int>>
106 auto multiplies(Vec
const& vec, Mat
const& mat)
108 static_assert(
int(Mat::rows) ==
int(Vec::dimension),
"");
109 using T = std::common_type_t<typename FieldTraits<Vec>::field_type,
typename FieldTraits<Mat>::field_type>;
110 FieldVector<T,Mat::cols> y;
116 template <
class T,
int L,
int M,
int N,
class S>
117 auto multiplies(FieldMatrix<T,L,M>
const& a, FieldMatrix<S,M,N>
const& b)
119 FieldMatrix<std::common_type_t<T,S>,L,N> C;
121 for (
int i = 0; i < L; ++i) {
122 for (
int j = 0; j < N; ++j) {
124 for (
int k = 0; k < M; ++k)
125 C[i][j] += a[i][k]*b[k][j];
131 template <
class T,
int SIZE>
132 class MatrixView<DiagonalMatrix<T,SIZE>>
134 using Matrix = DiagonalMatrix<T,SIZE>;
135 using size_type =
typename Matrix::size_type;
139 T operator[](size_type c)
const
141 assert(0 <= c && c < mat_->M());
142 return c == r_ ? mat_->diagonal(r_) : T(0);
150 MatrixView(Matrix
const& mat)
154 RowView operator[](size_type r)
const
156 assert(0 <= r && r < mat_.N());
180FieldVector<T, 2> cross(FieldVector<T, 2>
const& a)
182 return {{ a[1], -a[0] }};
187FieldVector<T, 3> cross(FieldVector<T, 3>
const& a, FieldVector<T, 3>
const& b)
189 return {{ a[1]*b[2] - a[2]*b[1],
190 a[2]*b[0] - a[0]*b[2],
191 a[0]*b[1] - a[1]*b[0] }};
195template <
class T,
class S,
int N>
196auto dot(FieldVector<T,N>
const& vec1, FieldVector<S,N>
const& vec2)
198 return vec1.dot(vec2);
201template <
class T,
class S,
int N>
202auto dot(FieldMatrix<T,1,N>
const& vec1, FieldMatrix<S,1,N>
const& vec2)
204 return vec1[0].dot(vec2[0]);
211 template <
class T,
int N,
class Operation>
212 T accumulate(FieldVector<T, N>
const& x, T init, Operation op)
214 for (
int i = 0; i < N; ++i)
215 init = op(init, x[i]);
219 template <
class T,
int N,
class Operation>
220 T accumulate(FieldMatrix<T, 1, N>
const& x, T init, Operation op)
222 for (
int i = 0; i < N; ++i)
223 init = op(init, x[0][i]);
230template <
class T,
int N>
231T sum(FieldVector<T, N>
const& x)
236template <
class T,
int N>
237T sum(FieldMatrix<T, 1, N>
const& x)
245 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
246auto unary_dot(T
const& x)
249 return AMDiS::Math::sqr(abs(x));
252template <
class T,
int N>
253auto unary_dot(FieldVector<T, N>
const& x)
255 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + AMDiS::Math::sqr(abs(b)); };
256 return Impl::accumulate(x, T(0), op);
259template <
class T,
int N>
260auto unary_dot(FieldMatrix<T, 1, N>
const& x)
262 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + AMDiS::Math::sqr(abs(b)); };
263 return Impl::accumulate(x, T(0), op);
267template <
class T,
int N>
268auto max(FieldVector<T, N>
const& x)
273template <
class T,
int N>
274auto max(FieldMatrix<T, 1, N>
const& x)
280template <
class T,
int N>
281auto min(FieldVector<T, N>
const& x)
286template <
class T,
int N>
287auto min(FieldMatrix<T, 1, N>
const& x)
293template <
class T,
int N>
294auto abs_max(FieldVector<T, N>
const& x)
299template <
class T,
int N>
300auto abs_max(FieldMatrix<T, 1, N>
const& x)
306template <
class T,
int N>
307auto abs_min(FieldVector<T, N>
const& x)
312template <
class T,
int N>
313auto abs_min(FieldMatrix<T, 1, N>
const& x)
323template <
class T,
int N>
324auto one_norm(FieldVector<T, N>
const& x)
326 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + abs(b); };
327 return Impl::accumulate(x, T(0), op);
330template <
class T,
int N>
331auto one_norm(FieldMatrix<T, 1, N>
const& x)
333 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + abs(b); };
334 return Impl::accumulate(x, T(0), op);
341 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
342auto two_norm(T
const& x)
348template <
class T,
int N>
349auto two_norm(FieldVector<T, N>
const& x)
352 return sqrt(unary_dot(x));
355template <
class T,
int N>
356auto two_norm(FieldMatrix<T, 1, N>
const& x)
359 return sqrt(unary_dot(x));
365template <
int p,
class T,
int N>
366auto p_norm(FieldVector<T, N>
const& x)
368 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + AMDiS::Math::pow<p>(abs(b)); };
369 return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
372template <
int p,
class T,
int N>
373auto p_norm(FieldMatrix<T, 1, N>
const& x)
375 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + AMDiS::Math::pow<p>(abs(b)); };
376 return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
382template <
class T,
int N>
383auto infty_norm(FieldVector<T, N>
const& x)
388template <
class T,
int N>
389auto infty_norm(FieldMatrix<T, 1, N>
const& x)
398 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
399T distance(T
const& lhs, T
const& rhs)
402 return abs(lhs - rhs);
405template <
class T,
int N>
406T distance(FieldVector<T, N>
const& lhs, FieldVector<T, N>
const& rhs)
410 for (
int i = 0; i < N; ++i)
411 result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
418template <
class T,
class S,
int N,
int M,
int K>
419auto outer(FieldMatrix<T,N,K>
const& vec1, FieldMatrix<S,M,K>
const& vec2)
421 using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
423 for (
int i = 0; i < N; ++i)
424 for (
int j = 0; j < M; ++j)
425 mat[i][j] = vec1[i].dot(vec2[j]);
430template <
class T,
class S,
int N,
int M>
431auto outer(FieldVector<T,N>
const& vec1, FieldVector<S,M>
const& vec2)
433 using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
435 for (
int i = 0; i < N; ++i)
436 for (
int j = 0; j < M; ++j)
437 mat[i][j] = vec1[i] * vec2[j];
444T det(FieldMatrix<T, 0, 0>
const& )
451T det(FieldMatrix<T, 1, 1>
const& mat)
458T det(FieldMatrix<T, 2, 2>
const& mat)
460 return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
465T det(FieldMatrix<T, 3, 3>
const& mat)
467 return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1]
468 - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]);
472template <
class T,
int N>
473T det(FieldMatrix<T, N, N>
const& mat)
475 return mat.determinant();
479template <
class T,
int N>
480auto inv(FieldMatrix<T, N, N> mat)
487template <
class T,
int N>
488void solve(FieldMatrix<T, N, N>
const& A, FieldVector<T, N>& x, FieldVector<T, N>
const& b)
495template <
class T,
int N,
int M>
496T gramian(FieldMatrix<T,N,M>
const& DF)
499 return sqrt( det(outer(DF, DF)) );
503template <
class T,
int M>
504T gramian(FieldMatrix<T, 1, M>
const& DF)
507 return sqrt(dot(DF[0], DF[0]));
513template <
class T,
int M,
int N>
514FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N>
const& A)
516 FieldMatrix<T,N,M> At;
517 for (
int i = 0; i < M; ++i)
518 for (
int j = 0; j < N; ++j)
525template <
class T1,
class T2,
int M,
int N,
int L>
526FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M>
const& A, FieldMatrix<T2, N, L>
const& B)
528 FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
530 for (
int m = 0; m < M; ++m) {
531 for (
int n = 0; n < N; ++n) {
533 for (
int l = 0; l < L; ++l)
534 C[m][n] += A[l][m] * B[n][l];
540template <
class T1,
class T2,
int M,
int N,
int L>
541FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L>
const& A, FieldMatrix<T2, N, L>
const& B)
543 FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
544 return multiplies_ABt(A,B,C);
547template <
class T1,
class T2,
class T3,
int M,
int N,
int L>
548FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L>
const& A, FieldMatrix<T2, N, L>
const& B, FieldMatrix<T3,M,N>& C)
550 for (
int m = 0; m < M; ++m) {
551 for (
int n = 0; n < N; ++n) {
553 for (
int l = 0; l < L; ++l)
554 C[m][n] += A[m][l] * B[n][l];
560template <
class T1,
class T2,
class T3,
int N,
int L>
561FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L>
const& A, FieldMatrix<T2, N, L>
const& B, FieldVector<T3,N>& C)
563 for (
int n = 0; n < N; ++n) {
565 for (
int l = 0; l < L; ++l)
566 C[n] += A[0][l] * B[n][l];
571template <
class T1,
class T2,
class T3,
int N,
int L>
572FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L>
const& A, FieldMatrix<T2, N, L>
const& B, FieldVector<T3,N>& C)
574 for (
int n = 0; n < N; ++n) {
576 for (
int l = 0; l < L; ++l)
577 C[n] += A[l] * B[n][l];
582template <
class T1,
class T2,
class T3,
int N,
int L>
583FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L>
const& A, FieldMatrix<T2, N, L>
const& B, FieldMatrix<T3,1,N>& C)
585 for (
int n = 0; n < N; ++n) {
587 for (
int l = 0; l < L; ++l)
588 C[0][n] += A[l] * B[n][l];
595template <
class T1,
class T2,
class T3,
int M,
int N>
596FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N>
const& A, DiagonalMatrix<T2, N>
const& B, FieldMatrix<T3,M,N>& C)
598 for (
int m = 0; m < M; ++m) {
599 for (
int n = 0; n < N; ++n) {
600 C[m][n] = A[m][n] * B.diagonal(n);
606template <
class T1,
class T2,
class T3,
int N>
607FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N>
const& A, DiagonalMatrix<T2, N>
const& B, FieldVector<T3,N>& C)
609 for (
int n = 0; n < N; ++n) {
610 C[n] = A[0][n] * B.diagonal(n);
615template <
class T1,
class T2,
class T3,
int N>
616FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N>
const& A, DiagonalMatrix<T2, N>
const& B, FieldVector<T3,N>& C)
618 for (
int n = 0; n < N; ++n) {
619 C[n] = A[n] * B.diagonal(n);
624template <
class T1,
class T2,
class T3,
int N>
625FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N>
const& A, DiagonalMatrix<T2, N>
const& B, FieldMatrix<T3,1,N>& C)
627 for (
int n = 0; n < N; ++n) {
628 C[0][n] = A[n] * B.diagonal(n);
634template <
class T,
int N>
635T
const& at(FieldMatrix<T,N,1>
const& vec, std::size_t i)
640template <
class T,
int M>
641T
const& at(FieldMatrix<T,1,M>
const& vec, std::size_t i)
647T
const& at(FieldMatrix<T,1,1>
const& vec, std::size_t i)
652template <
class T,
int N>
653T
const& at(FieldVector<T,N>
const& vec, std::size_t i)
Operation that represents max(|A|,|B|)
Definition: MaxMin.hpp:33
Operation that represents min(|A|,|B|)
Definition: MaxMin.hpp:46
Operation that represents max(A,B)
Definition: MaxMin.hpp:9
Operation that represents min(A,B)
Definition: MaxMin.hpp:21
Functor that represents A+B.
Definition: Arithmetic.hpp:20