ISTLSeqILU0Step: Overhaul. Now works again (mostly)

IgnoreNodes are still not handled.

Judging by all the errors that have kept this class from compiling for
quite a while, I doubt that the interface changes made in this commit
will affect anyone.

dune/solvers/iterationsteps/istlseqilu0step.hh

@@ -3,7 +3,7 @@
 #ifndef ISTL_SEQILU0_STEP_HH
 #define ISTL_SEQILU0_STEP_HH
 
-/** \file 
+/** \file
     \brief A wrapper class for the ISTL SeqILU0 implementation
  */
 
@@ -12,69 +12,56 @@
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
 
 #include <dune/istl/preconditioners.hh>
-#include <dune/istl/operators.hh>
 
 /** \brief A wrapper class for the ISTL SeqILU0 implementation
  */
 template <class MatrixType, class VectorType>
 class ISTLSeqILU0Step
+    // FIXME: ignoreNodes are not handled
     : public LinearIterationStep<MatrixType, VectorType>
 {
-    typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
-
     typedef Dune::SeqILU0<MatrixType,VectorType,VectorType> SeqILU0;
 
 public:
-
-    /** \brief Constructor which initializes and sets up an algebraic hierarchy
-     */
-    ISTLSeqILU0Step (const MatrixType* stiffnessMatrix, 
-                     VectorType& x, 
-                     VectorType& rhs,
-                     double relaxationFactor)
+    ISTLSeqILU0Step (double relaxationFactor)
         : relaxationFactor_(relaxationFactor)
-    {
-        setProblem(*stiffnessMatrix, x, rhs);
-
-        seqILU0_ = std::auto_ptr<SeqILU0>(new SeqILU0(*stiffnessMatrix, relaxationFactor));
-        seqILU0_->pre(*this->x_,*this->rhs_);
-    }
-
-    /** \brief Initialize the iteration step but not carry out any preprocessing */
-    ISTLSeqILU0Step (const MatrixType* stiffnessMatrix, 
-             VectorType& x, 
-             VectorType& rhs)
-    {
-        setProblem(*stiffnessMatrix, x, rhs);
-    }
+    {}
 
     /** \brief Sets up an algebraic hierarchy
      */
-    virtual void preprocess();
+    virtual void preprocess() override {
+        if (needReconstruction_) {
+            seqILU0_ = std::make_unique<SeqILU0>(*this->mat_, relaxationFactor_);
+            needReconstruction_ = false;
+        }
+        // Note: as of now, pre() is a dummy
+        mutable_rhs = *this->rhs_;
+        seqILU0_->pre(*this->x_,mutable_rhs);
+    }
 
     /** \brief Perform one iteration */
-    virtual void iterate();
+    virtual void iterate() override {
+        mutable_rhs = *this->rhs_;
+        seqILU0_->apply(*this->x_, mutable_rhs);
+    }
+
+    virtual void setMatrix(const MatrixType& mat) override {
+        this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+        needReconstruction_ = true;
+    }
 
 private:
 
     /** \brief The dune-istl sequential ILU0 preconditioner */
-    std::auto_ptr<SeqILU0> seqILU0_;
+    std::unique_ptr<SeqILU0> seqILU0_;
 
-    /** \brief Some magic relaxation factor */
-    double relaxationFactor_;
-};
+    VectorType mutable_rhs;
 
-template <class MatrixType, class VectorType>
-void ISTLSeqILU0Step<MatrixType,VectorType>::preprocess()
-{
-    seqILU0_ = std::auto_ptr<SeqILU0>(new SeqILU0(*this->mat_, relaxationFactor_));
-    seqILU0_->pre(*this->x_,*this->rhs_);
-}
-
-template <class MatrixType, class VectorType>
-void ISTLSeqILU0Step<MatrixType,VectorType>::iterate()
-{
-    seqILU0_->apply(*this->x_,*this->rhs_);
-}
+    double const relaxationFactor_;
 
+    /** \brief Calling the Dune::SeqILU0 constructor is expensive; we
+               thus remember if setMatrix has been called. And only
+               call the constructor if necessary */
+    bool needReconstruction_;
+};
 #endif