Store internal matrix hierarchy as a vector of shared_ptrs

Instead of a vector of raw pointers.  This fixes a memory
which bites you when you call 'preprocess' many many times.

This patch does not change any interface.

[[Imported from SVN: r6426]]

dune/solvers/iterationsteps/lineariterationstep.hh

@@ -21,7 +21,7 @@ public:
     
     //! Constructor being given linear operator, solution and right hand side
     LinearIterationStep(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
-        mat_     = &mat;
+        mat_     = Dune::stackobject_to_shared_ptr(mat);
         this->x_ = &x;
         rhs_     = &rhs;
     }
@@ -30,7 +30,7 @@ public:
     virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
         this->x_ = &x;
         rhs_     = &rhs;
-        mat_     = &mat;
+        mat_     = Dune::stackobject_to_shared_ptr(mat);
     }
     
     //! Do the actual iteration
@@ -40,7 +40,7 @@ public:
     virtual VectorType getSol() = 0;
     
     //! Return linear operator
-    virtual const MatrixType* getMatrix() {return mat_;}
+    virtual const MatrixType* getMatrix() {return mat_.get();}
     
     /** \brief Checks whether all relevant member variables are set
      * \exception SolverError if the iteration step is not set up properly
@@ -60,7 +60,7 @@ public:
     const VectorType* rhs_;
     
     //! The linear operator
-    const MatrixType* mat_;
+    Dune::shared_ptr<const MatrixType> mat_;
     
 };
 

dune/solvers/iterationsteps/multigridstep.cc

@@ -21,7 +21,7 @@ inline
 const MatrixType* MultigridStep<MatrixType, VectorType, BitVectorType>::
 getMatrix()
 {
-    return this->mat_[this->mat_.size()-1];
+    return this->mat_.back().get();
 }
 
 template<class MatrixType, class VectorType, class BitVectorType>
@@ -60,20 +60,14 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
     // /////////////////////////////////////////////////////
     for (int i=this->numLevels_-2; i>=0; i--) {
 
-        // Delete coarse grid matrix if it exists already
-        if (this->mat_[i]) {
-            delete(this->mat_[i]);
-            this->mat_[i] = NULL;
-        }
-
-        this->mat_[i] = new MatrixType;
+        this->mat_[i] = Dune::shared_ptr<MatrixType>(new MatrixType);
 
         // Compute which entries are present in the (sparse) coarse grid stiffness
         // matrices. 
-        this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->mat_[i+1], *(const_cast<MatrixType*>(this->mat_[i])));
+        this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
 
         // setup matrix
-        this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *(const_cast<MatrixType*>(this->mat_[i])));
+        this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
 
         // Set solution vector sizes for the lower levels
         GenericVector::resize(this->x_[i],*this->mat_[i]);
@@ -167,7 +161,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
     int& level = this->level_;
 
     // Define references just for ease of notation
-    std::vector<const MatrixType*>& mat = this->mat_;
+    std::vector<Dune::shared_ptr<const MatrixType> >& mat = this->mat_;
     std::vector<VectorType>& x   = this->x_;
     std::vector<VectorType>& rhs = this->rhs_;
 

dune/solvers/iterationsteps/multigridstep.hh

@@ -86,12 +86,8 @@
 
 
         virtual ~MultigridStep() {
-            // Delete all matrices but the top one, which has been supplied from the outside
             for (size_t i=0; i<mat_.size()-1; i++)
-            {
-                delete(mat_[i]);
                 delete(ignoreNodesHierarchy_[i]);
-            }
         }
 
         virtual void setProblem(const MatrixType& mat, 
@@ -107,10 +103,6 @@
         {
             numLevels_ = numLevels;
 
-            for (int i=0; i<int(mat_.size())-1; i++)
-                if (mat_[i])
-                    delete(mat_[i]);
-
             for (int i=0; i<int(ignoreNodesHierarchy_.size())-1; i++)
                 if (ignoreNodesHierarchy_[i])
                     delete(ignoreNodesHierarchy_[i]);
@@ -137,7 +129,7 @@
         {
             level_ = numLevels_-1;
 
-            mat_[level_] = &mat;
+            mat_[level_] = Dune::stackobject_to_shared_ptr(mat);
             x_[level_]   = x;
             rhs_[level_] = rhs;
 
@@ -225,7 +217,7 @@
 
     public:
         //! The linear operators on each level
-        std::vector<const MatrixType*> mat_;
+        std::vector<Dune::shared_ptr<const MatrixType> > mat_;
 
     protected:
         //! Flags specifying the dirichlet nodes on each level