Skip to content

GitLab

Commit 2494d6df authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/fixmatvec_operations' into 'master' 

densematrix view for sparse matrices

See merge request amdis/amdis-core!15
parents 557efb43 9e0fc1aa
Hide whitespace changes
Inline Side-by-side
src/amdis/common/FieldMatVec.hpp
View file @ 2494d6df
......@@ -82,51 +82,88 @@ namespace Dune
/// Convert pseudo-scalar to real scalar type
/// @{
template <class T>
decltype(auto) simplify(T&& t)
decltype(auto) as_scalar(T&& t)
{
return FWD(t);
}
template <class T>
T simplify(FieldVector<T,1> const& t)
T as_scalar(FieldVector<T,1> const& t)
{
return t[0];
}
template <class T>
T simplify(FieldMatrix<T,1,1> const& t)
T as_scalar(FieldMatrix<T,1,1> const& t)
{
return t[0][0];
}
template <class T>
T simplify(DiagonalMatrix<T,1> const& t)
T as_scalar(DiagonalMatrix<T,1> const& t)
{
return t.diagonal(0);
}
/// @}
/// Convert pseudo-vector to real vector type
/// @{
template <class T>
decltype(auto) as_vector(T&& t)
{
return FWD(t);
}
template <class T, int N>
FieldVector<T,N> const& as_vector(FieldMatrix<T,1,N> const& t)
{
return t[0];
}
template <class T, int N>
FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
{
return t[0];
}
/// @}
/// Convert pseudo-matrix to real matrix type with proper operator[][]
/// @{
template <class T>
decltype(auto) as_matrix(T&& t)
{
return FWD(t);
}
template <class Mat>
class MatrixView;
template <class T, int N>
MatrixView<DiagonalMatrix<T,N>> as_matrix(DiagonalMatrix<T,N> const& mat)
{
return {mat};
}
// returns -a
template <class A>
auto negate(A const& a);
// returns a+b
template <class A, class B>
auto plus(A a, B const& b)
{
return a += b;
}
auto plus(A const& a, B const& b);
// returns a-b
template <class A, class B>
auto minus(A a, B const& b)
{
return a -= b;
}
auto minus(A const& a, B const& b);
// returns a*b
template <class A, class B,
std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int> = 0>
std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int> = 0>
auto multiplies(A const& a, B const& b);
template <class A, class B,
std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int> = 0>
auto multiplies(A const& a, B const& b);
template <class T, int N, class S>
......@@ -158,35 +195,35 @@ namespace Dune
std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
auto operator-(A const& a)
{
return MatVec::negate(MatVec::simplify(a));
return MatVec::negate(MatVec::as_scalar(a));
}
template <class A, class B,
std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
auto operator+(A const& a, B const& b)
{
return MatVec::plus(MatVec::simplify(a), MatVec::simplify(b));
return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
}
template <class A, class B,
std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
auto operator-(A const& a, B const& b)
{
return MatVec::minus(MatVec::simplify(a), MatVec::simplify(b));
return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
}
template <class A, class B,
std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
auto operator*(A const& a, B const& b)
{
return MatVec::multiplies(MatVec::simplify(a), MatVec::simplify(b));
return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
}
template <class A, class B,
std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
auto operator/(A const& a, B const& b)
{
return MatVec::divides(MatVec::simplify(a), MatVec::simplify(b));
return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
}
// ----------------------------------------------------------------------------
......@@ -346,19 +383,28 @@ namespace Dune
template <class T, int M, int N>
FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
template <class T, int N>
DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
{
return A;
}
// -----------------------------------------------------------------------------
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B);
template <class T1, class T2, int M, int N, int L>
FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A, FieldMatrix<T2, N, L> const& B);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B);
template <class T1, class T2, int M, int N, int L>
FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C);
template <class T1, class T2, class T3, int M, int N, int L>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
template <class T, int M, int N>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C);
template <class T1, class T2, class T3, int M, int N>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);
template <class T1, class T2, class T3, int N>
FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
// -----------------------------------------------------------------------------
......
src/amdis/common/FieldMatVec.inc.hpp
View file @ 2494d6df
......@@ -3,6 +3,8 @@
#include <algorithm>
#include <limits>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <amdis/common/Math.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>
......@@ -19,8 +21,47 @@ namespace MatVec {
return multiplies(a, -1);
}
// returns a+b
template <class A, class B>
auto plus(A const& a, B const& b)
{
using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
typename MakeMatVec<A,T>::type c{a};
auto b_ = Dune::Functions::flatVectorView(b);
auto c_ = Dune::Functions::flatVectorView(c);
assert(int(b_.size()) == int(c_.size()));
for(int i = 0; i < int(b_.size()); ++i)
c_[i] += b_[i];
return c;
}
// returns a-b
template <class A, class B>
auto minus(A const& a, B const& b)
{
using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
typename MakeMatVec<A,T>::type c{a};
auto b_ = Dune::Functions::flatVectorView(b);
auto c_ = Dune::Functions::flatVectorView(c);
assert(int(b_.size()) == int(c_.size()));
for(int i = 0; i < int(b_.size()); ++i)
c_[i] -= b_[i];
return c;
}
template <class A, class B,
std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int>>
std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int>>
auto multiplies(A const& a, B const& b)
{
return a * b;
}
template <class A, class B,
std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int>>
auto multiplies(A const& a, B const& b)
{
using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
......@@ -85,6 +126,49 @@ namespace MatVec {
return C;
}
template <class T, int SIZE>
class MatrixView<DiagonalMatrix<T,SIZE>>
{
using Matrix = DiagonalMatrix<T,SIZE>;
using size_type = typename Matrix::size_type;
struct RowView
{
T operator[](size_type c) const
{
assert(0 <= c && c < mat_->M());
return c == r_ ? mat_->diagonal(r_) : T(0);
}
Matrix const* mat_;
size_type r_;
};
public:
MatrixView(Matrix const& mat)
: mat_(mat)
{}
RowView operator[](size_type r) const
{
assert(0 <= r && r < mat_.N());
return {&mat_,r};
}
size_type N() const
{
return mat_.N();
}
size_type M() const
{
return mat_.M();
}
private:
Matrix const& mat_;
};
} // end namespace MatVec
// ----------------------------------------------------------------------------
......@@ -158,14 +242,14 @@ T sum(FieldMatrix<T, 1, N> const& x)
template <class T, int N>
auto unary_dot(FieldVector<T, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::sqr(std::abs(b)); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
return Impl::accumulate(x, T(0), op);
}
template <class T, int N>
auto unary_dot(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::sqr(std::abs(b)); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
return Impl::accumulate(x, T(0), op);
}
......@@ -229,14 +313,14 @@ auto abs_min(FieldMatrix<T, 1, N> const& x)
template <class T, int N>
auto one_norm(FieldVector<T, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
return Impl::accumulate(x, T(0), op);
}
template <class T, int N>
auto one_norm(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
return Impl::accumulate(x, T(0), op);
}
......@@ -261,14 +345,14 @@ auto two_norm(FieldMatrix<T, 1, N> const& x)
template <int p, class T, int N>
auto p_norm(FieldVector<T, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::pow<p>(std::abs(b)); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}
template <int p, class T, int N>
auto p_norm(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::pow<p>(std::abs(b)); };
auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}
......@@ -293,10 +377,11 @@ auto infty_norm(FieldMatrix<T, 1, N> const& x)
template <class T, int N>
T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
{
using std::sqrt;
T result = 0;
for (int i = 0; i < N; ++i)
result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
return std::sqrt(result);
return sqrt(result);
}
// ----------------------------------------------------------------------------
......@@ -397,10 +482,10 @@ FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
}
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B)
template <class T1, class T2, int M, int N, int L>
FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A, FieldMatrix<T2, N, L> const& B)
{
FieldMatrix<T,M,N> C;
FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
for (int m = 0; m < M; ++m) {
for (int n = 0; n < N; ++n) {
......@@ -412,15 +497,15 @@ FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T,
return C;
}
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B)
template <class T1, class T2, int M, int N, int L>
FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B)
{
FieldMatrix<T,M,N> C;
FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
return multiplies_ABt(A,B,C);
}
template <class T, int M, int N, int L>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C)
template <class T1, class T2, class T3, int M, int N, int L>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C)
{
for (int m = 0; m < M; ++m) {
for (int n = 0; n < N; ++n) {
......@@ -432,8 +517,8 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T
return C;
}
template <class T, int M, int N>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C)
template <class T1, class T2, class T3, int M, int N>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C)
{
for (int m = 0; m < M; ++m) {
for (int n = 0; n < N; ++n) {
......@@ -443,6 +528,15 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatri
return C;
}
template <class T1, class T2, class T3, int N>
FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
{
for (int n = 0; n < N; ++n) {
C[n] = A[0][n] * B.diagonal(n);
}
return C;
}
template <class T, int N>
T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
{
......
src/amdis/gridfunctions/DiscreteFunction.inc.hpp
View file @ 2494d6df
......@@ -3,6 +3,7 @@
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/utility/LocalBasisCache.hpp>
#include <dune/common/ftraits.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
namespace AMDiS {
......@@ -197,7 +198,7 @@ GradientLocalFunction::operator()(Domain const& x) const
{
// TODO: may DOFVectorView::Range to FieldVector type if necessary
using LocalDerivativeTraits
= Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
= Dune::Functions::DefaultDerivativeTraits<VT(Domain)>;
using GradientBlock = typename LocalDerivativeTraits::Range;
auto&& fe = node.finiteElement();
......@@ -211,8 +212,8 @@ GradientLocalFunction::operator()(Domain const& x) const
// Compute the shape function gradients on the real element
std::vector<GradientBlock> gradients(referenceGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
for (std::size_t i = 0; i < gradients.size(); ++i) // J^(-T) * D[phi] -> grad^T
Dune::MatVec::as_vector(gradients[i]) = jacobian * Dune::MatVec::as_vector(referenceGradients[i]);
// Get range entry associated to this node
auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
......
src/amdis/linearalgebra/mtl/Constraints.hpp
View file @ 2494d6df
......@@ -130,8 +130,7 @@ namespace AMDiS
ins[t.row][t.col] += t.value;
if (setDiagonal) {
using mtl::size;
for (std::size_t i = 0; i < size(left); ++i) {
for (std::size_t i = 0; i < mtl::size(left); ++i) {
if (left[i]) {
ins[i][i] = T(1);
ins[i][at(left2right,i)] = T(-1);
......
src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
View file @ 2494d6df
......@@ -69,6 +69,7 @@ namespace AMDiS
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
assert(jacobian.M() == CG::dim);
// The multiplicative factor in the integral transformation formula
......@@ -85,9 +86,9 @@ namespace AMDiS
thread_local std::vector<RangeFieldType> colPartial;
colPartial.resize(shapeGradients.size());
for (std::size_t i = 0; i < colPartial.size(); ++i) {
colPartial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.M(); ++j)
colPartial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
colPartial[i] = jacobian_mat[comp_][0] * shapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
colPartial[i] += jacobian_mat[comp_][j] * shapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colSize; ++j) {
......
src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
View file @ 2494d6df
......@@ -67,6 +67,7 @@ namespace AMDiS
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
......@@ -80,17 +81,17 @@ namespace AMDiS
thread_local std::vector<RowFieldType> rowPartial;
rowPartial.resize(rowShapeGradients.size());
for (std::size_t i = 0; i < rowPartial.size(); ++i) {
rowPartial[i] = jacobian[compTest_][0] * rowShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.cols; ++j)
rowPartial[i] += jacobian[compTest_][j] * rowShapeGradients[i][0][j];
rowPartial[i] = jacobian_mat[compTest_][0] * rowShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
rowPartial[i] += jacobian_mat[compTest_][j] * rowShapeGradients[i][0][j];
}
thread_local std::vector<ColFieldType> colPartial;
colPartial.resize(colShapeGradients.size());
for (std::size_t i = 0; i < colPartial.size(); ++i) {
colPartial[i] = jacobian[compTrial_][0] * colShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.cols; ++j)
colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
colPartial[i] = jacobian_mat[compTrial_][0] * colShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
colPartial[i] += jacobian_mat[compTrial_][j] * colShapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colSize; ++j) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment