base interface for linear solvers defined, started to implement ITL solvers

dune/amdis/linear_algebra/LinearSolverInterface.hpp

+#pragma once
+
+/** \file LinearSolverInterface.h */
+
+/**
+ * \defgroup Solver Solver module
+ * @{ <img src="solver.png"> @}
+ *
+ * \brief
+ * Contains all classes needed for solving linear and non linear equation
+ * systems.
+ */
+
+namespace AMDiS
+{
+  /// Non-templated base-class for Preconditioner types
+  struct PreconditionerInterface
+  {
+    virtual ~PreconditionerInterface() {}
+  };
+
+  /// Non-templates base-class for Runner / Worker types
+  struct RunnerInterface
+  {
+    // virtual destructor
+    virtual ~RunnerInterface() {}
+
+    virtual PreconditionerInterface* getLeftPrecon()
+    {
+      return NULL;
+    }
+
+    virtual PreconditionerInterface* getRightPrecon()
+    {
+      return NULL;
+    }
+  };
+
+  // A general solver class for istl solvers
+  // ---------------------------------------------------------------------------
+
+  template <class Matrix, class VectorX, class VectorB = VectorX>
+  class LinearSolverInterface
+  {    
+    /// The constructor reads needed parameters and sets solvers \p prefix.
+    /**
+     * Reads parameters for a solver with name 'NAME':
+     *   NAME->absolute tolerance   \ref aTol
+     *   NAME->relative tolerance   \ref rTol
+     *   NAME->max iteration        \ref maxTter
+     *   NAME->print cycle          \ref printCycle
+     *   NAME->info                 \ref info
+     *   NAME->break if tolerance not reached \ref breakTolNotReached
+    **/
+    explicit LinearSolverInterface(std::string prefix);
+
+    /// Destructor.
+    virtual ~LinearSolverInterface() {}
+    
+    
+    /// Public method to call in order to solve a linear system Ax = b.
+    /**
+     * The method redirects to the specific linear solver and prints statistics
+     * and error estimations at the end.
+     * 
+     * The parameters correspond to
+     *  \p A     A block-matrix that represents the system-matrix.
+     *  \p x     A block-vector for the unknown components.
+     *  \p b     A block-vector for the right-hand side of the linear system.
+     *  \p createMatrixData  If true, the matrix will be initialized and the 
+     *                       corresponding runner of the system receives the 
+     *                       matrix in the init() method.
+     *  \p storeMatrixData   If true, the exit() method of the runner will be 
+     *                       called.
+     **/
+    void solve(Matrix const& A, VectorX& x, VectorB const& b,
+               bool createMatrixData,
+               bool storeMatrixData)
+    {      
+      solveImpl(A, x, b, createMatrixData, storeMatrixData);
+    }
+
+    /** \name getting methods
+     * \{
+     */
+
+    /// Returns \ref aTol
+    double getAbsTolerance() const
+    {
+      return aTol;
+    }
+
+    /// Returns \ref rTol
+    double getRelTolerance() const
+    {
+      return rTol;
+    }
+
+    /// Returns \ref maxIter
+    int getMaxIterations() const
+    {
+      return maxIter;
+    }
+
+    /// Returns number of iterations in last run of an iterative solver
+    int getIterations() const
+    {
+      return iterations;
+    }
+
+    /// Returns error code in last run of an iterative solver
+    int getErrorCode() const
+    {
+      return error;
+    }
+
+    /// Returns info
+    int getInfo() const
+    {
+      return info;
+    }
+
+    /// Returns \ref printCycle
+    int getPrintCycle() const
+    {
+      return printCycle;
+    }
+
+    /// Returns \ref residual
+    double getResidual() const
+    {
+      return residual;
+    }
+    
+    /// Returns the initfile \ref prefix
+    std::string getPrefix() const
+    {
+      return prefix;
+    }
+
+    /** \} */
+
+
+    /** \name setting methods
+     * \{
+     */
+
+    /// Sets \ref aTol
+    void setAbsTolerance(double tol)
+    {
+      aTol = tol;
+    }
+
+    /// Sets \ref rTol
+    void setAbsTolerance(double tol)
+    {
+      rTol = tol;
+    }
+
+    /// Sets \ref maxIter
+    void setMaxIterations(int i)
+    {
+      maxIter = i;
+    }
+
+    /// Sets \ref info
+    void setInfo(int i)
+    {
+      info = i;
+    }
+
+    /// Sets \ref info
+    void setError(int e)
+    {
+      error = e;
+    }
+
+    /** \} */
+
+  private:
+    /// main methods that all solvers must implement
+    virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
+                           bool createMatrixData,
+                           bool storeMatrixData) = 0;
+    
+  private:
+    /// The initfile prefix to read parameters
+    std::string prefix;
+    
+    /// The name of the solver
+    std::string solverName = "cg";
+    
+    /// The name of the preconditioner (not used by every solver)
+    std::string preconName = "no";
+    
+  private:
+    /// The abolute tolerance
+    double aTol = 1.e-6;
+    
+    /// The relative tolerance
+    double rTol = 0.0;
+    
+    /// The maximal number of iterations
+    size_t maxIter = 1000;
+    
+    /// Print solver residuals every \ref printCycles 's iteration
+    size_t print_cycle = 100;
+    
+    /// Throw an error if tolerance could not be reached
+    bool breakTolNotReached = true;
+    
+    /// The solver verbosity level
+    int info = 0;
+    
+    // some parameters used internally  
+    double absRresidual = -1.0;
+    double relResidual = -1.0;
+    int iterations = -1;
+    int error = -1;
+  };
+
+} // end namespace AMDiS
+
+#include "LinearSolverInterface.inc.hpp"

dune/amdis/linear_algebra/LinearSolverInterface.inc.hpp

+#pragma once
+
+#include <dune/amdis/Initfile.hpp>
+
+namespace AMDiS {
+
+template <class M, class VX, class VB>
+LinearSolverInterface<M, VX, VB>::LinearSolverInterface(std::string prefix)
+  : prefix(prefix)
+{
+  // get creator string for solver
+  Parameters::get(prefix + "->name", solverName);
+  
+  // get creator string for preconditioner
+  Parameters::get(prefix + "->left precon", preconName);
+  if (preconName.empty() || preconName == "no")
+    Parameters::get(prefix + "->right precon", preconName);
+  if (preconName.empty() || preconName == "no")
+    Parameters::get(prefix + "->precon->name", preconName);
+      
+  Parameters::get(prefix + "->absolute tolerance", aTol);
+  Parameters::get(prefix + "->relative tolerance", rTol);
+  Parameters::get(prefix + "->max iteration", maxIter);
+  Parameters::get(prefix + "->print cycle", printCycle);
+  Parameters::get(prefix + "->break if tolerance not reached", breakTolNotReached);
+  Parameters::get(prefix + "->info", info);
+}
+
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp

@@ -19,9 +19,14 @@ namespace AMDiS
     using Self = BlockMTLMatrix;
     
   public:
+    /// The index/size - type
     using size_type  = typename MTLMatrix::size_type;
+    
+    /// The type of the elements of the MTLMatrix
     using value_type = typename MTLMatrix::value_type;
 
+  public:
+    /// Return the (R,C)'th matrix block
     template <size_t R, size_t C>
     auto& operator()(const index_<R>, const index_<C>)
     {
@@ -29,7 +34,7 @@ namespace AMDiS
       return std::get<C>(std::get<R>(*this));
     }
     
-    /// Return a shared pointer to the i'th underlaying base vector
+    /// Return the (R,C)'th matrix block
     template <size_t R, size_t C>
     auto const& operator()(const index_<R>, const index_<C>) const
     {
@@ -37,7 +42,10 @@ namespace AMDiS
       return std::get<C>(std::get<R>(*this));
     }
     
+    /// Return the number of row blocks
     static constexpr size_t N() { return _N; }
+    
+    /// Return the number of column blocks
     static constexpr size_t M() { return _M; }
 
     /// perform blockwise multiplication A*b -> x
@@ -51,23 +59,28 @@ namespace AMDiS
       
       For<0, _N>::loop([this, &r_rows, &r_cols, &b, &x](const auto _i) {
         bool first = true;
+        
+        // a reference to the i'th block of x
         VectorOut x_i(x[r_rows[_i]]);
         For<0, _M>::loop([this, &b, &x_i, &r_cols, &first, _i](const auto _j) {
-          auto const& block = this->operator()(_i, _j);
+          auto const& A_ij = this->operator()(_i, _j);
           if (num_rows(block) > 0) {
+            // a reference to the j'th block of b
             const VectorIn b_j(b[r_cols[_j]]);
             if (first) {
-              Assign::first_update(x_i, block * b_j);
+              Assign::first_update(x_i, A_ij * b_j);
               first = false;
             }
             else {
-              Assign::update(x_i, block * b_j);
+              Assign::update(x_i, A_ij * b_j);
             }
           }
         });
       });
     }
 
+    /// A Multiplication operator returns a multiplication-expresssion.
+    /// Calls \ref mult internally.
     template <class VectorIn>
     mtl::vector::mat_cvec_multiplier<Self, VectorIn> 
     operator*(VectorIn const& v) const
@@ -75,6 +88,8 @@ namespace AMDiS
       return {*this, v};
     }
     
+    /// Fill an array of irange corresponding to the row-sizes, used 
+    /// to access sub-vectors
     void getRowRanges(std::array<mtl::irange, _N>& r_rows) const
     {      
       size_t start = 0;
@@ -85,6 +100,8 @@ namespace AMDiS
       });
     }
     
+    /// Fill an array of irange corresponding to the column-sizes, used 
+    /// to access sub-vectors
     void getColRanges(std::array<mtl::irange, _M>& r_cols) const
     {      
       size_t start = 0;
@@ -95,6 +112,8 @@ namespace AMDiS
       });
     }
     
+    /// Fill two arrays of irange corresponding to row and column sizes.
+    /// \see getRowRanges() and \see getColRanges()
     void getRanges(std::array<mtl::irange, _N>& r_rows, std::array<mtl::irange, _M>& r_cols) const
     {      
       getRowRanges(r_rows);
@@ -102,6 +121,7 @@ namespace AMDiS
     }
   };
 
+  /// Return the number of overall rows of a BlockMTLMatrix
   template <class MTLMatrix, size_t _N, size_t _M>
   inline size_t num_rows(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
   {
@@ -112,6 +132,7 @@ namespace AMDiS
     return nRows;
   }
 
+  /// Return the number of overall columns of a BlockMTLMatrix
   template <class MTLMatrix, size_t _N, size_t _M>
   inline size_t num_cols(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
   {
@@ -122,12 +143,14 @@ namespace AMDiS
     return nCols;
   }
 
+  /// Return the size, i.e. rows*columns of a BlockMTLMatrix
   template <class MTLMatrix, size_t _N, size_t _M>
   inline size_t size(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
   {
     return num_rows(A) * num_cols(A);
   }
 
+  /// Nullify a BlockMTLMatrix, i.e. nullify each block.
   template <class MTLMatrix, size_t _N, size_t _M>
   inline void set_to_zero(BlockMTLMatrix<MTLMatrix, _N, _M>& A) 
   {
@@ -140,6 +163,8 @@ namespace AMDiS
 
 } // end namespace AMDiS
 
+
+/// \cond HIDDEN_SYMBOLS
 namespace mtl
 {
   template <class MTLMatrix, size_t _N, size_t _M>
@@ -160,3 +185,4 @@ namespace mtl
   } // end namespace ashape
 
 } // end namespace mtl
+/// \endcond

dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp

+#pragma once
+
+namespace AMDiS
+{
+  /// A simple block-vector class
+  template <class MTLVector, size_t _N>
+  class BlockMTLVector
+    : public std::array<MTLVector, _N>
+  {
+  public:
+    /// The type of the blocks
+    using BaseVector = MTLVector;
+    
+    /// The index/size - type
+    using size_type  = typename MTLVector::size_type;
+    
+    /// The type of the elements of the MTLVector
+    using value_type = typename MTLVector::value_type;
+  };
+
+  /// Return the number of overall rows of a BlockMTLVector
+  template <class MTLVector, size_t _N>
+  inline size_t num_rows(BlockMTLVector<MTLVector, _N> const& vec)
+  {
+    size_t nRows = 0;
+    For<0, _N>::loop([&](const auto _i) {
+      nRows += num_rows(std::get<_i>(vec));
+    });
+    return nRows;
+  }
+
+  /// Return the number of overall columns of a BlockMTLVector
+  template <class MTLVector, size_t _N,>
+  inline size_t num_cols(BlockMTLVector<MTLVector, _N> const& vec)
+  {
+    return 1;
+  }
+
+  /// Return the size, i.e. rows*columns of a BlockMTLVector
+  template <class MTLVector, size_t _N>
+  inline size_t size(BlockMTLVector<MTLVector, _N> const& vec)
+  {
+    return num_rows(vec);
+  }
+
+  /// Nullify a BlockMTLVector, i.e. nullify each block.
+  template <class MTLVector, size_t _N>
+  inline void set_to_zero(BlockMTLVector<MTLVector, _N>& vec) 
+  {
+    For<0, _N>::loop([&](const auto _i) {
+        set_to_zero(std::get<_i>(vec);
+    });
+  }
+  
+  
+  /// A wrapper, that creates a contiguos vector corresponding to a block-vector
+  /// and copy the value on construction ans possibly back on destruction.
+  template <class BlockVector, class Vector>
+  struct BlockVectorWrapper
+  {
+    BlockVectorWrapper(blockVector& blockVector, 
+                       bool copyBack = std::is_const<BlockVector>::value)
+      : blockVector(blockVector)
+      , vector(num_rows(blockVector))
+      , copyBack(copyBack)
+    {
+      assignTo();
+    }
+    
+    ~BlockVectorWrapper()
+    {
+      if (copyBack)
+        assignFrom(bool_<!std::is_const<BlockVector>::value>());
+    }
+    
+    /// Return a reference to the block-vector
+    BlockVector const&  getBlockVector() const  { return blockVector; }
+    
+    /// Return a reference to the contiguose-vector
+    Vector&       getVector()       { return vector; }
+    Vector const& getVector() const { return vector; }
+
+  private:
+    /// copy from block-vector to vector
+    void assignTo()
+    {
+      size_t start = 0;
+      for (size_t r = 0; r < blockVector.size(); ++r) {
+        size_t finish = start + num_rows(blockVector[r]);
+        mtl::irange range(start, finish);
+        
+        vector[range] = blockVector[r];
+      }
+    }
+    
+    /// do not copy vector to block-vector since this is const
+    template <bool Copy> 
+    std::enable_if_t<!Copy> assignFrom(bool_<Copy>) {}
+    
+    /// copy from vector to block-vector
+    template <bool Copy> 
+    std::enable_if_t<Copy> assignFrom(bool_<Copy>)
+    {
+      size_t start = 0;
+      for (size_t r = 0; r < blockVector.size(); ++r) {
+        size_t finish = start + num_rows(blockVector[r]);
+        mtl::irange range(start, finish);
+        
+        blockVector[r] = vector[range];
+      }
+    }
+    
+  private:
+    BlockVector& blockVector;
+    Vector vector;
+    
+    /// Copy data back to block-vector of destruction
+    bool copyBack;
+  };
+
+  
+  template <class BlockVector, class Vector = mtl::dense_vector<typename BlockVector::value_type>>
+  BlockVectorWrapper<BlockVector, Vector>
+  blockWrapper(BlockVector& bvec)
+  {
+    return {bvec};
+  }
+    
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/DOFMatrix.hpp

@@ -13,34 +13,49 @@
 
 namespace AMDiS
 {
+  /// \brief The basic container that stores a base matrix and a corresponding 
+  /// row/column feSpace
   template <class RowFeSpaceType, class ColFeSpaceType, 
             class ValueType = double>
   class DOFMatrix
   {
   public:
+    /// The type of the finite element space / basis of the row
     using RowFeSpace = RowFeSpaceType;
+    
+    /// The type of the finite element space / basis of the column
     using ColFeSpace = ColFeSpaceType;
+    
+    /// The vector type of the underlying base matrix
     using BaseMatrix = mtl::compressed2D<ValueType>;
     
+    /// The index/size - type
     using size_type  = typename BaseMatrix::size_type;
     
+    /// The type of the elements of the DOFMatrix
     using value_type = ValueType;
+    
+    /// The underlying field-type (typically the same as \ref value_type)
     using field_type = typename BaseMatrix::value_type;
     
+    /// The type of the inserter used when filling the matrix. NOTE: need not be public
     using Inserter = mtl::matrix::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
     
-    static constexpr int MIN_NNZ_PER_ROW = 20;
-    
-    /// Constructor. Constructs new BaseVector.
-    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
+  public:
+    /// Constructor. Constructs new BaseMatrix.
+    DOFMatrix(RowFeSpace const& rowFeSpace, 
+              ColFeSpace const& colFeSpace, 
+              std::string name)
       : rowFeSpace(rowFeSpace)
       , colFeSpace(colFeSpace)
       , name(name)
       , matrix(ClonablePtr<BaseMatrix>::make())
     {}
     
-    /// Constructor. Takes pointer of data-vector.
-    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name, 
+    /// Constructor. Takes pointer of data-matrix.
+    DOFMatrix(RowFeSpace const& rowFeSpace, 
+              ColFeSpace const& colFeSpace, 
+              std::string name, 
               BaseMatrix& matrix_ref)
       : rowFeSpace(rowFeSpace)
       , colFeSpace(colFeSpace)
@@ -60,37 +75,40 @@ namespace AMDiS
       return colFeSpace;
     }
     
-    /// Return the data-vector \ref vector
-    BaseMatrix const& getMatrix() const
+    /// Return a reference to the data-matrix \ref matrix
+    BaseMatrix& getMatrix()
     {
       return *matrix;
     }
     
-    /// Return the data-vector \ref vector
-    BaseMatrix& getMatrix()
+    /// Return a reference to the data-matrix \ref matrix
+    BaseMatrix const& getMatrix() const
     {
       return *matrix;
     }
     
-    /// Return the size of the \ref feSpace
+    /// Return the size of the row \ref feSpace
     size_type N() const
     {
       return rowFeSpace.size();
     }
     
+    /// Return the size of the column \ref feSpace
     size_type M() const
     {
       return colFeSpace.size();
     }
     
-    /// Return the \ref name of this vector
+    /// Return the \ref name of this matrix
     std::string getName() const
     {
       return name;