bugfix: compute matrix occupation pattern on demand

[[Imported from SVN: r4063]]

src/geodesicfeassembler.hh

@@ -43,7 +43,8 @@ public:
     /** \brief Assemble the tangent stiffness matrix
      */
     virtual void assembleMatrix(const std::vector<TargetSpace>& sol,
-                        Dune::BCRSMatrix<MatrixBlock>& matrix) const;
+                                Dune::BCRSMatrix<MatrixBlock>& matrix,
+                                bool computeOccupationPattern=true) const;
     
     /** \brief Assemble the gradient */
     virtual void assembleGradient(const std::vector<TargetSpace>& sol,
@@ -94,13 +95,19 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
 template <class GridView, class TargetSpace>
 void GeodesicFEAssembler<GridView,TargetSpace>::
 assembleMatrix(const std::vector<TargetSpace>& sol,
-               Dune::BCRSMatrix<MatrixBlock>& matrix) const
+               Dune::BCRSMatrix<MatrixBlock>& matrix,
+               bool computeOccupationPattern) const
 {
     const typename GridView::IndexSet& indexSet = gridView_.indexSet();
 
-    Dune::MatrixIndexSet neighborsPerVertex;
-    getNeighborsPerVertex(neighborsPerVertex);
-    
+    if (computeOccupationPattern) {
+
+        Dune::MatrixIndexSet neighborsPerVertex;
+        getNeighborsPerVertex(neighborsPerVertex);
+        neighborsPerVertex.exportIdx(matrix);
+
+    }
+
     matrix = 0;
     
     ElementIterator it    = gridView_.template begin<0>();

src/riemanniantrsolver.cc

@@ -183,7 +183,10 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
         corr = 0;
 
         assembler_->assembleGradient(x_, rhs);
-        assembler_->assembleMatrix(x_, *hessianMatrix_);
+        assembler_->assembleMatrix(x_, 
+                                   *hessianMatrix_, 
+                                   i==0    // assemble occupation pattern only for the first call
+                                   );
 
         //gradientFDCheck(x_, rhs, *rodAssembler_);
         //hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_);
@@ -278,17 +281,8 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
             
         }
 
-        if (this->verbosity_ == NumProc::FULL) {
-            double translationMax = 0;
-            double rotationMax    = 0;
-            for (size_t j=0; j<corr.size(); j++) {
-                for (int k=0; k<3; k++) {
-                    translationMax = std::max(translationMax, corr[j][k]);
-                    rotationMax    = std::max(rotationMax, corr[j][k+3]);
-                }
-            }
-            printf("infinity norm of the correction: %g %g\n", translationMax, rotationMax);
-        }
+        if (this->verbosity_ == NumProc::FULL)
+            std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
 
         if (corr.infinity_norm() < this->tolerance_) {
             if (this->verbosity_ == NumProc::FULL)