Compute energy gradient together with the Hesse matrix

Because the gradient is needed anyway when computing the Hesse matrix.
Not computing it separately a second time shaves off about 8% wall time
of each iteration.

[[Imported from SVN: r9555]]

dune/gfe/geodesicfeassembler.hh

@@ -43,11 +43,15 @@ public:
           localStiffness_(localStiffness)
     {}
 
-    /** \brief Assemble the tangent stiffness matrix
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This is more efficient than computing them separately, because you need the gradient
+     * anyway to compute the Riemannian Hessian.
      */
-    virtual void assembleMatrix(const std::vector<TargetSpace>& sol,
-                                Dune::BCRSMatrix<MatrixBlock>& matrix,
-                                bool computeOccupationPattern=true) const;
+    virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                                            Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                                            Dune::BCRSMatrix<MatrixBlock>& hessian,
+                                            bool computeOccupationPattern=true) const;
 
     /** \brief Assemble the gradient */
     virtual void assembleGradient(const std::vector<TargetSpace>& sol,
@@ -97,19 +101,23 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
 
 template <class Basis, class TargetSpace>
 void GeodesicFEAssembler<Basis,TargetSpace>::
-assembleMatrix(const std::vector<TargetSpace>& sol,
-               Dune::BCRSMatrix<MatrixBlock>& matrix,
-               bool computeOccupationPattern) const
+assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                           Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                           Dune::BCRSMatrix<MatrixBlock>& hessian,
+                           bool computeOccupationPattern) const
 {
     if (computeOccupationPattern) {
 
         Dune::MatrixIndexSet neighborsPerVertex;
         getNeighborsPerVertex(neighborsPerVertex);
-        neighborsPerVertex.exportIdx(matrix);
+        neighborsPerVertex.exportIdx(hessian);
 
     }
 
-    matrix = 0;
+    hessian = 0;
+
+    gradient.resize(sol.size());
+    gradient = 0;
 
     ElementIterator it    = basis_.getGridView().template begin<0>();
     ElementIterator endit = basis_.getGridView().template end<0>  ();
@@ -124,8 +132,10 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
         for (int i=0; i<numOfBaseFct; i++)
             localSolution[i] = sol[basis_.index(*it,i)];
 
-        // setup matrix
-        localStiffness_->assembleHessian(*it, basis_.getLocalFiniteElement(*it), localSolution);
+        std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
+
+        // setup local matrix and gradient
+        localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
 
         // Add element matrix to global stiffness matrix
         for(int i=0; i<numOfBaseFct; i++) {
@@ -135,11 +145,15 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
             for (int j=0; j<numOfBaseFct; j++ ) {
 
                 int col = basis_.index(*it,j);
-                matrix[row][col] += localStiffness_->A_[i][j];
+                hessian[row][col] += localStiffness_->A_[i][j];
 
             }
         }
 
+        // Add local gradient to global gradient
+        for (int i=0; i<numOfBaseFct; i++)
+            gradient[basis_.index(*it,i)] += localGradient[i];
+
     }
 
 }

dune/gfe/localgeodesicfeadolcstiffness.hh

@@ -61,9 +61,10 @@ public:
 
     This uses the automatic differentiation toolbox ADOL_C.
     */
-    virtual void assembleHessian(const Entity& e,
+    virtual void assembleGradientAndHessian(const Entity& e,
                          const LocalFiniteElement& localFiniteElement,
-                         const std::vector<TargetSpace>& localSolution);
+                         const std::vector<TargetSpace>& localSolution,
+                         std::vector<typename TargetSpace::TangentVector>& localGradient);
 
     const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* localEnergy_;
 
@@ -162,13 +163,16 @@ assembleGradient(const Entity& element,
 
 
 // ///////////////////////////////////////////////////////////
-//   Compute gradient by finite-difference approximation
+//   Compute gradient and Hessian together
+//   To compute the Hessian we need to compute the gradient anyway, so we may
+//   as well return it.  This saves assembly time.
 // ///////////////////////////////////////////////////////////
 template <class GridType, class LocalFiniteElement, class TargetSpace>
 void LocalGeodesicFEADOLCStiffness<GridType, LocalFiniteElement, TargetSpace>::
-assembleHessian(const Entity& element,
+assembleGradientAndHessian(const Entity& element,
                 const LocalFiniteElement& localFiniteElement,
-                const std::vector<TargetSpace>& localSolution)
+                const std::vector<TargetSpace>& localSolution,
+                std::vector<typename TargetSpace::TangentVector>& localGradient)
 {
     // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
     energy(element, localFiniteElement, localSolution);
@@ -199,6 +203,12 @@ assembleHessian(const Entity& element,
         for (size_t j=0; j<embeddedBlocksize; j++)
             localEmbeddedGradient[i][j] = g[idx++];
 
+    // Express gradient in local coordinate system
+    for (size_t i=0; i<nDofs; i++) {
+        Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
+        orthonormalFrame.mv(localEmbeddedGradient[i],localGradient[i]);
+    }
+
     /////////////////////////////////////////////////////////////////
     // Compute Hessian
     /////////////////////////////////////////////////////////////////

dune/gfe/localgeodesicfestiffness.hh

@@ -39,9 +39,10 @@ public:
     We compute that using a finite difference approximation.
 
     */
-    virtual void assembleHessian(const Entity& e,
+    virtual void assembleGradientAndHessian(const Entity& e,
                                  const LocalFiniteElement& localFiniteElement,
-                                 const std::vector<TargetSpace>& localSolution);
+                                 const std::vector<TargetSpace>& localSolution,
+                                 std::vector<typename TargetSpace::TangentVector>& localGradient);
 
     /** \brief Compute the energy at the current configuration */
     virtual RT energy (const Entity& e,
@@ -114,9 +115,10 @@ assembleGradient(const Entity& element,
 // ///////////////////////////////////////////////////////////
 template <class GridType, class LocalFiniteElement, class TargetSpace>
 void LocalGeodesicFEStiffness<GridType, LocalFiniteElement, TargetSpace>::
-assembleHessian(const Entity& element,
+assembleGradientAndHessian(const Entity& element,
                 const LocalFiniteElement& localFiniteElement,
-                const std::vector<TargetSpace>& localSolution)
+                const std::vector<TargetSpace>& localSolution,
+                std::vector<typename TargetSpace::TangentVector>& localGradient)
 {
     // Number of degrees of freedom for this element
     size_t nDofs = localSolution.size();
@@ -161,9 +163,22 @@ assembleHessian(const Entity& element,
 
     }
 
-    // finite-difference approximation
-    // we loop over the lower left triangular half of the matrix.
-    // The other half follows from symmetry
+    //////////////////////////////////////////////////////////////
+    //   Compute gradient by finite-difference approximation
+    //////////////////////////////////////////////////////////////
+
+    localGradient.resize(localSolution.size());
+
+    for (size_t i=0; i<localSolution.size(); i++)
+        for (int j=0; j<blocksize; j++)
+            localGradient[i][j] = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    //   Compute Riemannian Hesse matrix by finite-difference approximation.
+    //   We loop over the lower left triangular half of the matrix.
+    //   The other half follows from symmetry.
+    ///////////////////////////////////////////////////////////////////////////
     //#pragma omp parallel for schedule (dynamic)
     for (size_t i=0; i<localSolution.size(); i++) {
         for (size_t i2=0; i2<blocksize; i2++) {

dune/gfe/riemanniantrsolver.cc

@@ -240,19 +240,19 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
         if (recomputeGradientHessian) {
 
-            assembler_->assembleGradient(x_, rhs);
+            assembler_->assembleGradientAndHessian(x_,
+                                                   rhs,
+                                                   *hessianMatrix_,
+                                                   i==0    // assemble occupation pattern only for the first call
+                                                   );
+
             rhs *= -1;        // The right hand side is the _negative_ gradient
-            if (this->verbosity_ == Solver::FULL)
-              std::cout << "gradient assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
-            gradientTimer.reset();
 
-            assembler_->assembleMatrix(x_,
-                                       *hessianMatrix_,
-                                       i==0    // assemble occupation pattern only for the first call
-                                       );
             if (this->verbosity_ == Solver::FULL)
-              std::cout << "hessian assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+              std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
             recomputeGradientHessian = false;
+            
         }
 
 /*        std::cout << "rhs:\n" << rhs << std::endl;