interfaces for runner and preconditioner cleaned up

dune/amdis/CreatorMap.hh

@@ -3,9 +3,7 @@ namespace AMDiS
   template <class BaseClass>
   void CreatorMap<BaseClass>::clear()
   {
-    typename std::map<std::string, CreatorInterface<BaseClass>*>::iterator it;
-
-    for (it = creatorMap.begin(); it != creatorMap.end(); ++it)
+    for (auto it = creatorMap.begin(); it != creatorMap.end(); ++it)
       delete it->second;
   }
 

dune/amdis/CreatorMap.hpp

+#pragma once
+
+// std c++ headers
+#include <map>
+
+// AMDiS includes
+#include "CreatorInterface.hpp"
+
+namespace AMDiS
+{
+  /** \ingroup Common
+   * \brief
+   * A CreatorMap is used to construct objects, which types depends on key words
+   * determined at run time. For example the LinearSolverInterfaceMap can create the different
+   * solver types depending on the solver parameter of the init file. The benefit
+   * of such creator maps is, that you can extend them only by writing an creator
+   * class for your own new class and give the creator together with a key word
+   * to the map.
+   */
+  template <class BaseClass>
+  class CreatorMap
+  {
+  public:
+    using CreatorMapType = std::map<std::string, CreatorInterface<BaseClass>*>;
+
+  public:
+    /// Adds a new creator together with the given key to the map.
+    static void addCreator(std::string key, CreatorInterface<BaseClass>* creator)
+    {
+      init();
+      AMDIS_TEST_EXIT(creatorMap[key] == NULL,
+        "There is already a creator for key " << key);
+      creatorMap[key] = creator;
+    }
+
+    static void addCreator(std::string backend, 
+                           std::string key, 
+                           CreatorInterface<BaseClass>* creator)
+    {
+      addCreator(backend + "_" + key, creator);
+    }
+
+    /// Creates a object of the type corresponding to key.
+    static CreatorInterface<BaseClass>* getCreator(std::string key, 
+                                                   std::string initFileStr)
+    {
+      init();
+      
+      CreatorInterface<BaseClass>* creator = creatorMap[key];
+      AMDIS_TEST_EXIT(creator,
+        "No creator for key \"" << key << "\" defined in init file for parameter \"" << initFileStr << "\"");
+
+      return creator;
+    }
+
+    static void clear()
+    {
+      for (auto entries : creatorMap)
+        delete it->second;
+    }
+
+    static void addDefaultCreators();
+
+  protected:
+    /// Constructor is protected because derived maps should be singleton.
+    static void init()
+    {
+      if (!initialized)
+      {
+        initialized = true;
+        NullCreator<BaseClass>* nullCreator = new NullCreator<BaseClass>;
+        addCreator("0", nullCreator);
+        addDefaultCreators();
+      }
+    }
+
+  protected:
+    /// STL map containing the pairs of keys and creators.
+    static CreatorMapType creatorMap;
+
+    static bool initialized;
+  };
+
+  template <class BaseClass>
+  typename CreatorMap<BaseClass>::CreatorMapType CreatorMap<BaseClass>::creatorMap;
+
+  template <class BaseClass>
+  bool CreatorMap<BaseClass>::initialized = false;
+
+} // end namespace AMDiS
+
+#include "CreatorMap.hh"

dune/amdis/Log.hpp

@@ -54,7 +54,7 @@
 #define AMDIS_WARNING(text) AMDIS_MSG("Warning: " << text)
 
   /// if test is false, an error message is printed
-#define AMDIS_TEST(test, msg) if (!(test)) { AMDIS_ERROR(msg); }
+#define AMDIS_TEST(test, msg) if (!(test)) { AMDIS_WARNING(msg); }
 
   /// if test is false, an error message is printed and the program exits
 #define AMDIS_TEST_EXIT(test, msg) if (!(test)) { AMDIS_ERROR_EXIT(msg); }
@@ -62,10 +62,12 @@
   /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
 #ifdef NDEBUG
   #define AMDIS_TEST_EXIT_DBG(test, msg)
+  #define AMDIS_TEST_DBG(test, msg)
   #define AMDIS_MSG_DBG(text)
   #define AMDIS_DBG_VAR(var)
 #else
   #define AMDIS_TEST_EXIT_DBG(test, msg) AMDIS_TEST_EXIT(test, msg)
+  #define AMDIS_TEST_DBG(test, msg) AMDIS_TEST(test, msg)
   #define AMDIS_MSG_DBG(text) AMDIS_MSG(text)
   #define AMDIS_DBG_VAR(var) var
 #endif

dune/amdis/linear_algebra/LinearSolverInterface.hpp

@@ -11,53 +11,18 @@
  * 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
-  // ---------------------------------------------------------------------------
+#include <dune/amdis/linear_algebra/Preconditioner.hpp>
+#include <dune/amdis/linear_algebra/SolverStatistics.hpp>
 
+namespace AMDiS
+{  
+  /// Abstract base class for linear 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
@@ -67,158 +32,19 @@ namespace AMDiS
      *  \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)
+               SolverStatistics& statistics)
     {      
-      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;
+      solveImpl(A, x, b, statistics);
     }
     
-    /// 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;
+                           SolverStatistics& statistics) = 0;
   };
 
 } // 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/PreconditionerInterface.hpp

+#pragma once
+
+namespace AMDiS
+{
+  /// base-class for Preconditioner types
+  template <class Matrix, class Vector>
+  struct PreconditionerInterface
+  {
+    /// Virtual destructor.
+    virtual ~PreconditionerInterface() {}
+    
+    /// Is called a the beginning of a solution procedure
+    virtual void init(Matrix const& fullMatrix) = 0;
+    
+    /// Is called at the end of a solution procedure
+    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;
+    
+    /// 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;
+  };
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/RunnerInterface.hpp

+#pragma once
+
+#include <dune/amdis/linear_algebra/SolverStatistics.hpp>
+
+namespace AMDiS
+{
+  /// base-class for Runner / Worker types
+  template <class Matrix, class Vector>
+  struct RunnerInterface
+  {
+    // virtual destructor
+    virtual ~RunnerInterface() {}
+
+    /// Is called at the beginning of a solution procedure
+    virtual void init(Matrix const& A) = 0;
+
+    /// Is called at the end of a solution procedure
+    virtual void exit() = 0;
+
+    /// Solve the system A*x = b
+    virtual int solve(Matrix const& A, Vector& x, Vector const& b, 
+                      SolverStatistics& stat) = 0;
+
+    /// Solve the adjoint system A^T*x = b
+    virtual int adjoint_solve(Matrix const& A, Vector& x, Vector const& b, 
+                              SolverStatistics& stat) = 0;
+  };
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/SolverStatistics.hpp

+#pragma once
+
+#include <string>
+
+namespace AMDiS
+{
+
+  class SolverStatistics
+  {
+  public:
+    /// 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->info                 \ref info
+     *   NAME->break if tolerance not reached \ref breakTolNotReached
+    **/
+    explicit SolverStatistics(std::string prefix);
+
+
+    /** \name getting methods
+     * \{
+     */
+
+    /// Returns \ref aTol
+    double getAbsTolerance() const
+    {
+      return aTol;
+    }
+
+    /// Returns \ref rTol
+    double getRelTolerance() const
+    {
+      return rTol;
+    }
+
+    /// Returns error code in last run of an iterative solver
+    int getErrorCode() const
+    {
+      return error;
+    }
+
+    /// Returns info
+    int getInfo() const
+    {
+      return info;
+    }
+
+    /// 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 info
+    void setInfo(int i)
+    {
+      info = i;
+    }
+
+    /// Sets \ref info
+    void setError(int e)
+    {
+      error = e;
+    }
+
+    /** \} */
+    
+  private:
+    /// The initfile prefix to read parameters
+    std::string prefix;
+    
+  private:
+    /// The abolute tolerance
+    double aTol = 1.e-6;
+    
+    /// The relative tolerance
+    double rTol = 0.0;
+    
+    /// 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 error = -1;
+  };
+
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp

+/** \file ITL_Preconditioner.h */
+
+#pragma once
+
+// MTL4 headers
+#include <boost/numeric/itl/itl.hpp>
+#include <boost/numeric/itl/pc/ilu_0.hpp>
+#include <boost/numeric/itl/pc/ic_0.hpp>
+#include <boost/numeric/mtl/vector/assigner.hpp>
+
+#include "itl/masslumping.hpp"
+
+#include <dune/amdis/linear_algebra/mtl/Preconditioner.hpp>
+
+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);
+  };
+
+  /**
+   * \ingroup Solver
+   * \class AMDiS::DiagonalPreconditioner
+   * \brief ITL_Preconditioner implementation of diagonal (jacobi) preconditioner,
+   * \implements ITL_Preconditioner
+   *
+   * Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M=diag(A) \f$.
+   */
+  template <class Matrix>
+  using DiagonalPreconditioner = itl::pc::diagonal<Matrix>;
+
+  /**
+   * \ingroup Solver
+   * \class AMDiS::DiagonalPreconditioner
+   * \brief ITL_Preconditioner implementation of diagonal (jacobi) preconditioner,
+   * \implements ITL_Preconditioner
+   *
+   * Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M_ii=sum_j(A_ij) \f$.
+   */
+  template <class Matrix>
+  using MassLumpingPreconditioner = itl::pc::masslumping<Matrix>;
+
+
+  /**
+   * \ingroup Solver
+   * \class AMDiS::IdentityPreconditioner
+   * \brief ITL_Preconditioner implementation of identity preconditioner,
+   * \implements ITL_Preconditioner
+   *
+   * Identity preconditioner. Behaves like no preconditioning.
+   */
+  template <class Matrix>
+  using IdentityPreconditioner = itl::pc::identity<Matrix>;
+
+
+  /**
+   * \ingroup Solver
+   * \class AMDiS::ILUPreconditioner