solver interface compiles and heat runs correctly

dune/amdis/LinearAlgebra.hpp

 #pragma once
 
+#include "linear_algebra/LinearSolverInterface.hpp"
+#include "linear_algebra/SolverInfo.hpp"
+
 #if defined(AMDIS_BACKEND_ISTL)
 
 #include "linear_algebra/istl/SystemVector.hpp"
@@ -11,6 +14,8 @@
 #include "linear_algebra/mtl/SystemVector.hpp"
 #include "linear_algebra/mtl/SystemMatrix.hpp"
 #include "linear_algebra/mtl/LinearSolver.hpp"
+#include "linear_algebra/mtl/ITL_Solver.hpp"
+#include "linear_algebra/mtl/BITL_Solver.hpp"
 
 #elif defined(AMDIS_BACKEND_PETSC)
 

dune/amdis/ProblemStat.hpp

@@ -61,8 +61,8 @@ namespace AMDiS
     using SystemMatrixType = SystemMatrix<FeSpaces>;
     
     using OperatorType = Operator<MeshView>;
-    using LinearSolverType = LinearSolver<typename SystemMatrixType::MultiMatrix,
-                                          typename SystemVectorType::MultiVector>;
+    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
+                                                   typename SystemVectorType::MultiVector>;
     
   public:
     /// \brief Constructor. Takes the name of the problem that is used to 
@@ -213,7 +213,12 @@ namespace AMDiS
     //
     void createSolver()
     {
-      linearSolver = std::make_shared<LinearSolverType>(name + "->solver");
+      using Creator = SolverCreator<typename SystemMatrixType::MultiMatrix,
+                                    typename SystemVectorType::MultiVector>;
+                                    
+      std::string solverName = "cg";
+      Parameters::get(name + "->solver->name", solverName);
+      linearSolver = Creator::create(solverName, name + "->solver");
     }
     
     //

dune/amdis/ProblemStat.inc.hpp

@@ -161,14 +161,39 @@ namespace AMDiS
                                      bool storeMatrixData)
   {
     Timer t;
+    
+    SolverInfo solverInfo(name + "->solver");
+    solverInfo.setCreateMatrixData(createMatrixData);
+    solverInfo.setStoreMatrixData(storeMatrixData);
+    
     solution->compress();
-    linearSolver->solve(adaptInfo, systemMatrix->getMatrix(), 
-                        solution->getVector(), rhs->getVector());
+    linearSolver->solve(systemMatrix->getMatrix(), 
+                        solution->getVector(), rhs->getVector(),
+                        solverInfo);
     
     AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
-
-//     adaptInfo.setSolverIterations(statistics.iterations);
-//     adaptInfo.setSolverResidual(statistics.reduction);
+    
+    if (solverInfo.getInfo() > 0) {
+      if (solverInfo.getAbsResidual() >= 0.0) {
+        if (solverInfo.getRelResidual() >= 0.0)
+          AMDIS_MSG("Residual norm: ||b-Ax|| = " << solverInfo.getAbsResidual()
+                        << ", ||b-Ax||/||b|| = " << solverInfo.getRelResidual());
+        else
+          AMDIS_MSG("Residual norm: ||b-Ax|| = " << solverInfo.getAbsResidual());
+      }
+    }
+    
+    if (solverInfo.doBreak()) {
+      std::stringstream tol_str;
+      if (solverInfo.getAbsTolerance() > 0
+          && solverInfo.getAbsResidual() > solverInfo.getAbsTolerance()) 
+        tol_str << "absTol = " << solverInfo.getAbsTolerance() << " ";
+      if (solverInfo.getRelTolerance() > 0
+          && solverInfo.getRelResidual() > solverInfo.getRelTolerance()) 
+        tol_str << "relTol = " << solverInfo.getRelTolerance() << " ";
+      AMDIS_ERROR_EXIT("Tolerance " << tol_str.str() << " could not be reached!");
+    }
+    
   }
   
   

dune/amdis/linear_algebra/LinearAlgebraBase.hpp

+#pragma once
+
+namespace AMDiS
+{
+  namespace Impl 
+  {
+    // Type-Trait needs to be specialized by linear algebra libraries
+    template <class Vector>
+    struct BaseVector { using type = Vector; };
+  }
+
+  template <class Vector>
+  using BaseVector = typename Impl::BaseVector<Vector>::type;
+    
+} // end namespace AMDiS

dune/amdis/linear_algebra/LinearSolverInterface.hpp

@@ -11,8 +11,11 @@
  * systems.
  */
 
-#include <dune/amdis/linear_algebra/Preconditioner.hpp>
-#include <dune/amdis/linear_algebra/SolverStatistics.hpp>
+#include <dune/amdis/linear_algebra/LinearAlgebraBase.hpp>
+
+#include <dune/amdis/linear_algebra/PreconditionerInterface.hpp>
+#include <dune/amdis/linear_algebra/RunnerInterface.hpp>
+#include <dune/amdis/linear_algebra/SolverInfo.hpp>
 
 namespace AMDiS
 {  
@@ -20,6 +23,11 @@ namespace AMDiS
   template <class Matrix, class VectorX, class VectorB = VectorX>
   class LinearSolverInterface
   {
+  protected:
+    using RunnerBase = RunnerInterface<Matrix, BaseVector<VectorX>, 
+                                               BaseVector<VectorB>>;
+    
+  public:
     /// Destructor.
     virtual ~LinearSolverInterface() {}
     
@@ -29,21 +37,23 @@ namespace AMDiS
      * 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 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.
      **/
     void solve(Matrix const& A, VectorX& x, VectorB const& b,
-               SolverStatistics& statistics)
+               SolverInfo& solverInfo)
     {      
-      solveImpl(A, x, b, statistics);
+      solveImpl(A, x, b, solverInfo);
     }
     
+    // return the runner/worker corresponding to this solver (optional)
+    virtual std::shared_ptr<RunnerBase> getRunner() {};
 
   private:
     /// main methods that all solvers must implement
     virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
-                           SolverStatistics& statistics) = 0;
+                           SolverInfo& solverInfo) = 0;
   };
 
 } // end namespace AMDiS

dune/amdis/linear_algebra/PreconditionerInterface.hpp

@@ -2,8 +2,8 @@
 
 namespace AMDiS
 {
-  /// base-class for Preconditioner types
-  template <class Matrix, class Vector>
+  /// Interface for Preconditioner types
+  template <class Matrix, class VectorX, class VectorB = VectorX>
   struct PreconditionerInterface
   {
     /// Virtual destructor.
@@ -16,10 +16,13 @@ namespace AMDiS
     virtual void exit() = 0;
 
     /// Apply the preconditioner to a vector \p x and store the result in \p y
-    virtual void solve(Vector const& x, Vector& y) const = 0;
+    virtual void solve(VectorX const& x, VectorX& y) const = 0;
     
     /// Apply the transposed preconditioner to a vector \p x and store the result in \p y
-    virtual void adjoint_solve(Vector const& x, Vector& y) const = 0;
+    virtual void adjointSolve(VectorX const& x, VectorX& y) const
+    {
+      AMDIS_ERROR_EXIT("Must be implemented by derived class.");
+    }
   };
 
 } // end namespace AMDiS

dune/amdis/linear_algebra/RunnerInterface.hpp

 #pragma once
 
-#include <dune/amdis/linear_algebra/SolverStatistics.hpp>
+#include <dune/amdis/linear_algebra/SolverInfo.hpp>
 
 namespace AMDiS
 {
-  /// base-class for Runner / Worker types
-  template <class Matrix, class Vector>
-  struct RunnerInterface
+  /// Interface for Runner / Worker types used in solver classes
+  template <class Matrix, class VectorX, class VectorB = VectorX>
+  class RunnerInterface
   {
-    // virtual destructor
+  protected:
+    using PreconBase = PreconditionerInterface<Matrix, VectorX, VectorB>;
+    
+  public:
+    /// virtual destructor
     virtual ~RunnerInterface() {}
 
     /// Is called at the beginning of a solution procedure
@@ -18,12 +22,19 @@ namespace AMDiS
     virtual void exit() = 0;
 
     /// Solve the system A*x = b
-    virtual int solve(Matrix const& A, Vector& x, Vector const& b, 
-                      SolverStatistics& stat) = 0;
+    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b, 
+                      SolverInfo& solverInfo) = 0;
 
     /// Solve the adjoint system A^T*x = b
-    virtual int adjoint_solve(Matrix const& A, Vector& x, Vector const& b, 
-                              SolverStatistics& stat) = 0;
+    virtual int adjointSolve(Matrix const& A, VectorX& x, VectorB const& b, 
+                             SolverInfo& solverInfo)
+    {
+      AMDIS_ERROR_EXIT("Must be implemented by derived class.");
+      return 0;
+    }
+                      
+    virtual std::shared_ptr<PreconBase> getLeftPrecon() {}
+    virtual std::shared_ptr<PreconBase> getRightPrecon() {}
   };
 
 } // end namespace AMDiS

dune/amdis/linear_algebra/SolverStatistics.hpp → dune/amdis/linear_algebra/SolverInfo.hpp

@@ -4,11 +4,12 @@
 
 namespace AMDiS
 {
-
-  class SolverStatistics
+  /// Class that stores information about the solution process, like tolerances 
+  /// and iteration counts and is passed to the solver and filled there.
+  class SolverInfo
   {
   public:
-    /// The constructor reads needed parameters and sets solvers \p prefix.
+    /// The constructor reads needed parameters and sets solver \p prefix.
     /**
      * Reads parameters for a solver with name 'NAME':
      *   NAME->absolute tolerance   \ref aTol
@@ -16,8 +17,13 @@ namespace AMDiS
      *   NAME->info                 \ref info
      *   NAME->break if tolerance not reached \ref breakTolNotReached
     **/
-    explicit SolverStatistics(std::string prefix);
-
+    explicit SolverInfo(std::string prefix)
+    {
+      Parameters::get(prefix + "->abolute tolerance", aTol);
+      Parameters::get(prefix + "->relative tolerance", rTol);
+      Parameters::get(prefix + "->info", info);
+      Parameters::get(prefix + "->break if tolerance not reached", breakTolNotReached);
+    }
 
     /** \name getting methods
      * \{
@@ -47,10 +53,16 @@ namespace AMDiS
       return info;
     }
 
-    /// Returns \ref residual
-    double getResidual() const
+    /// Returns \ref absResidual
+    double getAbsResidual() const
+    {
+      return absResidual;
+    }
+
+    /// Returns \ref relResidual
+    double getRelResidual() const
     {
-      return residual;
+      return relResidual;
     }
     
     /// Returns the initfile \ref prefix
@@ -59,6 +71,25 @@ namespace AMDiS
       return prefix;
     }
 
+    /// Returns \ref createMatrixData
+    bool doCreateMatrixData() const
+    {
+      return createMatrixData;
+    }
+
+    /// Returns \ref storeMatrixData
+    bool doStoreMatrixData() const
+    {
+      return storeMatrixData;
+    }
+    
+    bool doBreak() const
+    {
+      return breakTolNotReached && (
+          (aTol > 1.e-30 && absResidual > aTol) || 
+          (rTol > 1.e-30 && relResidual > rTol) );
+    }
+
     /** \} */
 
 
@@ -73,46 +104,83 @@ namespace AMDiS
     }
 
     /// Sets \ref rTol
-    void setAbsTolerance(double tol)
+    void setRelTolerance(double tol)
     {
       rTol = tol;
     }
 
+    /// Sets \ref aTol
+    void setAbsResidual(double r)
+    {
+      absResidual = r;
+    }
+
+    /// Sets \ref rTol
+    void setRelResidual(double r)
+    {
+      relResidual = r;
+    }
+
     /// Sets \ref info
     void setInfo(int i)
     {
       info = i;
     }
 
-    /// Sets \ref info
+    /// Sets \ref error
     void setError(int e)
     {
       error = e;
     }
 
+    /// Sets \ref createMatrixData
+    void setCreateMatrixData(bool b)
+    {
+      createMatrixData = b;
+    }
+
+    /// Sets \ref storeMatrixData
+    void setStoreMatrixData(bool b)
+    {
+      storeMatrixData = b;
+    }
+    
+    
+
     /** \} */
     
   private:
     /// The initfile prefix to read parameters
     std::string prefix;
     
-  private:
     /// The abolute tolerance
-    double aTol = 1.e-6;
+    double aTol = 0;
     
     /// The relative tolerance
-    double rTol = 0.0;
+    double rTol = 1.e-6;
     
     /// Throw an error if tolerance could not be reached
-    bool breakTolNotReached = true;
+    bool breakTolNotReached = false;
     
     /// The solver verbosity level
     int info = 0;
     
-    // some parameters used internally  
-    double absRresidual = -1.0;
+    /// The absolute residual, default (-1): not set
+    double absResidual = -1.0;
+    
+    /// The relative residual, relative to the two_norm(b), default (-1): not set
     double relResidual = -1.0;
+    
+    /// The error-code, default (-1): not set
     int error = -1;
+    
+    /// If true, the matrix will be initialized and the 
+    /// corresponding runner of the system receives the 
+    /// matrix in the init() method.
+    bool createMatrixData = true;
+    
+    /// If false, the exit() method of the runner will be called.
+    bool storeMatrixData = false;
   };
 
 } // end namespace AMDiS

dune/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp

+#pragma once
+
+#include <dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
+#include "itl/block_diagonal.hpp"
+
+namespace AMDiS
+{
+
+  template <class SubMatrix, size_t N, size_t M, class Vector>
+  struct PreconditionerCreator< BlockMTLMatrix<SubMatrix, N, M>, Vector >
+  {
+    using Matrix = BlockMTLMatrix<SubMatrix, N, M>;
+    using PreconBase = PreconditionerInterface<Matrix, Vector>;
+    
+    template <template<class> class ITLPrecon>
+    using Precon = Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>;
+    
+    static std::unique_ptr<PreconBase> create(std::string name)
+    { 
+      if (name == "diag" || name == "jacobi")
+        return std::make_unique<Precon<DiagonalPreconditioner>>();
+      else // by default, do nothing
+        return std::make_unique<Precon<IdentityPreconditioner>>();
+    }
+  };
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/BITL_Solver.hpp

+#pragma once
+
+#include <dune/amdis/linear_algebra/mtl/ITL_Solver.hpp>
+
+namespace AMDiS
+{
+
+  template <class SubMatrix, size_t N, size_t M, class Vector>
+  struct SolverCreator< BlockMTLMatrix<SubMatrix, N, M>, Vector >
+  {
+    using Matrix = BlockMTLMatrix<SubMatrix, N, M>;
+    using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
+    
+    template <class ITLSolver>
+    using Solver = LinearSolver<Matrix, Vector, 
+                      KrylovRunner<ITLSolver, Matrix, BaseVector<Vector>>>;
+    
+    /// Instantiate a new linear solver
+    static std::unique_ptr<SolverBase> create(std::string name, std::string prefix)
+    { 
+      if (name == "cg")
+        return std::make_unique<Solver<cg_solver_type>>(prefix);
+      else if (name == "cgs")
+        return std::make_unique<Solver<cgs_solver_type>>(prefix);
+      else if (name == "bicgstab" || name == "bcgs")
+        return std::make_unique<Solver<bicgstab_type>>(prefix);
+      else if (name == "tfqmr")
+        return std::make_unique<Solver<tfqmr_solver_type>>(prefix);
+      else if (name == "bicgstab_ell" || name == "bcgsl")
+        return std::make_unique<Solver<bicgstab_ell_type>>(prefix);
+      else if (name == "gmres")
+        return std::make_unique<Solver<gmres_type>>(prefix);
+      else if (name == "fgmres")
+        return std::make_unique<Solver<fgmres_type>>(prefix);
+      else if (name == "minres")
+        return std::make_unique<Solver<minres_solver_type>>(prefix);
+      else if (name == "gcr")
+        return std::make_unique<Solver<gcr_type>>(prefix);
+//       else if (name == "preonly")
+//         return std::make_unique<Solver<preonly_type>>(prefix);
+      else
+        AMDIS_ERROR_EXIT("Unknown Solver-name!");
+    }
+  };
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp

@@ -64,7 +64,7 @@ namespace AMDiS
         VectorOut x_i(x[r_rows[_i]]);
         For<0, _M>::loop([this, &b, &x_i, &r_cols, &first, _i](const auto _j) {
           auto const& A_ij = this->operator()(_i, _j);
-          if (num_rows(block) > 0) {
+          if (num_rows(A_ij) > 0) {
             // a reference to the j'th block of b
             const VectorIn b_j(b[r_cols[_j]]);
             if (first) {

dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp

 #pragma once
 
+#include <dune/amdis/linear_algebra/LinearAlgebraBase.hpp>
+
 namespace AMDiS
 {
   /// A simple block-vector class
@@ -17,6 +19,18 @@ namespace AMDiS
     /// The type of the elements of the MTLVector
     using value_type = typename MTLVector::value_type;
   };
+  
+  
+  namespace Impl
+  {
+    /// Specialization of Impl::BaseVector from \file LinearSolverInterface.hpp
+    template <class MTLVector, size_t N>
+    struct BaseVector<BlockMTLVector<MTLVector, N>>
+    {
+      using type = MTLVector;
+    };
+  }
+  
 
   /// Return the number of overall rows of a BlockMTLVector
   template <class MTLVector, size_t _N>
@@ -30,7 +44,7 @@ namespace AMDiS
   }
 
   /// Return the number of overall columns of a BlockMTLVector
-  template <class MTLVector, size_t _N,>
+  template <class MTLVector, size_t _N>
   inline size_t num_cols(BlockMTLVector<MTLVector, _N> const& vec)
   {
     return 1;
@@ -48,7 +62,7 @@ namespace AMDiS
   inline void set_to_zero(BlockMTLVector<MTLVector, _N>& vec) 
   {
     For<0, _N>::loop([&](const auto _i) {
-        set_to_zero(std::get<_i>(vec);
+        set_to_zero(std::get<_i>(vec));
     });
   }
   
@@ -58,7 +72,7 @@ namespace AMDiS
   template <class BlockVector, class Vector>
   struct BlockVectorWrapper
   {
-    BlockVectorWrapper(blockVector& blockVector, 
+    BlockVectorWrapper(BlockVector& blockVector, 
                        bool copyBack = std::is_const<BlockVector>::value)
       : blockVector(blockVector)
       , vector(num_rows(blockVector))

dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp

@@ -15,12 +15,7 @@
 namespace AMDiS
 {
   template <class Matrix, class Vector>
-  struct PreconditionerCreator
-  {
-    using PreconBase = PreconditionerInterface<Matrix, Vector>;
-  
-    static std::unique_ptr<PreconBase> create(std::string name, std::string prefix);
-  };
+  struct PreconditionerCreator;
 
   /**
    * \ingroup Solver
@@ -82,27 +77,30 @@ namespace AMDiS
   template <class Matrix>
   using ICPreconditioner = itl::pc::ic_0<Matrix>;
 
-                       
-                       
-  template <class Matrix, class Vector>
-  std::unique_ptr<typename PreconditionerCreator<Matrix, Vector>::PreconBase> 
-  PreconditionerCreator<Matrix, Vector>::create(std::string name)
-  {    
-    template <template<class> class ITLPrecon>
-    using Type = typename Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>;
+               
+  template <class T, class Param, class Vector>
+  struct PreconditionerCreator<mtl::compressed2D<T, Param>, Vector>
+  {
+    using Matrix = mtl::compressed2D<T, Param>;
+    using PreconBase = PreconditionerInterface<Matrix, Vector>;
     
-    if (name == "diag" || name == "jacobi")
-      return std::make_unique<Type<DiagonalPreconditioner>>();
-    else if (name == "masslumping")
-      return std::make_unique<Type<cMassLumpingPreconditioner>>();
-    else if (name == "ilu" || name == "ilu0")
-      return std::make_unique<Type<ILUPreconditioner>>();
-    else if (name == "ic" || name == "ic0")
-      return std::make_unique<Type<ICPreconditioner>>();
-    else // by default, do nothing
-      return std::make_unique<Type<IdentityPreconditioner>>();
-  }
-
+    template <template<class> class ITLPrecon>
+    using Precon = Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>;
+  
+    static std::unique_ptr<PreconBase> create(std::string name)
+    { 
+      if (name == "diag" || name == "jacobi")
+        return std::make_unique<Precon<DiagonalPreconditioner>>();
+      else if (name == "masslumping")
+        return std::make_unique<Precon<MassLumpingPreconditioner>>();
+      else if (name == "ilu" || name == "ilu0")
+        return std::make_unique<Precon<ILUPreconditioner>>();
+      else if (name == "ic" || name == "ic0")
+        return std::make_unique<Precon<ICPreconditioner>>();
+      else // by default, do nothing
+        return std::make_unique<Precon<IdentityPreconditioner>>();
+    }
+  };