dune/solvers/iterationsteps/projectedblockgsstep.cc

-template<class OperatorType, class DiscFuncType>
-inline
-DiscFuncType ProjectedBlockGSStep<OperatorType, DiscFuncType>::getSol()
-{
-    return *(this->x_);
-}
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 
-template<class OperatorType, class DiscFuncType>
+template<class MatrixType, class VectorType>
 inline
-void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
+void ProjectedBlockGSStep<MatrixType, VectorType>::iterate()
 {
+    assert(hasObstacle_!=nullptr);
+
     if (hasObstacle_->size()!= (unsigned int)this->x_->size())
-        DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size() 
+        DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size()
                    << ") doesn't match solution vector (" << this->x_->size() << ")");
 
-    const OperatorType& mat = *this->mat_;
+    const MatrixType& mat = *this->mat_;
 
     for (size_t i=0; i<this->x_->size(); i++) {
-
-        /** \todo Handle more general boundary conditions */
-        if ((*this->ignoreNodes_)[i][0])
+        if (this->ignore()[i].all())
             continue;
+        if (this->ignore()[i].any())
+          DUNE_THROW(Dune::NotImplemented,
+                     "Individual blocks must be either ignored completely, or not at all");
 
         bool zeroDiagonal = false;
-        for (int j=0; j<BlockSize; j++) {
+        for (size_t j=0; j<BlockSize; j++) {
             // When using this solver as part of a truncated multigrid solver,
             // the diagonal entries of the matrix may get completely truncated
             // away.  In this case we just do nothing here.
@@ -35,44 +34,43 @@ void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
         VectorBlock r;
         this->residual(i, r);
 
-        // Compute x_i += A_{i,i}^{-1} r[i] 
+        // Compute x_i += A_{i,i}^{-1} r[i]
         VectorBlock v;
         VectorBlock& x = (*this->x_)[i];
 
 
-        if ((*hasObstacle_)[i] == false) {
+        if ((*hasObstacle_)[i].none()) {
 
             // No obstacle --> solve linear problem
             mat[i][i].solve(v, r);
-        
+
         } else {
-            
+
             // Solve the local constraint minimization problem
             // We use a projected Gauss-Seidel, for lack of anything better
-            BoxConstraint<field_type,BlockSize> defectObstacle = (*obstacles_)[i];
+            assert(obstacles_!=nullptr);
+
+            Obstacle defectObstacle = (*obstacles_)[i];
             defectObstacle -= x;
 
             // Initial correction
             v = 0;
 
-            for (int j=0; j< ((BlockSize==1) ? 1 : 20); j++) {
+            for (size_t j=0; j< ((BlockSize==1) ? 1 : 20); j++) {
 
-                for (int k=0; k<BlockSize; k++) {
+                for (size_t k=0; k<BlockSize; k++) {
 
                     // Compute residual
-                    field_type sr = 0;
-                    for (int l=0; l<BlockSize; l++)
+                    Field sr = 0;
+                    for (size_t l=0; l<BlockSize; l++)
                         sr += mat[i][i][k][l] * v[l];
 
                     sr = r[k] - sr;
                     v[k] += sr / mat[i][i][k][k];
 
                     // Project
-                    if (v[k] < defectObstacle.lower(k))
-                        v[k] = defectObstacle.lower(k);
-                    else if (v[k] > defectObstacle.upper(k))
-                        v[k] = defectObstacle.upper(k);
-
+                    if ((*hasObstacle_)[i][k])
+                        v[k] = defectObstacle[k].projectIn(v[k]);
                 }
 
             }

dune/solvers/iterationsteps/projectedblockgsstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
 #define DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
 
 #include <vector>
 
-#include "blockgsstep.hh"
+#include <dune/matrix-vector/axpy.hh>
+
 #include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
 
-template<class OperatorType, class DiscFuncType>
-class ProjectedBlockGSStep : public BlockGSStep<OperatorType, DiscFuncType>
+template<class MatrixType, class VectorType>
+class ProjectedBlockGSStep : public LinearIterationStep<MatrixType, VectorType>
 {
-    
-    typedef typename DiscFuncType::block_type VectorBlock;
-    
-    typedef typename DiscFuncType::field_type field_type;
+    using Base = LinearIterationStep<MatrixType, VectorType>;
+    using VectorBlock = typename VectorType::block_type;
+    using Field = typename VectorType::field_type;
 
     enum {BlockSize = VectorBlock::dimension};
-    
+
 public:
-    
     //! Default constructor.  Doesn't init anything
     ProjectedBlockGSStep() {}
-    
+
     //! Constructor with a linear problem
-    ProjectedBlockGSStep(const OperatorType& mat, DiscFuncType& x, const DiscFuncType& rhs)
-        : BlockGSStep<OperatorType,DiscFuncType>(mat, x, rhs)
+    ProjectedBlockGSStep(const MatrixType& mat, VectorType& x, const VectorType& rhs)
+        : Base(mat, x, rhs)
     {}
 
     //! Perform one iteration
     virtual void iterate();
-    
-    virtual DiscFuncType getSol();
-    
-    Dune::BitSetVector<1>* hasObstacle_;
-    
-    const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
+
+    /** \brief Compute the residual of the current iterate of a (block) degree of freedom
+     *
+     *  \param index Global index of the (block) degree of freedom
+     *  \param r Write residual in this vector
+     */
+    virtual void residual(int index, VectorBlock& r) const {
+        r = (*this->rhs_)[index];
+
+        const auto& row = (*this->mat_)[index];
+        for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
+            // r_i -= A_ij x_j
+            Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
+        }
+    }
+
+    using HasObstacle = Dune::BitSetVector<BlockSize>;
+    using Obstacle = BoxConstraint<Field, BlockSize>;
+    const HasObstacle* hasObstacle_ = nullptr;
+    const std::vector<Obstacle>* obstacles_ = nullptr;
 };
 
 #include "projectedblockgsstep.cc"

dune/solvers/iterationsteps/projectedgradientstep.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#include <dune/solvers/computeenergy.hh>
+
+template<class MatrixType, class VectorType>
+inline void ProjectedGradientStep<MatrixType, VectorType>::computeGeneralizedCP(const VectorType& dir)
+{
+
+    // Compute the first local minimum along the piecewise path of the projected direction
+    VectorType projGrad = dir;
+
+    ////////////////////////////////////////////////
+    // First compute the projection of the direction
+    // ////////////////////////////////////////////
+
+    // Compute all the break points and sort them
+    std::set<field_type> bPValues;
+    field_type infinity = std::numeric_limits<field_type>::max();
+
+    VectorType bPCoords(rhs_->size());
+    bPCoords = 0;
+
+    const VectorType& oldIt = *this->x_;
+
+    for (size_t i=0; i<obstacles_->size(); i++)
+        for (int j=0; j<blocksize; j++) {
+
+            if (this->ignore()[i][j]) {
+                projGrad[i][j]=0;
+                continue;
+            }
+
+            bPCoords[i][j] = infinity;
+
+            if (projGrad[i][j]>0 && (*obstacles_)[i].upper(j)<infinity)
+                bPCoords[i][j] = ((*obstacles_)[i].upper(j)-oldIt[i][j])/projGrad[i][j];
+
+            if (projGrad[i][j]<0 && (*obstacles_)[i].lower(j)>-infinity)
+                bPCoords[i][j] = ((*obstacles_)[i].lower(j)-oldIt[i][j])/projGrad[i][j];
+
+            if (bPCoords[i][j] <1e-14)
+                projGrad[i][j] = 0;
+            else if (bPCoords[i][j] < infinity)
+                bPValues.insert(bPCoords[i][j]);
+        }
+
+    // Walk along the piecewise linear path until we find a minimizer
+    VectorType cauchyPoint(rhs_->size());
+    cauchyPoint = *this->x_;
+
+    field_type oldBreakPoint(0);
+    for ( auto bP : bPValues) {
+
+        // if there are multiple breakpoints at the same position then skip them
+        if (bP-oldBreakPoint<1e-14)
+            continue;
+
+        // Compute minimal value on the line segment
+        VectorType dummy(projGrad.size());
+        mat_->mv(projGrad,dummy);
+
+        field_type linTerm = -((*rhs_)*projGrad) + (dummy*cauchyPoint);
+        field_type localMin = -linTerm/(dummy*projGrad);
+
+        if (dummy*projGrad <0)
+            DUNE_THROW(Dune::Exception,"Not positive definite!\n");
+
+        // Check if we really have a local minimum in this segment
+        if (localMin < 0) {
+            //field_type min = computeEnergy(*mat_,cauchyPoint,*rhs_);
+            //std::cout<<"Generalized Cauchy Point "<<min<<" at previous breakpoint "<<oldBreakPoint<<std::endl;
+
+            *this->x_ = cauchyPoint;
+
+            return;
+
+        } else if (localMin<=(bP-oldBreakPoint)) {
+
+            cauchyPoint.axpy(localMin,projGrad);
+            //field_type min = computeEnergy(*mat_,cauchyPoint,*rhs_);
+            //std::cout<<"Generalized Cauchy Point found  "<<min<<" at local min "<<localMin
+            //         <<" BP "<<bP+localMin<<std::endl;
+
+            *this->x_ = cauchyPoint;
+
+            return;
+        }
+
+        // Move to next interval
+        cauchyPoint.axpy(bP-oldBreakPoint,projGrad);
+        oldBreakPoint = bP;
+
+        // Truncate projected gradient
+        for (size_t i=0; i<bPCoords.size(); i++)
+            for (int j=0; j<blocksize; j++)
+                if (std::fabs(bPCoords[i][j]-bP)<1e-14)
+                    projGrad[i][j] = 0;
+    }
+    *this->x_ += cauchyPoint;
+}

dune/solvers/iterationsteps/projectedgradientstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_PROJECTED_GRADIENT_STEP_HH
+#define DUNE_SOLVERS_PROJECTED_GRADIENT_STEP_HH
+
+/** \file
+    \brief Computes the exact generalized Cauchy point for a possibly nonconvex,
+           box-constraints quadratic problem.
+*/
+
+#include <memory>
+
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+
+/*
+   \brief Computes the exact generalized Cauchy point for a possibly nonconvex,
+   box-constraints quadratic problem, which is the first local minimizer along the
+   projected gradient path.
+ */
+template<class MatrixType, class VectorType>
+class ProjectedGradientStep : public IterationStep<VectorType,Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+
+    typedef typename VectorType::block_type VectorBlock;
+
+    enum {blocksize = VectorBlock::dimension};
+
+    typedef typename VectorType::field_type field_type;
+    typedef IterationStep<VectorType,Dune::BitSetVector<blocksize> > Base;
+
+public:
+
+    //! Default constructor.  Doesn't init anything
+    ProjectedGradientStep() {}
+
+    //! Constructor that gets the quadratic energy functional x^T mat x - rhs^T*x, and the initial iterate
+    ProjectedGradientStep(const MatrixType& mat, VectorType& x, const VectorType& rhs)
+        : Base(x), rhs_(&rhs)
+    {
+        mat_     = Dune::stackobject_to_shared_ptr(mat);
+    }
+
+    /** \brief Compute the generalized Cauchy point along a given direction. This is sometimes needed in
+     *         nonlinear programming where this direction then corresponds to the gradient of the nonlinear function
+     *         and not the gradient of the quadratic approximation. This can also be used as a criticality measure.
+     **/
+    void computeGeneralizedCP(const VectorType& dir);
+
+    //! Perform one iteration
+    virtual void iterate() {
+        VectorType negativeGradient=*rhs_;
+        mat_->mmv(*this->x_,negativeGradient);
+        computeGeneralizedCP(negativeGradient);
+    }
+
+    //! The obstacles
+    const std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
+private:
+    //! The container for the right hand side
+    const VectorType* rhs_;
+
+    //! The linear operator
+    std::shared_ptr<const MatrixType> mat_;
+};
+
+#include "projectedgradientstep.cc"
+
+#endif

dune/solvers/iterationsteps/projectedlinegsstep.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
 #include <dune/istl/scaledidmatrix.hh>
 
 #include <cmath>
@@ -19,38 +22,38 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
     // /////////////////////////////
     x = 0;
 
-    
+
     // //////////////////////////////////////
     //   The TNNMG loop
     // //////////////////////////////////////
 
     // stupid:
-    Dune::BitSetVector<1> hasObstacle(rhs.size(), true);
+    Dune::BitSetVector<BlockSize> hasObstacle(rhs.size(), true);
     // The nodes to ignore are already written into the matrix/rhs
     Dune::BitSetVector<1> ignoreNodes(rhs.size(), false);
 
     for (int iteration=0; iteration<10; iteration++) {
 
         VectorType oldX = x;
-        
+
         // ///////////////////////////
         //   Nonlinear presmoothing
         // ///////////////////////////
         VectorType rhsCopy = rhs;  // need a non-const version of rhs
-        
+
         ProjectedBlockGSStep<LocalMatrixType, VectorType> nonlinearSmootherStep(matrix, x, rhsCopy);
 
-        nonlinearSmootherStep.ignoreNodes_ = &ignoreNodes;
+        nonlinearSmootherStep.setIgnore(ignoreNodes);
         nonlinearSmootherStep.hasObstacle_ = &hasObstacle;
         nonlinearSmootherStep.obstacles_   = &localObstacles;
 
         for (int i=0; i<3; i++)
             nonlinearSmootherStep.iterate();
-        
+
         // Compute residual
         VectorType residual = rhs;
         matrix.mmv(x,residual);
-        
+
         // ///////////////////////////
         //   Truncation
         // ///////////////////////////
@@ -67,10 +70,10 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
                 // left off-diagonal:
                 if (i>0)
                     truncatedMatrix[i][i-1] = 0;
-                
+
                 // diagonal
                 truncatedMatrix[i][i] = Dune::ScaledIdentityMatrix<double,BlockSize>(1);
-                
+
                 // right off-diagonal:
                 if (i<truncatedMatrix.N()-1)
                     truncatedMatrix[i][i+1] = 0;
@@ -81,11 +84,11 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
             }
 
         }
-        
+
         // ///////////////////////////
         //   Linear correction
         // ///////////////////////////
-        
+
         VectorType linearCorrection(residual.size());
         truncatedMatrix.solve(linearCorrection, residual);   // the direct linear solution algorithm
 
@@ -93,56 +96,47 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
         //   Line search in the direction of the projected coarse
         //   grid correction to ensure monotonicity.
         // //////////////////////////////////////////////////////////
-        
+
         // L2-projection of the correction onto the defect obstacle
         for (size_t i=0; i<linearCorrection.size(); i++) {
-            
+
             for (int j=0; j<BlockSize; j++) {
-                
+
                 linearCorrection[i][j] = std::max(linearCorrection[i][j], localObstacles[i].lower(j) - x[i][j]);
                 linearCorrection[i][j] = std::min(linearCorrection[i][j], localObstacles[i].upper(j) - x[i][j]);
-                
+
             }
-            
+
         }
-        
+
         // Construct obstacles in the direction of the projected correction
         BoxConstraint<field_type,1> lineSearchObs;
         for (size_t i=0; i<linearCorrection.size(); i++) {
-            
+
             for (int j=0; j<BlockSize; j++) {
-                
+                // fmin and fmax ignore NAN values
                 if (linearCorrection[i][j] > 0) {
-                    
-                    // This division can cause nan on some platforms...
-                    if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]) )
-                        lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
-                                                          (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
-                    
-                    if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]) )
-                        lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0), 
-                                                          (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
+                    lineSearchObs.lower(0) = std::fmax(lineSearchObs.lower(0),
+                                                       (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
+                    lineSearchObs.upper(0) = std::fmin(lineSearchObs.upper(0),
+                                                       (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
                 }
-                
                 if (linearCorrection[i][j] < 0) {
-                    if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]) )
-                        lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0), 
-                                                          (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
-                    if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]) )
-                        lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0), 
-                                                          (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
+                    lineSearchObs.lower(0) = std::fmax(lineSearchObs.lower(0),
+                                                       (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
+                    lineSearchObs.upper(0) = std::fmin(lineSearchObs.upper(0),
+                                                       (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
                 }
-                
             }
-            
+
         }
 
         // abort when the linear correction has a small norm
         field_type correctionNormSquared = EnergyNorm<LocalMatrixType,VectorType>::normSquared(linearCorrection,truncatedMatrix);
-        
+
         // Line search
         field_type alpha = (residual*linearCorrection) / correctionNormSquared;
-        alpha = std::max(std::min(alpha, lineSearchObs.upper(0)), lineSearchObs.lower(0));
+        alpha = lineSearchObs[0].projectIn(alpha);
 
         if (isnan(alpha))
             alpha = 1;
@@ -184,13 +178,13 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
 
         const int current_block_size = this->blockStructure_[b_num].size();
 
-        //! compute and save the residuals for the curent block: 
+        //! compute and save the residuals for the current block:
         // determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)]
         // this means that we determine the residuals for the current block
         VectorType permuted_r_i(current_block_size);
         for (int k=0; k<current_block_size; k++)
           {
-           if ( (*this->ignoreNodes_)[this->blockStructure_[b_num][k]][0])
+           if ( this->ignore()[this->blockStructure_[b_num][k]][0])
                permuted_r_i[k] = 0.0;
            else
                this->residual( this->blockStructure_[b_num][k], permuted_r_i[k]); // get r[p(i+k)]
@@ -203,22 +197,22 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
         // (Later: x_now[p(i+k)] = x_now[p(i+k)] + v[p(i+k)] )
 
         // //////////////////////////////////////////////////////////////////////////////////////////////////
-        // Copy the linear system for the current line/block into a tridiagonal matrix 
+        // Copy the linear system for the current line/block into a tridiagonal matrix
         // //////////////////////////////////////////////////////////////////////////////////////////////////
 
         Dune::BTDMatrix<typename MatrixType::block_type> tridiagonalMatrix(current_block_size);
 
         for (int j=0; j<current_block_size; j++) {
 
-            if ( (*this->ignoreNodes_)[this->blockStructure_[b_num][j]][0] ) {
+            if ( this->ignore()[this->blockStructure_[b_num][j]][0] ) {
 
                 // left off-diagonal:
                 if (j>0)
                     tridiagonalMatrix[j][j-1] = 0;
-                
+
                 // diagonal
                 tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<field_type,BlockSize>(1);
-                
+
                 // left off-diagonal:
                 if (j<current_block_size-1)
                     tridiagonalMatrix[j][j+1] = 0;
@@ -228,10 +222,10 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
                 // left off-diagonal:
                 if (j>0)
                     tridiagonalMatrix[j][j-1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j-1]];
-                
+
                 // diagonal
                 tridiagonalMatrix[j][j] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j]];
-                
+
                 // left off-diagonal:
                 if (j<current_block_size-1)
                     tridiagonalMatrix[j][j+1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j+1]];

dune/solvers/iterationsteps/projectedlinegsstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_PROJECTED_LINE_GAUSS_SEIDEL_STEP_HH
 #define DUNE_PROJECTED_LINE_GAUSS_SEIDEL_STEP_HH
 
@@ -13,34 +15,34 @@ template<class MatrixType,
          class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
 class ProjectedLineGSStep : public LineGSStep<MatrixType, VectorType, BitVectorType>
 {
-    
+
     typedef typename VectorType::block_type VectorBlock;
-    
+
     typedef typename VectorType::field_type field_type;
-    
+
     enum {BlockSize = VectorBlock::dimension};
-    
+
 public:
-    
+
     //! Default constructor.  Doesn't init anything
     ProjectedLineGSStep() {}
-    
+
     //! Constructor for usage of Permutation Manager
     ProjectedLineGSStep( const std::vector<std::vector<unsigned int> >& blockStructure)
         : LineGSStep<MatrixType, VectorType, BitVectorType>( blockStructure )
     {}
-    
+
     /** \brief Solve one tridiagonal system */
     void solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
                           VectorType& x,
                           const VectorType& rhs,
                           const std::vector<BoxConstraint<field_type,BlockSize> >& localObstacles) const;
-    
+
     //! Perform one iteration
     virtual void iterate();
-    
+
     const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
-    
+
 };
 
 

dune/solvers/iterationsteps/richardsonstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_RICHARDSON_STEP_HH
 #define DUNE_RICHARDSON_STEP_HH
 
@@ -27,22 +29,12 @@ public:
     //! Perform one iteration
     virtual void iterate();
 
-    /** \brief Retrieve the solution */
-    virtual VectorType getSol();
-    
     const Preconditioner<VectorType>* preconditioner_;
 
     double damping_;
     
 };
 
-template<class VectorType, class BitVectorType>
-inline
-VectorType RichardsonStep<VectorType, BitVectorType>::getSol()
-{
-    return *(this->x_);
-}
-
 template<class VectorType, class BitVectorType>
 inline
 void RichardsonStep<VectorType, BitVectorType>::iterate()

dune/solvers/iterationsteps/truncatedblockgsstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef TRUNCATED_BLOCK_GAUSS_SEIDEL_STEP_HH
 #define TRUNCATED_BLOCK_GAUSS_SEIDEL_STEP_HH
 
@@ -5,11 +7,11 @@
 
 #include <dune/common/bitsetvector.hh>
 
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/promote.hh>
+
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
 #include <dune/solvers/operators/sumoperator.hh>
-#include <dune/solvers/common/staticmatrixtools.hh>
-#include <dune/solvers/common/arithmetic.hh>
-
 
 namespace TruncatedBlockGSStepNS {
 template<int blocklevel> struct RecursiveGSStep;
@@ -32,20 +34,15 @@ public:
     {}
 
     //! Constructor with a linear problem
-    TruncatedBlockGSStep(MatrixType& mat, VectorType& x, VectorType& rhs, int innerLoops=1)
+    TruncatedBlockGSStep(const MatrixType& mat, VectorType& x, const VectorType& rhs, int innerLoops=1)
         : LinearIterationStep<MatrixType,VectorType>(mat, x, rhs),
         innerLoops_(innerLoops)
     {}
 
-    virtual VectorType getSol()
-    {
-        return *(this->x_);
-    }
-
     //! Perform one iteration
     virtual void iterate()
     {
-        TruncatedBlockGSStepNS::RecursiveGSStep<MatrixType::blocklevel>::apply(*this->mat_, *this->rhs_, *this->ignoreNodes_, *this->x_, innerLoops_);
+        TruncatedBlockGSStepNS::RecursiveGSStep<Dune::blockLevel<MatrixType>()>::apply(*this->mat_, *this->rhs_, this->ignore(), *this->x_, innerLoops_);
     }
 
 protected:
@@ -99,7 +96,7 @@ template<>
 struct RecursiveGSStep<1>
 {
     template<class MType, class VType, class BVType>
-    static void apply(const MType& mat, const VType& rhs, const BVType& ignore, VType& x, int innerLoops)
+    static void apply(const MType& mat, const VType& rhs, const BVType& ignore, VType& x, [[maybe_unused]] int innerLoops)
     {
         typedef typename MType::block_type MBlock;
 
@@ -143,7 +140,7 @@ class TruncatedBlockGSStep<SumOperator<SparseMatrixType, LowRankMatrixType>, Vec
 
     typedef typename SparseMatrixType::block_type SpBlock;
     typedef typename LowRankMatrixType::block_type LRBlock;
-    typedef typename StaticMatrix::Promote<SpBlock, LRBlock>::Type MBlock;
+    typedef typename Dune::MatrixVector::Promote<SpBlock, LRBlock>::Type MBlock;
     typedef typename VectorType::block_type VBlock;
 
 public:
@@ -159,11 +156,6 @@ public:
         localSolver(localSolver)
     {}
 
-    virtual VectorType getSol()
-    {
-        return *(this->x_);
-    }
-
     //! Perform one iteration
     virtual void iterate()
     {
@@ -174,7 +166,7 @@ public:
         const LowRankMatrixType& lowrank_mat = mat.lowRankMatrix();
         const VectorType& rhs = *this->rhs_;
         VectorType& x   = *this->x_;
-        const BitVectorType& ignore = *this->ignoreNodes_;
+        const BitVectorType& ignore = this->ignore();
 
         // compute m*x once and only add/subtract single summands in the iterations later
         VectorType mx(lowrank_mat.lowRankFactor().N());
@@ -194,9 +186,7 @@ public:
             {
                 size_t col = it.index();
                 if (col == row)
-                {
-                    StaticMatrix::axpy(Aii, 1.0, *it);
-                }
+                    Dune::MatrixVector::addProduct(Aii, 1.0, *it);
                 else
                     it->mmv(x[col],r);
             }
@@ -208,7 +198,7 @@ public:
                 const LRBlock& mki(lowrank_mat.lowRankFactor()[k][row]);
                 // Aii += m^tm[i][i]
                 // this should actually be \sum(lr_block^t*lr_block) if the blocks are nonsymmetric
-                Arithmetic::addProduct(Aii,mki,mki);
+                Dune::MatrixVector::addProduct(Aii,mki,mki);
 
                 // r -= m^tm*x, where x[row] were zero
                 mki.mmv(mx[k], r);
@@ -293,19 +283,30 @@ public:
         {
             for(typename MBlock::size_type i=0; i<A.N(); ++i)
             {
+                assert(A.exists(i,i)); // No implementation for matrices w/o explicit diagonal elements is available here.
+
+                bool rowIsZero = true;
                 if (ignore[i])
                 {
                     A[i] = 0.0;
                 }
+                else
+                {
+                    typename MBlock::row_type::Iterator inner_it = A[i].begin();
+                    typename MBlock::row_type::Iterator inner_end = A[i].end();
+                    for(; inner_it!=inner_end; ++inner_it)
+                        if (*inner_it != 0.0)
+                        {
+                            rowIsZero = false;
+                            break;
+                        }
+                }
 
-                typename MBlock::row_type::Iterator inner_it = A[i].begin();
-                typename MBlock::row_type::Iterator inner_end = A[i].end();
-                for(; inner_it!=inner_end; ++inner_it)
-                    if (inner_it.index()==i and *inner_it==0.0)
-                    {
-                        *inner_it = 1.0;
-                        b[i] = x[i];
-                    }
+                if (rowIsZero)
+                {
+                    A[i][i] = 1.0;
+                    b[i] = x[i];
+                }
             }
 
             A.solve(x,b);

dune/solvers/iterationsteps/truncatedsaddlepointgsstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef TRUNCATED_SADDLE_POINT_GAUSS_SEIDEL_STEP_HH
 #define TRUNCATED_SADDLE_POINT_GAUSS_SEIDEL_STEP_HH
 
@@ -24,17 +26,11 @@ public:
         : LinearIterationStep<MatrixType,VectorType>(mat, x, rhs)
     {}
 
-    virtual VectorType getSol()
-    {
-        return *(this->x_);
-    }
-
     //! Perform one iteration
     virtual void iterate()
     {
         const MatrixType& mat = *this->mat_;
         const VectorType& rhs = *this->rhs_;
-        const BitVectorType& ignore = *this->ignoreNodes_;
         VectorType& x = *this->x_;
 
         typedef typename VectorType::block_type VectorBlock;

dune/solvers/iterationsteps/trustregiongsstep.cc

-template<class MatrixType, class VectorType>
-inline
-VectorType TrustRegionGSStep<MatrixType, VectorType>::getSol()
-{
-    return *(this->x_);
-}
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 
 template<class MatrixType, class VectorType>
 inline
-typename VectorType::field_type TrustRegionGSStep<MatrixType, VectorType>::
-residual(unsigned int index, int blockIndex) const
+void TrustRegionGSStep<MatrixType, VectorType>::iterate()
 {
     const MatrixType& mat = *this->mat_;
+    VectorType& x = *this->x_;
+    const std::vector<BoxConstraint<Field,BlockSize> >& obstacles = *this->obstacles_;
 
-    typedef typename MatrixType::row_type RowType;
-    typedef typename RowType::ConstIterator ColumnIterator;
-
-    /* The following loop computes
-     * \f[ sum_i = \sum_{j \ne i} A_{ij}x_j \f]
-     */
-    field_type r = (*this->rhs_)[index][blockIndex];
-    ColumnIterator cIt    = mat[index].template begin();
-    ColumnIterator cEndIt = mat[index].template end();
-
-    for (; cIt!=cEndIt; ++cIt) {
-        // r_i -= A_ij x_j
-        // r -= (*cIt)[blockIndex] * (*this->x_)[cIt.index()];
-        for (int i=0; i<blocksize; i++) {
-
-            // Omit diagonal element
-            if (cIt.index()==index && i==blockIndex)
-                continue;
-
-            r -= (*cIt)[blockIndex][i] * (*this->x_)[cIt.index()][i];
-
-        }
-
-    }
-
-    //std::cout << "+++ index " << index << "  rhs: " << (*this->rhs_)[index] << std::endl;
-    return r;
+    for (size_t i=0; i<x.size(); i++) {
 
-}
+        // Dirichlet nodes in this block
+        std::bitset<BlockSize> ignoreNodes(0);
+        if (this->hasIgnore())
+            ignoreNodes = this->ignore()[i];
 
-template<class MatrixType, class VectorType>
-inline
-void TrustRegionGSStep<MatrixType, VectorType>::iterate()
-{
-    const MatrixType& mat = *this->mat_;
-    const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_;
+        if (ignoreNodes.all())
+            continue;
 
-    // Debug: check for energy decrease
-    //std::cout << "----------------------------------\n";
-    //double oldEnergy = computeEnergy(mat, *this->x_, *this->rhs_);
+        VectorBlock blockResidual;
+        this->residual(i, blockResidual);
 
-    for (size_t i=0; i<this->x_->size(); i++) {
+        mat[i][i].umv(x[i], blockResidual);
 
         // scalar gauss-seidel
-        for (int j=0; j<blocksize; j++) {
+        for (size_t j=0; j<BlockSize; j++) {
 
-            if (this->ignoreNodes_ and (*this->ignoreNodes_)[i][j])
+            if (ignoreNodes[j])
                 continue;
 
             // Compute residual
-            field_type r = residual(i, j);
-            const field_type& diag = mat[i][i][j][j];
+            Field r = blockResidual[j] - mat[i][i][j] * x[i];
+            r += mat[i][i][j][j] * x[i][j];
 
-            // Find local minimum
-            field_type& x = (*this->x_)[i][j];
+            const Field& diag = mat[i][i][j][j];
 
+            // Find local minimum
             if (diag > 0) {
-                //printf("%d %d konvex\n", i, j);
-                //double oldEnergy = computeEnergy(mat, *this->x_, *this->rhs_);
                 // 1d problem is convex
-#ifndef NDEBUG
-                field_type oldX = x;
-#endif
-                x = r / diag;
-#ifndef NDEBUG
-                field_type energyX = 0.5*diag*x*x - r*x;
-                field_type energyOldX = 0.5*diag*oldX*oldX - r*oldX;
-                assert(energyX <= energyOldX + 1e-6);
-#endif
-                if ( x < obstacles[i].lower(j)) {
-                    // assert that projection increases functional value
-//                     double trueMin = 0.5*diag*x*x - r*x;
-//                     double proMin  = 0.5*diag*obstacles[i].lower(j)*obstacles[i].lower(j) - r*obstacles[i].lower(j);
-                    //assert(trueMin < proMin);
-                    x = obstacles[i].lower(j);
-                }
-                else if ( x > obstacles[i].upper(j)) {
-                    // assert that projection increases functional value
-//                     double trueMin = 0.5*diag*x*x - r*x;
-//                     double proMin  = 0.5*diag*obstacles[i].upper(j)*obstacles[i].upper(j) - r*obstacles[i].upper(j);
-                    //assert(trueMin < proMin);
-                    x = obstacles[i].upper(j);
-                }
-
-                // Debug: check for energy decrease
-                //double energy = computeEnergy(mat, *this->x_, *this->rhs_);
-                //assert(oldEnergy+1e-6>=energy);
+                x[i][j] = r / diag;
+
+                x[i][j] = obstacles[i][j].projectIn(x[i][j]);
 
             } else  {
-                //printf("%d %d konkav\n", i, j);
-                // 1d problem is concave or linear.
-                // Minimum is attained at one of the boundaries
-                field_type lBound = obstacles[i].lower(j);
-                field_type uBound = obstacles[i].upper(j);
 
-                field_type lValue = 0.5*diag*lBound*lBound - r*lBound;
-                field_type uValue = 0.5*diag*uBound*uBound - r*uBound;
+                // If the corresponding dof was truncated, then compute no correction
+                if (diag>-1e-13 && std::fabs(r)<1e-13)
+                    continue;
 
-#if 0
-                double max = 0.5*diag*(r/diag)*(r/diag) - r*r/diag;
-                assert(max >= lValue);
-                assert(max >= uValue);
-#endif
+                // 1d problem is concave or linear.
+                // Minimum is attained at one of the boundaries
+                Field lBound = obstacles[i].lower(j);
+                Field uBound = obstacles[i].upper(j);
 
-                x = (lValue < uValue) ? lBound : uBound;
+                Field lValue = 0.5*diag*lBound*lBound - r*lBound;
+                Field uValue = 0.5*diag*uBound*uBound - r*uBound;
 
-                // Debug: check for energy decrease
-                //double energy = computeEnergy(mat, *this->x_, *this->rhs_);
-                //assert(oldEnergy+1e-6>=energy);
+                x[i][j] = (lValue < uValue) ? lBound : uBound;
 
             }
 
@@ -128,7 +65,4 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
 
     }
 
-    // Debug: check for energy decrease
-    //double energy = computeEnergy(mat, *this->x_, *this->rhs_);
-    //assert(oldEnergy+1e-6>=energy);
 }

dune/solvers/iterationsteps/trustregiongsstep.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_TRUST_REGION_BLOCK_GAUSS_SEIDEL_STEP_HH
 #define DUNE_TRUST_REGION_BLOCK_GAUSS_SEIDEL_STEP_HH
 
@@ -15,15 +17,12 @@
     template<class MatrixType, class VectorType>
     class TrustRegionGSStep : public ProjectedBlockGSStep<MatrixType, VectorType>
     {
+        using VectorBlock = typename VectorType::block_type;
+        using Field = typename VectorType::field_type;
 
-        typedef typename VectorType::block_type VectorBlock;
-
-        enum {blocksize = VectorBlock::dimension};
-
-        typedef typename VectorType::field_type field_type;
+        enum {BlockSize = VectorBlock::dimension};
 
     public:
-
         //! Default constructor.  Doesn't init anything
         TrustRegionGSStep() {}
 
@@ -34,11 +33,6 @@
 
         //! Perform one iteration
         virtual void iterate();
-
-        virtual VectorType getSol();
-
-        field_type residual(unsigned int index, int blockIndex) const;
-
     };
 
 #include "trustregiongsstep.cc"

dune/solvers/norms/Makefile.am

-SUBDIRS =
-
-normsdir = $(includedir)/dune/solvers/norms
-norms_HEADERS = \
-	blocknorm.hh \
-	diagnorm.hh \
-	energynorm.hh \
-	fullnorm.hh \
-	h1seminorm.hh \
-	norm.hh \
-	pnorm.hh \
-	reorderedblocknorm.hh \
-	sumnorm.hh \
-	twonorm.hh
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules
\ No newline at end of file

dune/solvers/norms/blocknorm.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef BLOCK_NORM_HH
 #define BLOCK_NORM_HH
 
+#include <memory>
 #include <vector>
 #include <cmath>
 
-#include <dune/common/shared_ptr.hh>
-
 #include "norm.hh"
 
 //! A norm for blocks of vectors
-template <class VectorType>
-class BlockNorm: public Norm<VectorType>
+template <class V>
+class BlockNorm: public Norm<V>
 {
     public:
+        using VectorType = V;
+        using Base = Norm<V>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
+
 //        typedef typename std::vector<const Norm<typename VectorType::block_type>* > NormPointerVector;
-        typedef std::vector< Dune::shared_ptr<const Norm<typename VectorType::block_type> > > NormPointerVector;
+        typedef std::vector< std::shared_ptr<const Norm<typename VectorType::block_type> > > NormPointerVector;
 
         BlockNorm(const NormPointerVector& norms) :
             norms_(norms)
         {}
 
         //! Compute the norm of the given vector
-        double operator()(const VectorType &v) const
+        field_type operator()(const VectorType &v) const override
         {
-            double r = 0.0;
+            field_type r = 0.0;
             for (int i=0; i<norms_.size(); ++i)
             {
-                double ri = (*norms_[i])(v[i]);
+                auto ri = (*norms_[i])(v[i]);
                 r += ri*ri;
             }
             return std::sqrt(r);
         }
 
         //! Compute the norm of the difference of two vectors
-        double diff(const VectorType &v1, const VectorType &v2) const
+        field_type diff(const VectorType &v1, const VectorType &v2) const override
         {
-            double r = 0.0;
+            field_type r = 0.0;
             for (int i=0; i<norms_.size(); ++i)
             {
-                double ri = (*norms_[i]).diff(v1[i], v2[i]);
+                auto ri = (*norms_[i]).diff(v1[i], v2[i]);
                 r += ri*ri;
             }
             return std::sqrt(r);

dune/solvers/norms/diagnorm.hh

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DIAGNORM_HH
 #define DIAGNORM_HH
 
@@ -21,7 +21,10 @@ class DiagNorm:
 {
     public:
         typedef V VectorType;
-        typedef typename VectorType::field_type field_type;
+        using Base = Norm<V>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
     private:
         typedef typename VectorType::size_type SizeType;
 
@@ -32,13 +35,13 @@ class DiagNorm:
         {}
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType &f) const
+        field_type operator()(const VectorType &f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         //! Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             field_type r = 0.0;
 
@@ -53,7 +56,7 @@ class DiagNorm:
         }
 
         //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType &f1, const VectorType &f2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const override
         {
             field_type r = 0.0;
 
@@ -79,8 +82,11 @@ class DiagNorm<Dune::BlockVector<Dune::FieldVector <double,1> >, Dune::BlockVect
     public Norm<Dune::BlockVector<Dune::FieldVector <double,1> > >
 {
     public:
-        typedef double field_type;
-        typedef Dune::BlockVector<Dune::FieldVector <field_type,1> > VectorType;
+        typedef Dune::BlockVector<Dune::FieldVector <double,1> > VectorType;
+        using Base = Norm<VectorType>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
     private:
         typedef VectorType::size_type SizeType;
 
@@ -91,13 +97,13 @@ class DiagNorm<Dune::BlockVector<Dune::FieldVector <double,1> >, Dune::BlockVect
         {}
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType &f) const
+        field_type operator()(const VectorType &f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         //! Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             field_type r = 0.0;
 
@@ -108,7 +114,7 @@ class DiagNorm<Dune::BlockVector<Dune::FieldVector <double,1> >, Dune::BlockVect
         }
 
         //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType &f1, const VectorType &f2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const override
         {
             field_type r = 0.0;
 

dune/solvers/norms/energynorm.hh

-#ifndef ENERGY_NORM_HH
-#define ENERGY_NORM_HH
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_NORMS_ENERGY_NORM_HH
+#define DUNE_SOLVERS_NORMS_ENERGY_NORM_HH
 
 #include <cmath>
+#include <functional>
 
-#include "norm.hh"
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axy.hh>
+
+#include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
 
+namespace Dune {
+namespace Solvers {
+
     /** \brief Vector norm induced by linear operator
      *
      *  \f$\Vert u \Vert_A = (u, Au)^{1/2}\f$
@@ -23,97 +30,101 @@
      *
      *  \todo Elaborate documentation.
      */
-    template<class OperatorType, class V>
+    template<class MatrixType, class V>
     class EnergyNorm : public Norm<V>
     {
     public:
         typedef V VectorType;
+        using Base = Norm<V>;
 
         /** \brief The type used for the result */
-        typedef typename VectorType::field_type field_type;
+        using typename Base::field_type;
 
-        EnergyNorm(const double tol=1e-10 ) : iterationStep_(NULL), matrix_(NULL), tol_(tol) {}
+        EnergyNorm(const field_type tol=1e-10 ) : tol_(tol) {}
 
-        EnergyNorm(LinearIterationStep<OperatorType, VectorType>& it, const double tol=1e-10)
-            : iterationStep_(&it), matrix_(NULL), tol_(tol)
-        {}
+        template<class BV>
+        EnergyNorm(LinearIterationStep<MatrixType, VectorType, BV>& it, const field_type tol=1e-10)
+            : EnergyNorm(tol)
+        {
+            setIterationStep(&it);
+        }
 
-        EnergyNorm(const OperatorType& matrix, const double tol=1e-10)
-            : iterationStep_(NULL), matrix_(&matrix), tol_(tol)
-        {}
+        EnergyNorm(const MatrixType& matrix, const field_type tol=1e-10)
+            : EnergyNorm(tol)
+        {
+            setMatrix(&matrix);
+        }
 
-        void setMatrix(const OperatorType* matrix) {
-            matrix_ = matrix;
+        //! \brief sets the energy norm matrix
+        void setMatrix(const MatrixType* matrix) {
+            matrixProvider_ = [=]() -> const MatrixType& { return *matrix; };
         }
 
-        void setIterationStep(LinearIterationStep<OperatorType, VectorType>* iterationStep) {
-            iterationStep_ = iterationStep;
+        //! \brief get the energy norm matrix
+        const MatrixType& getMatrix() const {
+            return matrixProvider_();
         }
 
-        //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType& f1, const VectorType& f2) const {
-            if (iterationStep_ == NULL && matrix_ == NULL)
-                DUNE_THROW(Dune::Exception, "You have supplied neither a matrix nor an IterationStep to the EnergyNorm!");
+        //! \brief sets to use the current problem matrix of the linear iteration step
+        template<class BV>
+        void setIterationStep(LinearIterationStep<MatrixType, VectorType, BV>* step) {
+            matrixProvider_ = [=]() -> const MatrixType& { return *step->getMatrix(); };
+        }
 
+        //! Compute the norm of the difference of two vectors
+        field_type diff(const VectorType& f1, const VectorType& f2) const override {
             VectorType tmp_f = f1;
             tmp_f -= f2;
             return (*this)(tmp_f);
         }
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType& f) const
+        field_type operator()(const VectorType& f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         // \brief Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
-            if (iterationStep_ == NULL && matrix_ == NULL)
+            if (not matrixProvider_)
                 DUNE_THROW(Dune::Exception, "You have supplied neither a matrix nor an IterationStep to the EnergyNorm!");
-
-            const OperatorType& A = (iterationStep_)
-                ? *(iterationStep_->getMatrix())
-                : *matrix_;
-
-            const field_type ret = Arithmetic::Axy(A, f, f);
-
-            if (ret < 0)
-            {
-                if (ret < -tol_)
-                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << ret);
-                return 0.0;
-            }
-
-            return ret;
+            const field_type ret = Dune::MatrixVector::Axy(matrixProvider_(), f, f);
+            return checkedValue(ret, tol_);
         }
 
         // \brief Compute the squared norm for a given vector and matrix
-        DUNE_DEPRECATED static field_type normSquared(const VectorType& u,
-                                                      const OperatorType& A,
-                                                      const double tol=1e-10)
+        [[deprecated]] static field_type normSquared(const VectorType& u,
+                                                      const MatrixType& A,
+                                                      const field_type tol=1e-10)
         {
-            const field_type ret = Arithmetic::Axy(A, u, u);
-
-            if (ret < 0)
-            {
-                if (ret < -tol)
-                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << ret);
-                return 0.0;
-            }
-
-            return ret;
-
-        };
+            const field_type ret = Dune::MatrixVector::Axy(A, u, u);
+            return checkedValue(ret, tol);
+        }
 
     protected:
 
-        LinearIterationStep<OperatorType, VectorType>* iterationStep_;
+        //! yields access to energy matrix
+        std::function<const MatrixType&()> matrixProvider_;
 
-        const OperatorType* matrix_;
-
-        const double tol_;
+        const field_type tol_;
 
+        //! \brief throw an exception if value is below tolerance, project to R^+ otherwise.
+        static field_type checkedValue(field_type value, field_type tolerance)
+        {
+            if (value < 0) {
+                if (value < -tolerance)
+                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << value);
+                else return 0.0;
+            }
+            return value;
+        }
     };
 
+} /* namespace Solvers */
+} /* namespace Dune */
+
+// For backward compatibility: will be removed eventually
+using Dune::Solvers::EnergyNorm;
+
 #endif

dune/solvers/norms/fullnorm.hh

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef FULLNORM_HH
 #define FULLNORM_HH
 
@@ -23,7 +23,10 @@ class FullNorm: public Norm<V>
 {
     public:
         typedef V VectorType;
-        typedef typename VectorType::field_type field_type;
+        using Base = Norm<V>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
 
         FullNorm(const field_type alpha, const LowRankFactor &lowRankFactor) :
             lowRankFactor_(lowRankFactor),
@@ -31,13 +34,13 @@ class FullNorm: public Norm<V>
         {}
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType &f) const
+        field_type operator()(const VectorType &f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         //! Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             VectorType r(lowRankFactor_.N());
             r = 0.0;
@@ -48,7 +51,7 @@ class FullNorm: public Norm<V>
         }
 
         //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType &f1, const VectorType &f2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const override
         {
             VectorType r(lowRankFactor_.N());
             r = 0.0;
@@ -72,8 +75,11 @@ class FullNorm<Dune::BlockVector<Dune::FieldVector<double,1> >, Dune::BlockVecto
     public Norm<Dune::BlockVector<Dune::FieldVector<double,1> > >
 {
     public:
-        typedef double field_type;
-        typedef Dune::BlockVector<Dune::FieldVector<field_type,1> > VectorType;
+        typedef Dune::BlockVector<Dune::FieldVector<double,1> > VectorType;
+        using Base = Norm<VectorType>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
     private:
         typedef VectorType::size_type SizeType;
     public:
@@ -84,13 +90,13 @@ class FullNorm<Dune::BlockVector<Dune::FieldVector<double,1> >, Dune::BlockVecto
         {}
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType &f) const
+        field_type operator()(const VectorType &f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         //! Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             field_type r = 0.0;
 
@@ -101,7 +107,7 @@ class FullNorm<Dune::BlockVector<Dune::FieldVector<double,1> >, Dune::BlockVecto
         }
 
         //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType &f1, const VectorType &f2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const override
         {
             field_type r = 0.0;
 

dune/solvers/norms/h1seminorm.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_H1_SEMINORM_HH
 #define DUNE_H1_SEMINORM_HH
 
@@ -6,61 +8,79 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
 
+#include <dune/solvers/common/wrapownshare.hh>
 #include "norm.hh"
 
 //! Specialisation of the EnergyNorm class to identity blocks
-template<class VectorType>
-class H1SemiNorm : public Norm<VectorType>
+template<class V, class M = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >>
+class H1SemiNorm : public Norm<V>
 {
 public:
-    H1SemiNorm() : matrix_(NULL) {}
-    
-    H1SemiNorm(const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >& matrix) 
-        : matrix_(&matrix)
-    {}
+    typedef V VectorType;
+    using Base = Norm<V>;
+
+    using MatrixType = M;
+
+    /** \brief The type used for the result */
+    using typename Base::field_type;
+
+    H1SemiNorm() = default;
+
+    template<typename Matrix>
+    H1SemiNorm(Matrix&& matrix)
+        : matrix_(Dune::Solvers::wrap_own_share<const MatrixType>(std::forward<Matrix>(matrix)))
+    {
+      assert(matrix_->N() == matrix_->M());
+    }
 
     //! Compute the norm of the difference of two vectors
-    double diff(const VectorType& u1, const VectorType& u2) const 
-	{
-		double sum = 0;
-		
-		for (size_t i=0; i<matrix_->N(); i++) 
+    field_type diff(const VectorType& u1, const VectorType& u2) const override
+    {
+          assert(u1.size()==u2.size());
+          assert(u1.size()==matrix_->N());
+
+    field_type sum = 0;
+
+		for (size_t i=0; i<matrix_->N(); i++)
 		{
-			typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cIt    = (*matrix_)[i].begin();
-			typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cEndIt = (*matrix_)[i].end();
-			
-                        typename VectorType::block_type differ_i = u1[i] - u2[i];
+      auto cIt    = (*matrix_)[i].begin();
+      auto cEndIt  = (*matrix_)[i].end();
+
+      typename VectorType::block_type differ_i = u1[i] - u2[i];
 			for (; cIt!=cEndIt; ++cIt)
 				sum += differ_i * (u1[cIt.index()]-u2[cIt.index()]) * (*cIt);
         }
 
         return std::sqrt(sum);
     }
-    
+
     //! Compute the norm of the given vector
-    double operator()(const VectorType& u) const
+    field_type operator()(const VectorType& u) const override
     {
         if (matrix_ == NULL)
             DUNE_THROW(Dune::Exception, "You have not supplied neither a matrix nor an IterationStep to the EnergyNorm routine!");
+
+        assert(u.size()==matrix_->N());
+
         // we compute sqrt{uAu^t}.  We have implemented this by hand because the matrix is
         // always scalar but the vectors may not be
-        double sum = 0;
-        
+        field_type sum = 0;
+
         for (size_t i=0; i<matrix_->N(); i++) {
-            
-            typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cIt    = (*matrix_)[i].begin();
-            typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cEndIt = (*matrix_)[i].end();
-            
+
+            auto cIt    = (*matrix_)[i].begin();
+            auto cEndIt = (*matrix_)[i].end();
+
             for (; cIt!=cEndIt; ++cIt)
                 sum += u[i]*u[cIt.index()] * (*cIt);
-            
+
         }
-        
+
         return std::sqrt(sum);
     }
-    
-    const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* matrix_;
-    
+
+    std::shared_ptr<const MatrixType> matrix_;
+
 };
 
 #endif

dune/solvers/norms/norm.hh

-#ifndef NORM_HH
-#define NORM_HH
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_NORMS_NORM_HH
+#define DUNE_SOLVERS_NORMS_NORM_HH
+
+#include <dune/common/ftraits.hh>
+
+namespace Dune {
+namespace Solvers {
 
 //! Abstract base for classes computing norms of discrete functions
 template <class V>
@@ -8,10 +15,10 @@ class Norm {
     public:
         typedef V VectorType;
 
-        typedef typename VectorType::field_type field_type;
+        using field_type =  typename Dune::FieldTraits<VectorType>::field_type;
 
         /** \brief Destructor, doing nothing */
-        virtual ~Norm() {};
+        virtual ~Norm() {}
 
         //! Compute the norm of the given vector
         virtual field_type operator()(const VectorType& f) const = 0;
@@ -28,4 +35,10 @@ class Norm {
 
 };
 
+} /* namespace Solvers */
+} /* namespace Dune */
+
+// For backward compatibility: will be removed eventually
+using Dune::Solvers::Norm;
+
 #endif

dune/solvers/norms/pnorm.hh

-#ifndef __PNORM__HH__
-#define __PNORM__HH__
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_PNORM_HH
+#define DUNE_SOLVERS_PNORM_HH
 
 #include <cmath>
 
@@ -8,18 +10,23 @@
 
 #include "norm.hh"
 
-typedef Dune::BlockVector<Dune::FieldVector <double,1> > Vector;
-
-class PNorm: public Norm<Vector>
+template <class V>
+class PNorm: public Norm<V>
 {
 	public:
-		PNorm(int p=2, double alpha=1.0):
+  typedef V VectorType;
+  using Base = Norm<V>;
+
+  /** \brief The type used for the result */
+  using typename Base::field_type;
+
+  PNorm(int p=2, field_type alpha=1.0):
 			p(p),
 			alpha(alpha)
 		{}
 		
-		double operator()(const Vector &v) const
-		{
+    field_type operator()(const VectorType &v) const override
+    {
 			double r = 0.0;
 			
 			if (p<1)
@@ -50,17 +57,17 @@ class PNorm: public Norm<Vector>
 			}
 			
 			return alpha*r;
-		}
+    }
 
-		double diff(const Vector &v1, const Vector &v2) const
+    field_type diff(const VectorType &v1, const VectorType &v2) const override
 		{
-			double r = 0.0;
+      field_type r = 0.0;
 			
 			if (p<1)
 			{
 				for(int row = 0; row < v1.size(); ++row)
 				{
-					double z = std::abs(v1[row]-v2[row]);
+          field_type z = std::abs(v1[row]-v2[row]);
 					if (z>r)
 						r = z;
 				}
@@ -87,8 +94,8 @@ class PNorm: public Norm<Vector>
 		}
 
 	private:
-		double alpha;
-		double p;
+    field_type alpha;
+    int p;
 		
 };
 

dune/solvers/norms/reorderedblocknorm.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef REORDERED_BLOCK_NORM_HH
 #define REORDERED_BLOCK_NORM_HH
 
+#include <memory>
 #include <vector>
 #include <cmath>
 
-#include <dune/common/shared_ptr.hh>
 #include <dune/istl/bvector.hh>
 
 #include "../common/genericvectortools.hh"
 #include "norm.hh"
 
 //! A norm for blocks of interlaced vectors
-template <class VectorType, class ReorderedLocalVector>
-class ReorderedBlockNorm: public Norm<VectorType>
+template <class V, class ReorderedLocalVector>
+class ReorderedBlockNorm: public Norm<V>
 {
     public:
+        typedef V VectorType;
+        using Base = Norm<V>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
+
         typedef Dune::BlockVector<ReorderedLocalVector> ReorderedVector;
 
-        typedef std::vector< Dune::shared_ptr<const Norm<ReorderedLocalVector> > > NormPointerVector;
+        typedef std::vector< std::shared_ptr<const Norm<ReorderedLocalVector> > > NormPointerVector;
 
         ReorderedBlockNorm(const NormPointerVector& norms) :
             norms_(norms)
         {}
 
         //! Compute the norm of the given vector
-        double operator()(const VectorType &v) const
+        field_type operator()(const VectorType &v) const override
         {
             return std::sqrt(normSquared(v));
         }
 
         //! Compute the square of the norm of the given vector
-        double normSquared(const VectorType& v) const
+        field_type normSquared(const VectorType& v) const override
         {
             double r = 0.0;
 
@@ -46,7 +54,7 @@ class ReorderedBlockNorm: public Norm<VectorType>
         }
 
         //! Compute the norm of the difference of two vectors
-        double diff(const VectorType &v1, const VectorType &v2) const
+        field_type diff(const VectorType &v1, const VectorType &v2) const override
         {
             double r = 0.0;