From 924813d55f37f00589b903bf79109ddd3f4458f6 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 3 Jul 2014 19:15:52 +0000
Subject: [PATCH] Implement vector-mode AD for the Hesse matrix
This almost halves the assembly time for the Hesse matrix. Nice!
Nevertheless I still keep the scalar-mode code around, hidden
behind and #ifdef. Higher-order vector-mode in ADOL-C appears
to be not quite as reliable as the rest, so I want to be able to
easily compare results with the scalar code.
Also, unlike the scalar code, the vector-mode one again computes
the full Hessian in the surrounding Euclidean space. Only afterwards,
that full Hessian is multiplied by the tangent space basis vectors
to obtain a Hesse matrix on the tangent space. To obtain better
efficiency one could have ADOL-C compute the Hessian in the tangent
directions directly. However, ADOL-C sometimes segfaults... when
doing that.
[[Imported from SVN: r9815]]
---
dune/gfe/localgeodesicfeadolcstiffness.hh | 73 +++++++++++++++++++++++
1 file changed, 73 insertions(+)
diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index a1b7a878..3bbbcdd8 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -4,6 +4,7 @@
#include <adolc/adouble.h> // use of active doubles
#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
// gradient(.) and hessian(.)
+#include <adolc/interfaces.h> // use of "Easy to Use" drivers
#include <adolc/taping.h> // use of taping
#include <dune/gfe/adolcnamespaceinjections.hh>
@@ -13,6 +14,8 @@
#include <dune/gfe/localgeodesicfestiffness.hh>
+#define ADOLC_VECTOR_MODE
+
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class GridView, class LocalFiniteElement, class TargetSpace>
@@ -241,6 +244,7 @@ assembleGradientAndHessian(const Entity& element,
for (size_t i=0; i<nDofs; i++)
orthonormalFrame[i] = localSolution[i].orthonormalFrame();
+#ifndef ADOLC_VECTOR_MODE
Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
std::vector<double> v(nDoubles);
@@ -296,6 +300,75 @@ assembleGradientAndHessian(const Entity& element,
}
+#else
+ Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
+
+ int n = nDoubles;
+ //std::cout << "n: " << n << std::endl;
+ int nDirections = nDoubles;
+ double* tangent[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ tangent[i] = (double*)malloc(nDirections*sizeof(double));
+
+ double* rawHessian[nDoubles];
+ for(size_t i=0; i<nDoubles; i++)
+ rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
+
+ //std::cout << "nDirections: " << nDirections << std::endl;
+ for (int i=0; i<n; i++)
+ for (int j=0; j<nDirections; j++)
+ tangent[i][j] = (i==j);
+
+ hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+
+ //std::cout << "Hallo Welt!" << std::endl;
+
+ // Copy Hessian into Dune data type
+ for(size_t i=0; i<nDoubles; i++)
+ for (size_t j=0; j<nDoubles; j++)
+ embeddedHessian[i/embeddedBlocksize][j/embeddedBlocksize][i%embeddedBlocksize][j%embeddedBlocksize] = rawHessian[i][j];
+
+
+ for(size_t i=0; i<nDoubles; i++) {
+ free(rawHessian[i]);
+ free(tangent[i]);
+ }
+ //printmatrix(std::cout, embeddedHessian, "hessian", "--");
+
+ // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
+ // isometrically in a Euclidean space.
+ // For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
+ // at the Riemannian Hessian".
+
+ typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+ typedef typename TargetSpace::TangentVector TangentVector;
+
+ this->A_.setSize(nDofs,nDofs);
+
+ for (size_t col=0; col<nDofs; col++) {
+
+ for (size_t subCol=0; subCol<blocksize; subCol++) {
+
+ 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);
+
+ for (int subRow=0; subRow<blocksize; subRow++)
+ this->A_[row][col][subRow][subCol] = tmp2[subRow];
+ }
+
+ }
+
+ }
+#endif
+
//////////////////////////////////////////////////////////////////////////////////////
// Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f)
--
GitLab