From 28674731ab5f2bb47e2a0a43a5f7409c197b4934 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 23 Apr 2014 10:20:53 +0000
Subject: [PATCH] Hand-code one specific matrix-matrix multiplication
This has two advantages:
- The two matrices use different number types (adouble vs. double),
but the generic multiplication in dune/istl/matrix.hh doesn't support
that. Hence before this change I had to patch dune-istl.
- It avoids one extra copying operation
[[Imported from SVN: r9703]]
---
dune/gfe/localgeodesicfefunction.hh | 20 +++++++++++++-------
1 file changed, 13 insertions(+), 7 deletions(-)
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 9316ff48..1ac7a335 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -245,19 +245,25 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
// ////////////////////////////////////////////////////////////////////////
// the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
- std::vector<Dune::FieldMatrix<ctype,1,dim> > BNested(coefficients_.size());
- localFiniteElement_.localBasis().evaluateJacobian(local, BNested);
- Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > B(coefficients_.size(), dim);
- for (size_t i=0; i<coefficients_.size(); i++)
- for (size_t j=0; j<dim; j++)
- B[i][j] = BNested[i][0][j];
+ std::vector<Dune::FieldMatrix<ctype,1,dim> > B(coefficients_.size());
+ localFiniteElement_.localBasis().evaluateJacobian(local, B);
// compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
Dune::Matrix<Dune::FieldMatrix<RT,1,1> > dFdw = computeDFdw(q);
dFdw *= -1;
// multiply the two previous matrices: the result is the right hand side
- Dune::Matrix<Dune::FieldMatrix<RT,1,1> > RHS = dFdw * B;
+ // RHS = dFdw * B;
+ // We need to write this multiplication by hand, because B is actually an (#coefficients) times 1 times dim matrix,
+ // and we need to ignore the middle index.
+ Dune::Matrix<Dune::FieldMatrix<RT,1,1> > RHS(dFdw.N(), dim);
+
+ for (size_t i=0; i<RHS.N(); i++)
+ for (size_t j=0; j<RHS.M(); j++) {
+ RHS[i][j] = 0;
+ for (size_t k=0; k<dFdw.M(); k++)
+ RHS[i][j] += dFdw[i][k]*B[k][0][j];
+ }
// the actual system matrix
std::vector<Dune::FieldVector<ctype,1> > w;
--
GitLab