From aa9490a86514bb9239ae0da6b570f127067887d5 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 7 Jun 2011 08:03:08 +0000
Subject: [PATCH] Remove unnecessary multiplications with the orthonormal frame
matrices
[[Imported from SVN: r7378]]
---
dune/gfe/localgeodesicfestiffness.hh | 88 ++++++++++++++++++----------
1 file changed, 57 insertions(+), 31 deletions(-)
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 36c60d71..7a9995d1 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -380,6 +380,10 @@ public:
/** \brief Assemble the local stiffness matrix at the current position
This default implementation used finite-difference approximations to compute the second derivatives
+
+ The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
+ 'Optimization algorithms on matrix manifolds', page 107
+
*/
virtual void assembleHessian(const Entity& e,
const std::vector<TargetSpace>& localSolution);
@@ -454,7 +458,7 @@ public:
}
}
-
+
// assembled data
Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_;
@@ -496,43 +500,65 @@ assembleHessian(const Entity& element,
A_ = 0;
- // first compute the Hessian in the embedding space
- Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
- std::vector<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > embeddedGradient;
+
+
+ const double eps = 1e-4;
+
+ std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
+ for (size_t i=0; i<B.size(); i++)
+ B[i] = localSolution[i].orthonormalFrame();
+
+ // finite-difference approximation
for (size_t i=0; i<localSolution.size(); i++) {
+ for (size_t i2=0; i2<blocksize; i2++) {
+ for (size_t j=0; j<localSolution.size(); j++) {
+ for (size_t j2=0; j2<blocksize; j2++) {
- embeddedGradientOfEmbeddedGradient(element,localSolution, i, embeddedGradient);
-
- for (size_t j=0; j<localSolution.size(); j++)
- embeddedHessian[i][j] = embeddedGradient[j];
-
- }
-
- // transform to local tangent space bases
- std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > orthonormalFrames(nDofs);
- std::vector<Dune::FieldMatrix<double,embeddedBlocksize,blocksize> > orthonormalFramesTransposed(nDofs);
-
- for (size_t i=0; i<nDofs; i++) {
- orthonormalFrames[i] = localSolution[i].orthonormalFrame();
-
- for (int j=0; j<embeddedBlocksize; j++)
- for (int k=0; k<blocksize; k++)
- orthonormalFramesTransposed[i][j][k] = orthonormalFrames[i][k][j];
-
- }
-
- for (size_t i=0; i<nDofs; i++)
- for (size_t j=0; j<nDofs; j++) {
-
- Dune::FieldMatrix<double,blocksize,embeddedBlocksize> tmp;
- Dune::FMatrixHelp::multMatrix(orthonormalFrames[i],embeddedHessian[i][j],tmp);
- A_[i][j] = tmp.rightmultiplyany(orthonormalFramesTransposed[j]);
+ Dune::FieldVector<double,embeddedBlocksize> epsXi = B[i][i2]; epsXi *= eps;
+ Dune::FieldVector<double,embeddedBlocksize> epsEta = B[j][j2]; epsEta *= eps;
+
+ Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi; minusEpsXi *= -1;
+ Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta; minusEpsEta *= -1;
+
+ std::vector<TargetSpace> forwardSolutionXiEta = localSolution;
+ std::vector<TargetSpace> forwardSolutionXi = localSolution;
+ std::vector<TargetSpace> forwardSolutionEta = localSolution;
+ std::vector<TargetSpace> backwardSolutionXiEta = localSolution;
+ std::vector<TargetSpace> backwardSolutionXi = localSolution;
+ std::vector<TargetSpace> backwardSolutionEta = localSolution;
+ if (i==j)
+ forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
+ else {
+ forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi);
+ forwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],epsEta);
+ }
+ forwardSolutionXi[i] = TargetSpace::exp(localSolution[i],epsXi);
+ forwardSolutionEta[j] = TargetSpace::exp(localSolution[j],epsEta);
+
+ if (i==j)
+ backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi+minusEpsEta);
+ else {
+ backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
+ backwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
+ }
+
+ backwardSolutionXi[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
+ backwardSolutionEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
+ double forwardValue = energy(element, forwardSolutionXiEta) - energy(element, forwardSolutionXi) - energy(element, forwardSolutionEta);
+ double centerValue = -energy(element, localSolution);
+ double backwardValue = energy(element, backwardSolutionXiEta) - energy(element, backwardSolutionXi) - energy(element, backwardSolutionEta);
+
+ A_[i][j][i2][j2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
+
+ }
+ }
}
-
+ }
+
}
#endif
--
GitLab