Eigen solvers cleaned up

src/amdis/LinearAlgebra.hpp

 #pragma once
 
-#include <amdis/linear_algebra/LinearSolver.hpp>
-#include <amdis/linear_algebra/SolverInfo.hpp>
-
 #if HAVE_MTL
 
 #include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
 #include <amdis/linear_algebra/mtl/DOFVector.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
 
 #elif HAVE_EIGEN
 
 #include <amdis/linear_algebra/eigen/DOFMatrix.hpp>
 #include <amdis/linear_algebra/eigen/DOFVector.hpp>
-#include <amdis/linear_algebra/eigen/SolverCreator.hpp>
 
 #else // ISTL
 
 #include <amdis/linear_algebra/istl/DOFMatrix.hpp>
 #include <amdis/linear_algebra/istl/DOFVector.hpp>
-#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
-#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
 
 #endif
\ No newline at end of file

src/amdis/LinearSolvers.hpp

+#pragma once
+
+#include <amdis/linear_algebra/LinearSolver.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+
+#if HAVE_MTL
+
+#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
+#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
+
+#elif HAVE_EIGEN
+
+#include <amdis/linear_algebra/eigen/SolverCreator.hpp>
+
+#else // ISTL
+
+#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
+
+#endif
\ No newline at end of file

src/amdis/ProblemStat.cpp

@@ -6,5 +6,7 @@
 namespace AMDiS
 {
   // explicit template instatiation
-  template class ProblemStat<YaspGridBasis<2,1,2>>;
+  // template class ProblemStat<YaspGridBasis<2,1>>;
+  // template class ProblemStat<YaspGridBasis<2,2>>;
+  // template class ProblemStat<YaspGridBasis<2,1,2>>;
 } // end namespace AMDiS

src/amdis/ProblemStat.hpp

@@ -15,10 +15,12 @@
 #include <amdis/AdaptInfo.hpp>
 #include <amdis/Assembler.hpp>
 #include <amdis/CreatorInterface.hpp>
+#include <amdis/CreatorMap.hpp>
 #include <amdis/DirichletBC.hpp>
 #include <amdis/Flag.hpp>
 #include <amdis/Initfile.hpp>
 #include <amdis/LinearAlgebra.hpp>
+#include <amdis/LinearSolvers.hpp>
 #include <amdis/LocalAssemblerList.hpp>
 #include <amdis/Marker.hpp>
 #include <amdis/Mesh.hpp>
@@ -397,7 +399,9 @@ namespace AMDiS
 
 
 #ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
-  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
+  // extern template class ProblemStat<YaspGridBasis<2,1>>;
+  // extern template class ProblemStat<YaspGridBasis<2,2>>;
+  // extern template class ProblemStat<YaspGridBasis<2,1,2>>;
 #endif
 
 } // end namespace AMDiS

src/amdis/QuadratureFactory.hpp

@@ -69,7 +69,7 @@ namespace AMDiS
 
   struct PreQuadFactoryFromLocalFunction {};
 
-  auto makePreQuadratureFactory()
+  inline auto makePreQuadratureFactory()
   {
     return PreQuadFactoryFromLocalFunction{};
   }
@@ -125,7 +125,7 @@ namespace AMDiS
     Dune::QuadratureType::Enum qt;
   };
 
-  auto makePreQuadratureFactory(int order, Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
+  inline auto makePreQuadratureFactory(int order, Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
   {
     return PreQuadFactoryFromOrder{order, qt};
   }

src/amdis/linear_algebra/DOFMatrixBase.hpp

@@ -71,7 +71,6 @@ namespace AMDiS
     /// If \p setToZero is true, the matrix is set to 0
     void init(bool setToZero)
     {
-      timeInsertion_ = 0;
       backend_.init(*rowBasis_, *colBasis_, setToZero);
     }
 
@@ -79,7 +78,6 @@ namespace AMDiS
     void finish()
     {
       backend_.finish();
-      msg("  time(insertion) = {} sec", timeInsertion_);
     }
 
     /// Insert a block of values into the matrix (add to existing values)
@@ -88,7 +86,6 @@ namespace AMDiS
                 typename ColBasis::LocalView const& colLocalView,
                 ElementMatrix const& elementMatrix)
     {
-      Dune::Timer t;
       for (size_type i = 0; i < rowLocalView.size(); ++i) {
         size_type const row = flatMultiIndex(rowLocalView.index(i));
         for (size_type j = 0; j < colLocalView.size(); ++j) {
@@ -98,7 +95,6 @@ namespace AMDiS
           }
         }
       }
-      timeInsertion_ += t.elapsed();
     }
 
     /// Insert a single value into the matrix (add to existing value)

src/amdis/linear_algebra/RunnerInterface.hpp

 #pragma once
 
 #include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/SolverInfo.hpp>
+#include <amdis/linear_algebra/PreconditionerInterface.hpp>
 
 namespace AMDiS
 {
+  class SolverInfo;
+
   /// Interface for Runner / Worker types used in solver classes
   template <class Matrix, class VectorX, class VectorB = VectorX>
   class RunnerInterface

src/amdis/linear_algebra/eigen/DOFMatrix.hpp

@@ -63,11 +63,8 @@ namespace AMDiS
     void init(RowBasis const& rowBasis, ColBasis const& colBasis,
               bool prepareForInsertion)
     {
-      size_type nnzPerRow = matrix_.rows() == 0 ? 50 :
-        matrix_.nonZeros()/double(matrix_.rows());
       matrix_.resize(rowBasis.dimension(), colBasis.dimension());
       if (prepareForInsertion) {
-        tripletList_.reserve(2 * nnzPerRow * rowBasis.dimension());
         matrix_.setZero(); // maybe not necessary
       }
     }
@@ -76,13 +73,8 @@ namespace AMDiS
     /// final construction of the matrix.
     void finish()
     {
-      Dune::Timer t;
       matrix_.setFromTriplets(tripletList_.begin(), tripletList_.end());
-      msg("  time(setFromTriplets) = {} sec", t.elapsed());
-
-      t.reset();
       matrix_.makeCompressed();
-      msg("  time(makeCompressed) = {} sec", t.elapsed());
 
       tripletList_.clear(); // NOTE: this does not free the memory
     }
@@ -97,28 +89,23 @@ namespace AMDiS
     /// a one on the diagonal of the matrix.
     auto applyDirichletBC(std::vector<char> const& dirichletNodes)
     {
-      Dune::Timer t;
       for (std::size_t i = 0; i < dirichletNodes.size(); ++i) {
         if (dirichletNodes[i] == 1)
-          clearDirichletRow(i);
+          clearDirichletRow(i, int_<Orientation>);
       }
-      msg("  time(applyDirichletBC) = {} sec", t.elapsed());
       return std::list<Triplet<value_type>>{};
     }
 
 
   protected:
-    void clearDirichletRow(size_type idx)
-    {
-      clearDirichletRow(idx, int_<Orientation>);
-    }
-
     void clearDirichletRow(size_type idx, int_t<Eigen::ColMajor>)
     {
       for (size_type k = 0; k < matrix_.outerSize(); ++k) {
         for (typename BaseMatrix::InnerIterator it(matrix_, k); it; ++it) {
-          if (it.row() == idx)
+          if (it.row() == idx) {
             it.valueRef() = (it.row() == it.col() ? value_type(1) : value_type(0));
+            break;
+          }
         }
       }
     }

src/amdis/linear_algebra/eigen/DirectRunner.hpp

@@ -33,9 +33,12 @@ namespace AMDiS
     virtual void init(Matrix const& A) override
     {
       solver_.compute(A);
-      test_exit(solver_.info() == Eigen::Success, "Error in solver.compute(matrix)");
+
+      test_exit(solver_.info() == Eigen::Success,
+        "Error in solver.compute(matrix)");
     }
 
+
     /// Implementes \ref RunnerInterface::exit()
     virtual void exit() override {}
 

src/amdis/linear_algebra/eigen/IterativeRunner.hpp

@@ -29,7 +29,9 @@ namespace AMDiS
     virtual void init(Matrix const& A) override
     {
       solver_.compute(A);
-      test_exit(solver_.info() == Eigen::Success, "Error in solver.compute(matrix)");
+
+      test_exit(solver_.info() == Eigen::Success,
+        "Error in solver.compute(matrix)");
     }
 
     /// Implementes \ref RunnerInterface::exit()

src/amdis/linear_algebra/eigen/PreconConfig.hpp

@@ -19,14 +19,13 @@ namespace AMDiS
   {
     static void init(std::string const& prefix, Eigen::IncompleteLUT<T,I>& solver)
     {
-      double dropTol = -1;
-      Parameters::get(prefix + "->drop tolerance", dropTol);
-      if (dropTol >= 0)
-        solver.setDroptol(dropTol);
+      auto dropTol = Parameters::get<double>(prefix + "->drop tolerance");
+      if (dropTol)
+        solver.setDroptol(dropTol.value());
 
-      int fillFactor = 10;
-      Parameters::get(prefix + "->fill factor", fillFactor);
-      solver.setFillfactor(fillFactor);
+      auto fillFactor = Parameters::get<int>(prefix + "->fill factor");
+      if (fillFactor)
+        solver.setFillfactor(fillFactor.value());
     }
   };
 

src/amdis/linear_algebra/eigen/SolverConfig.hpp

@@ -2,9 +2,8 @@
 
 #include <string>
 
-#include <Eigen/SparseLU>
-//#include <Eigen/SuperLUSupport>
 #include <Eigen/UmfPackSupport>
+#include <unsupported/Eigen/IterativeSolvers>
 
 #include <amdis/Initfile.hpp>
 
@@ -22,7 +21,48 @@ namespace AMDiS
     }
 
     // fall-back default implementation
-    static void init(std::string const&, ...) {}
+    template <class D>
+    static void init(std::string const&, Eigen::SparseSolverBase<D>&) {}
+  };
+
+  template <class M, class P>
+  struct SolverConfig<Eigen::GMRES<M,P>>
+  {
+    using Base = Eigen::IterativeSolverBase<Eigen::GMRES<M,P>>;
+
+    static void init(std::string const& prefix, Eigen::GMRES<M,P>& solver)
+    {
+      SolverConfig<Base>::init(prefix, solver);
+
+      auto restart = Parameters::get<int>(prefix + "->restart");
+      if (restart)
+        solver.set_restart(restart.value());
+    }
+  };
+
+  template <class M, class P>
+  struct SolverConfig<Eigen::DGMRES<M,P>>
+  {
+    using Base = Eigen::IterativeSolverBase<Eigen::DGMRES<M,P>>;
+
+    static void init(std::string const& prefix, Eigen::DGMRES<M,P>& solver)
+    {
+      SolverConfig<Base>::init(prefix, solver);
+
+      auto restart = Parameters::get<int>(prefix + "->restart");
+      if (restart)
+        solver.set_restart(restart.value());
+
+      // number of eigenvalues to deflate at each restart
+      auto neigs = Parameters::get<int>(prefix + "->num eigenvalues");
+      if (neigs)
+        solver.setEigenv(neigs.value());
+
+      // maximum size of the deflation subspace
+      auto maxNeig = Parameters::get<int>(prefix + "->max num eigenvalues");
+      if (maxNeig)
+        solver.setMaxEigenv(maxNeig.value());
+    }
   };
 
 #if HAVE_SUITESPARSE_UMFPACK
@@ -37,20 +77,4 @@ namespace AMDiS
   };
 #endif
 
-  template <class M, class O>
-  struct SolverConfig<Eigen::SparseLU<M,O>>
-  {
-    static void init(std::string const& prefix, Eigen::SparseLU<M,O>& solver)
-    {
-      bool symmetric = false;
-      Parameters::get(prefix + "->symmetric pattern", symmetric);
-      solver.isSymmetric(symmetric);
-
-      double pivotThreshold = -1;
-      Parameters::get(prefix + "->pivot threshold", pivotThreshold);
-      if (pivotThreshold >= 0)
-        solver.setPivotThreshold(pivotThreshold);
-    }
-  };
-
 } // end namespace AMDiS

src/amdis/linear_algebra/eigen/SolverCreator.hpp

@@ -5,12 +5,14 @@
 #include <Eigen/IterativeLinearSolvers>
 #include <Eigen/MetisSupport>
 #include <Eigen/SparseLU>
-//#include <Eigen/SuperLUSupport>
+#include <Eigen/SuperLUSupport>
 #include <Eigen/UmfPackSupport>
+#include <unsupported/Eigen/IterativeSolvers>
 
 #include <amdis/CreatorMap.hpp>
 #include <amdis/Initfile.hpp>
 
+#include <amdis/linear_algebra/LinearSolver.hpp>
 #include <amdis/linear_algebra/eigen/DirectRunner.hpp>
 #include <amdis/linear_algebra/eigen/IterativeRunner.hpp>
 
@@ -31,27 +33,49 @@ namespace AMDiS
     virtual std::shared_ptr<SolverBase> create(std::string prefix) override
     {
       // get creator string for preconditioner
-      std::string preconName = "no";
-      Parameters::get(prefix + "->precon->name", preconName);
+      std::string precon = "no";
+      Parameters::get(prefix + "->precon->name", precon);
 
-      if (preconName == "diag" ||
-          preconName == "jacobi")
+      if (precon == "diag" ||
+          precon == "jacobi")
       {
         auto creator = SolverCreator<Eigen::DiagonalPreconditioner<Scalar>>{};
         return creator.create(prefix);
       }
-      else if (preconName == "ilu" ||
-               preconName == "ilu0" ||
-               preconName == "ilut")
+      else if (precon == "ilu")
       {
         auto creator = SolverCreator<Eigen::IncompleteLUT<Scalar>>{};
         return creator.create(prefix);
       }
+      else if (precon == "ic")
+      {
+        return createIncompleteCholesky(prefix);
+      }
       else {
         auto creator = SolverCreator<Eigen::IdentityPreconditioner>{};
         return creator.create(prefix);
       }
     }
+
+
+    template <class Ordering>
+    using IncompleteCholesky =
+      SolverCreator<Eigen::IncompleteCholesky<Scalar, Eigen::Lower|Eigen::Upper, Ordering>>;
+
+    std::shared_ptr<SolverBase> createIncompleteCholesky(std::string const& prefix) const
+    {
+      std::string ordering = "amd";
+      Parameters::get(prefix + "->precon->ordering", ordering);
+
+      if (ordering == "amd") {
+        using AMD = Eigen::AMDOrdering<typename Matrix::StorageIndex>;
+        return IncompleteCholesky<AMD>{}.create(prefix);
+      }
+      else {
+        using NATURAL = Eigen::NaturalOrdering<typename Matrix::StorageIndex>;
+        return IncompleteCholesky<NATURAL>{}.create(prefix);
+      }
+    }
   };
 
 
@@ -60,6 +84,9 @@ namespace AMDiS
    * Adds creators for full-matrix aware solvers.
    * - *cg*: conjugate gradient method, \see Eigen::ConjugateGradient
    * - *bcgs*: stabilized bi-conjugate gradient method, \see Eigen::BiCGSTAB
+   * - *minres*: Minimal Residual Algorithm (for symmetric matrices), \see Eigen::MINRES
+   * - *gmres*: Generalized Minimal Residual Algorithm, \see Eigen::GMRES
+   * - *dgmres*: stabilized bi-conjugate gradient method, \see Eigen::DGMRES
    * - *umfpack*: external UMFPACK solver, \see Eigen::UmfPackLU
    * - *superlu*: external SparseLU solver, \see Eigen::SparseLU
    **/
@@ -73,48 +100,64 @@ namespace AMDiS
     using SolverCreator
       = IterativeSolverCreator<Matrix, VectorX, VectorB, IterativeSolver>;
 
+    template <class DirectSolver>
+    using DirectSolverCreator
+      = typename LinearSolver<Matrix, VectorX, VectorB,
+          DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;
+
+    //@{
     template <class Precon>
     using CG = Eigen::ConjugateGradient<Matrix, Eigen::Lower|Eigen::Upper, Precon>;
 
     template <class Precon>
     using BCGS = Eigen::BiCGSTAB<Matrix, Precon>;
 
+    template <class Precon>
+    using MINRES = Eigen::MINRES<Matrix, Eigen::Lower|Eigen::Upper, Precon>;
 
-    template <class DirectSolver>
-    using DirectSolverCreator
-      = typename LinearSolver<Matrix, VectorX, VectorB,
-          DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;
+    template <class Precon>
+    using GMRES = Eigen::GMRES<Matrix, Precon>;
+
+    template <class Precon>
+    using DGMRES = Eigen::DGMRES<Matrix, Precon>;
+    // @}
 
     using Map = CreatorMap<SolverBase>;
 
   public:
     static void init()
     {
+      // conjugate gradient solver
       auto cg = new SolverCreator<CG>;
       Map::addCreator("cg", cg);
 
+      // bi-conjugate gradient stabilized solver
       auto bicgstab = new SolverCreator<BCGS>;
       Map::addCreator("bicgstab", bicgstab);
       Map::addCreator("bcgs", bicgstab);
 
+      // Minimal Residual Algorithm (for symmetric matrices)
+      auto minres = new SolverCreator<MINRES>;
+      Map::addCreator("minres", minres);
+
+      // Generalized Minimal Residual Algorithm
+      auto gmres = new SolverCreator<GMRES>;
+      Map::addCreator("gmres", gmres);
+
+      // Restarted GMRES with deflation.
+      auto dgmres = new SolverCreator<DGMRES>;
+      Map::addCreator("dgmres", dgmres);
+
 #if HAVE_SUITESPARSE_UMFPACK
+      // sparse LU factorization and solver based on UmfPack
       auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU<Matrix>>;
       Map::addCreator("umfpack", umfpack);
 #endif
 
-      using COLAMD = Eigen::COLAMDOrdering<typename Matrix::StorageIndex>;
-      auto sparselu = new DirectSolverCreator<Eigen::SparseLU<Matrix, COLAMD>>;
-      Map::addCreator("sparselu", sparselu);
-      Map::addCreator("sparselu_colamd", sparselu);
-
-      using AMD = Eigen::AMDOrdering<typename Matrix::StorageIndex>;
-      auto sparselu_amd = new DirectSolverCreator<Eigen::SparseLU<Matrix, AMD>>;
-      Map::addCreator("sparselu_amd", sparselu_amd);
-
-#if HAVE_METIS
-      using METIS = Eigen::MetisOrdering<typename Matrix::StorageIndex>;
-      auto sparselu_metis = new DirectSolverCreator<Eigen::SparseLU<Matrix, METIS>>;
-      Map::addCreator("sparselu_metis", sparselu_metis);
+#if HAVE_SUPERLU
+      // sparse direct LU factorization and solver based on the SuperLU library
+      auto superlu = new DirectSolverCreator<Eigen::SuperLU<Matrix>>;
+      Map::addCreator("superlu", superlu);
 #endif
 
       // default iterative solver
@@ -123,8 +166,8 @@ namespace AMDiS
       // default direct solvers
 #if HAVE_SUITESPARSE_UMFPACK
       Map::addCreator("direct", umfpack);
-#else
-      Map::addCreator("direct", sparselu);
+#elif HAVE_SUPERLU
+      Map::addCreator("direct", superlu);
 #endif
     }
   };

src/amdis/linear_algebra/istl/DOFMatrix.hpp

@@ -104,7 +104,6 @@ namespace AMDiS
 
     auto applyDirichletBC(std::vector<char> const& dirichletNodes)
     {
-      Dune::Timer t;
       // loop over the matrix rows
       for (std::size_t i = 0; i < matrix_.N(); ++i) {
         if (dirichletNodes[i]) {
@@ -115,8 +114,6 @@ namespace AMDiS
                 *cIt = (i == cIt.index() ? 1 : 0);
         }
       }
-
-      msg("  time(applyDirichletBC) = {} sec", t.elapsed());
       return std::list<Triplet<value_type>>{};
     }
 

src/amdis/linear_algebra/mtl/DOFMatrix.hpp

@@ -13,7 +13,6 @@
 #include <amdis/Output.hpp>
 #include <amdis/linear_algebra/Common.hpp>
 #include <amdis/linear_algebra/DOFMatrixBase.hpp>
-#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
@@ -85,10 +84,8 @@ namespace AMDiS
     /// final construction of the matrix.
     void finish()
     {
-      Dune::Timer t;
       delete inserter_;
       inserter_ = nullptr;
-      msg("  time(finish insertion) = {} sec", t.elapsed());
     }
 
     /// Return the number of nonzeros in the matrix
@@ -101,8 +98,6 @@ namespace AMDiS
     /// a one on the diagonal of the matrix.
     auto applyDirichletBC(std::vector<char> const& dirichletNodes)
     {
-      Dune::Timer t;
-
       // Define the property maps
       auto row   = mtl::mat::row_map(matrix_);
       auto col   = mtl::mat::col_map(matrix_);
@@ -118,7 +113,6 @@ namespace AMDiS
         }
       }
 
-      msg("  time(applyDirichletBC) = {} sec", t.elapsed());
       return std::list<Triplet<value_type>>{};
     }
 

src/amdis/utility/MultiIndex.hpp

@@ -5,17 +5,17 @@
 
 namespace AMDiS
 {
-  std::size_t flatMultiIndex(std::size_t idx)
+  inline std::size_t flatMultiIndex(std::size_t idx)
   {
     return idx;
   }
 
-  std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx)
+  inline std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx)
   {
     return idx[0];
   }
 
-  std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx)
+  inline std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx)
   {
     return idx[0];
   }