From 66383a0f068049e5f45d172f8c39c59b4db69322 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 20 Apr 2009 13:12:00 +0000
Subject: [PATCH] bugfix: compute matrix occupation pattern on demand
[[Imported from SVN: r4063]]
---
src/geodesicfeassembler.hh | 17 ++++++++++++-----
src/riemanniantrsolver.cc | 18 ++++++------------
2 files changed, 18 insertions(+), 17 deletions(-)
diff --git a/src/geodesicfeassembler.hh b/src/geodesicfeassembler.hh
index f6573046..10a06071 100644
--- a/src/geodesicfeassembler.hh
+++ b/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>();
diff --git a/src/riemanniantrsolver.cc b/src/riemanniantrsolver.cc
index 0e386667..9c362fe1 100644
--- a/src/riemanniantrsolver.cc
+++ b/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)
--
GitLab