Simplify using smoother as coarse solver

Is no base solver is provided, the MutliGridStep does now
do pre- and postsmoothing on the coarse level as it does
on any other level.

dune/solvers/iterationsteps/multigridstep.cc

@@ -146,49 +146,52 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
     //   Set up base solver
     // /////////////////////////////////////////////
 
-    if (this->basesolver_ == NULL)
-        DUNE_THROW(SolverError, "You have not provided a base solver!");
+    if (basesolver_)
+    {
+        if (this->basesolver_ == NULL)
+            DUNE_THROW(SolverError, "You have not provided a base solver!");
 
-    // If the base solver can ignore dofs give it the ignoreNodes field
-    if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_))
-        dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
+        // If the base solver can ignore dofs give it the ignoreNodes field
+        if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_))
+            dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
 
-    typedef ::LoopSolver<VectorType> DuneSolversLoopSolver;
+        typedef ::LoopSolver<VectorType> DuneSolversLoopSolver;
 
-    if (typeid(*this->basesolver_) == typeid(DuneSolversLoopSolver)) {
+        if (typeid(*this->basesolver_) == typeid(DuneSolversLoopSolver)) {
 
-        DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_);
+            DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_);
 
-        typedef LinearIterationStep<MatrixType, VectorType> SmootherType;
-        assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));
+            typedef LinearIterationStep<MatrixType, VectorType> SmootherType;
+            assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));
 
-        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
-        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
+            dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
+            dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
 
-    }
-    // TODO: The following two #if-clauses do the exactly the same thing: they call setProblem
-    // However, we cannot write the code generically because there is no general abstraction
-    // for solvers that I can call 'setProblem' for.
+        }
+        // TODO: The following two #if-clauses do the exactly the same thing: they call setProblem
+        // However, we cannot write the code generically because there is no general abstraction
+        // for solvers that I can call 'setProblem' for.
 #if HAVE_IPOPT
-    else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {
+        else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {
 
-        QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
+            QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
 
-        ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
-    }
+            ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
+        }
 #endif
 #if HAVE_UMFPACK
-    else if (typeid(*this->basesolver_) == typeid(Dune::Solvers::UMFPackSolver<MatrixType,VectorType>)) {
+        else if (typeid(*this->basesolver_) == typeid(Dune::Solvers::UMFPackSolver<MatrixType,VectorType>)) {
 
-        Dune::Solvers::UMFPackSolver<MatrixType,VectorType>* umfpackBaseSolver
-            = dynamic_cast<Dune::Solvers::UMFPackSolver<MatrixType,VectorType>*> (this->basesolver_);
+            Dune::Solvers::UMFPackSolver<MatrixType,VectorType>* umfpackBaseSolver
+                = dynamic_cast<Dune::Solvers::UMFPackSolver<MatrixType,VectorType>*> (this->basesolver_);
 
-        umfpackBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
-    }
+            umfpackBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
+        }
 #endif
-else {
-        DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
-                   << " as a base solver in a MultigridStep!");
+        else {
+            DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
+                       << " as a base solver in a MultigridStep!");
+        }
     }
 }
 
@@ -226,7 +229,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
     std::vector<VectorType>& rhs = this->rhsHierarchy_;
 
     // Solve directly if we're looking at the coarse problem
-    if (level == 0) {
+    if ((level == 0) and (basesolver_)) {
         basesolver_->solve();
         return;
     }
@@ -239,38 +242,41 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
     for (int i=0; i<this->nu1_; i++)
         presmoother_[level]->iterate();
 
-    // /////////////////////////
-    // Restriction
+    if (level>0)
+    {
+        // /////////////////////////
+        // Restriction
 
-    // compute residual
-    // fineResidual = rhs[level] - mat[level] * x[level];
-    VectorType fineResidual = rhs[level];
-    mat[level]->mmv(*x[level], fineResidual);
+        // compute residual
+        // fineResidual = rhs[level] - mat[level] * x[level];
+        VectorType fineResidual = rhs[level];
+        mat[level]->mmv(*x[level], fineResidual);
 
-    // restrict residual
-    this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
+        // restrict residual
+        this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
 
 
-    // Set Dirichlet values.
-    GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
+        // Set Dirichlet values.
+        GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
 
-    // Choose all zeros as the initial correction
-    *x[level-1] = 0;
+        // Choose all zeros as the initial correction
+        *x[level-1] = 0;
 
-    // ///////////////////////////////////////
-    // Recursively solve the coarser system
-    level--;
-    for (int i=0; i<mu_; i++)
-        iterate();
-    level++;
+        // ///////////////////////////////////////
+        // Recursively solve the coarser system
+        level--;
+        for (int i=0; i<mu_; i++)
+            iterate();
+        level++;
 
-    // ////////////////////////////////////////
-    // Prolong
+        // ////////////////////////////////////////
+        // Prolong
 
-    // add correction to the presmoothed solution
-    VectorType tmp;
-    this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
-    *x[level] += tmp;
+        // add correction to the presmoothed solution
+        VectorType tmp;
+        this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
+        *x[level] += tmp;
+    }
 
     // Postsmoothing
     postsmoother_[level]->setProblem(*(mat[level]), *x[level], rhs[level]);

[akbib]

Can't one remove these two lines now?

[graeser]

Done.