From 0a20cfb860831fc03f5a92477b7af7b6c8e830a1 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 13 Jun 2014 09:25:38 +0000
Subject: [PATCH] Do not use the 'hessian' driver of ADOL-C to compute the
Hessian
Instead, only let ADOL-C evaluate the Hessian in the directions
of the vectors of the orthonormal frames. The result is the same,
but we get away with fewer calls to ADOL-C, which, in my Cosserat
shell example, shaves off about 10% of the assembly time.
[[Imported from SVN: r9779]]
---
dune/gfe/localgeodesicfeadolcstiffness.hh | 83 +++++++++++++++--------
1 file changed, 56 insertions(+), 27 deletions(-)
diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index eaaaca33..2f154dab 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -214,22 +214,60 @@ assembleGradientAndHessian(const Entity& element,
/////////////////////////////////////////////////////////////////
// Compute Hessian
/////////////////////////////////////////////////////////////////
- double* rawHessian[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- rawHessian[i] = (double*)malloc((i+1)*sizeof(double));
- hessian(1,nDoubles,xp.data(),rawHessian);
-
- // Copy Hessian into Dune data type
- Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
- for(size_t i=0; i<nDoubles; i++) {
- for (size_t j=0; j<nDoubles; j++) {
- double value = (j<=i) ? rawHessian[i][j] : rawHessian[j][i];
- embeddedHessian[i/embeddedBlocksize][j/embeddedBlocksize][i%embeddedBlocksize][j%embeddedBlocksize] = value;
- }
- }
- for(size_t i=0; i<nDoubles; i++)
- free(rawHessian[i]);
+ // We compute the Hessian of the energy functional using the ADOL-C system.
+ // Since ADOL-C does not know about nonlinear spaces, what we get is actually
+ // the Hessian of a prolongation of the energy functional into the surrounding
+ // Euclidean space. To obtain the Riemannian Hessian from this we apply the
+ // formula described in Absil, Mahoney, Trumpf, "An extrinsic look at the Riemannian Hessian".
+ // This formula consists of two steps:
+ // 1) Remove all entries of the Hessian pertaining to the normal space of the
+ // manifold. In the aforementioned paper this is done by projection onto the
+ // tangent space. Since we want a matrix that is really smaller (but full rank again),
+ // we can achieve the same effect by multiplying the embedded Hessian from the left
+ // and from the right by the matrix of orthonormal frames.
+ // 2) Add a correction involving the Weingarten map.
+ //
+ // This works, and is easy to implement using the ADOL-C "hessian" driver.
+ // However, here we implement a small shortcut. Computing the embedded Hessian and
+ // multiplying one side by the orthonormal frame is the same as evaluating the Hessian
+ // (seen as an operator from R^n to R^n) in the directions of the vectors of the
+ // orthonormal frame. By luck, ADOL-C can compute the evaluations of the Hessian in
+ // a given direction directly (in fact, this is also how the "hessian" driver works).
+ // Since there are less frame vectors than the dimension of the embedding space,
+ // this reinterpretation allows to reduce the number of calls to ADOL-C.
+ // In my Cosserat shell tests this reduced assembly time by about 10%.
+
+ std::vector<Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> > orthonormalFrame(nDofs);
+
+ for (size_t i=0; i<nDofs; i++)
+ orthonormalFrame[i] = localSolution[i].orthonormalFrame();
+
+ Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
+
+ std::vector<double> v(nDoubles);
+ std::vector<double> w(nDoubles);
+
+ std::fill(v.begin(), v.end(), 0.0);
+
+ for (int i=0; i<nDofs; i++)
+ for (int ii=0; ii<blocksize; ii++)
+ {
+ // Evaluate Hessian in the direction of each vector of the orthonormal frame
+ for (size_t k=0; k<embeddedBlocksize; k++)
+ v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
+
+ int rc= 3;
+ MINDEC(rc, hess_vec(1, nDoubles, xp.data(), v.data(), w.data()));
+ if (rc < 0)
+ DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+ for (int j=0; j<nDoubles; j++)
+ embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
+
+ // Make v the null vector again
+ std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
+ }
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
@@ -241,11 +279,6 @@ assembleGradientAndHessian(const Entity& element,
this->A_.setSize(nDofs,nDofs);
- std::vector<Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> > orthonormalFrame(nDofs);
-
- for (size_t i=0; i<nDofs; i++)
- orthonormalFrame[i] = localSolution[i].orthonormalFrame();
-
for (size_t col=0; col<nDofs; col++) {
for (size_t subCol=0; subCol<blocksize; subCol++) {
@@ -253,16 +286,12 @@ assembleGradientAndHessian(const Entity& element,
EmbeddedTangentVector z = orthonormalFrame[col][subCol];
// P_x \partial^2 f z
- std::vector<EmbeddedTangentVector> semiEmbeddedProduct(nDofs);
-
for (size_t row=0; row<nDofs; row++) {
- embeddedHessian[row][col].mv(z,semiEmbeddedProduct[row]);
- //tmp1 = localSolution[row].projectOntoTangentSpace(tmp1);
- TangentVector tmp2;
- orthonormalFrame[row].mv(semiEmbeddedProduct[row],tmp2);
+ TangentVector semiEmbeddedProduct;
+ embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize; subRow++)
- this->A_[row][col][subRow][subCol] = tmp2[subRow];
+ this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
--
GitLab