linear algebra backend Eigen3 added

cmake/modules/AmdisMacros.cmake

 include(AMDiSCXXFeatures)
 
-# some additional packages and flags
-find_package(MTL
-  PATHS ${MTL_DIR}
-  HINTS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4
-)
+set(BACKEND "ISTL" CACHE STRING "LinearAlgebra backend. One of MTL, EIGEN, ISTL")
+set_property(CACHE BACKEND PROPERTY STRINGS "MTL" "EIGEN" "ISTL")
 
-if (MTL_FOUND)
-  set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
-  set(MTL_COMPILE_DEFINITIONS "")
-  foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
-    list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
-  endforeach()
+if (BACKEND STREQUAL "MTL")
+  find_package(MTL
+    PATHS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4)
 
-  find_package(SuiteSparse QUIET)
-  if (SuiteSparse_FOUND)
-    list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
-  endif (SuiteSparse_FOUND)
-endif (MTL_FOUND)
+  if (MTL_FOUND)
+    set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
+    set(MTL_COMPILE_DEFINITIONS "")
+    foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
+      list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
+    endforeach()
 
-set(HAVE_MTL MTL_FOUND)
+    find_package(SuiteSparse QUIET)
+    if (SuiteSparse_FOUND)
+      list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
+    endif (SuiteSparse_FOUND)
+  endif (MTL_FOUND)
 
-if (MTL_FOUND)
-  list(APPEND MTL_COMPILE_DEFINITIONS "ENABLE_MTL=1")
-  dune_register_package_flags(
-    COMPILE_DEFINITIONS ${MTL_COMPILE_DEFINITIONS}
-    INCLUDE_DIRS ${MTL_INCLUDE_DIRS})
-endif (MTL_FOUND)
\ No newline at end of file
+  set(HAVE_MTL MTL_FOUND)
+  if (MTL_FOUND)
+    list(APPEND MTL_COMPILE_DEFINITIONS "ENABLE_MTL=1")
+    dune_register_package_flags(
+      COMPILE_DEFINITIONS ${MTL_COMPILE_DEFINITIONS}
+      INCLUDE_DIRS ${MTL_INCLUDE_DIRS})
+  endif (MTL_FOUND)
+elseif (BACKEND STREQUAL "EIGEN")
+  find_package(Eigen3)
+
+  set(HAVE_EIGEN EIGEN_FOUND)
+  if (EIGEN3_FOUND)
+    message(STATUS "  Found Eigen3, version: ${EIGEN3_VERSION}")
+    list(APPEND EIGEN3_DEFINITIONS "ENABLE_EIGEN=1")
+    dune_register_package_flags(
+      COMPILE_DEFINITIONS ${EIGEN3_DEFINITIONS}
+      INCLUDE_DIRS ${EIGEN3_INCLUDE_DIRS})
+  endif (EIGEN3_FOUND)
+endif ()

config.h.cmake

@@ -40,9 +40,12 @@
 /* Define to the revision of amdis */
 #define AMDIS_VERSION_REVISION @AMDIS_VERSION_REVISION@
 
-/* Define to ENABLE_UMFPACK if the MTL library is available */
+/* Define to ENABLE_MTL if the MTL library is available */
 #cmakedefine HAVE_MTL ENABLE_MTL
 
+/* Define to ENABLE_EIGEN if the Eigen3 library is available */
+#cmakedefine HAVE_EIGEN ENABLE_EIGEN
+
 /* some detected compiler features may be used in AMDiS */
 #cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
 #cmakedefine AMDIS_HAS_CXX_CONSTEXPR_IF 1

examples/ellipt.cc

@@ -78,9 +78,9 @@ int main(int argc, char** argv)
   }
 
   std::cout << "\n";
-  msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n",
+  msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
       "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
-  msg("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
+  msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
 
   std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
   for (int i = 1; i < numLevels; ++i) {
@@ -89,7 +89,7 @@ int main(int argc, char** argv)
   }
 
   for (int i = 0; i < numLevels; ++i)
-    msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}\n",
+    msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
         i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
 
   AMDiS::finalize();

examples/init/ellipt.dat.2d

@@ -5,9 +5,10 @@ ellipt->mesh:                   elliptMesh
 ellipt->solver->name:           bicgstab
 ellipt->solver->max iteration:  10000
 ellipt->solver->relative tolerance: 1.e-10
-ellipt->solver->info:           1
+ellipt->solver->info:           10
 ellipt->solver->left precon:    ilu
 ellipt->solver->right precon:   no
+ellipt->solver->precon->name:   ilu
 
 ellipt->output[0]->filename:    ellipt.2d
 ellipt->output[0]->name:        u

src/amdis/LinearAlgebra.hpp

@@ -4,13 +4,23 @@
 #include <amdis/linear_algebra/SolverInfo.hpp>
 
 #if HAVE_MTL
-#include <amdis/linear_algebra/mtl/DOFVector.hpp>
+
 #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>
-#else
-#include <amdis/linear_algebra/istl/DOFVector.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/linear_algebra/CMakeLists.txt

-#install headers
-
 install(FILES
     Common.hpp
     DOFMatrixBase.hpp
@@ -12,4 +10,6 @@ install(FILES
     SolverInfo.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra)
 
+add_subdirectory("eigen")
+add_subdirectory("istl")
 add_subdirectory("mtl")

src/amdis/linear_algebra/DOFMatrixBase.hpp

 #pragma once
 
 #include <cmath>
+
+#include <dune/common/timer.hh>
+
 #include <amdis/common/Math.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
@@ -67,6 +71,7 @@ namespace AMDiS
     /// If \p setToZero is true, the matrix is set to 0
     void init(bool setToZero)
     {
+      timeInsertion_ = 0;
       backend_.init(*rowBasis_, *colBasis_, setToZero);
     }
 
@@ -74,6 +79,7 @@ namespace AMDiS
     void finish()
     {
       backend_.finish();
+      msg("  time(insertion) = {} sec", timeInsertion_);
     }
 
     /// Insert a block of values into the matrix (add to existing values)
@@ -82,6 +88,7 @@ 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) {
@@ -91,6 +98,7 @@ namespace AMDiS
           }
         }
       }
+      timeInsertion_ += t.elapsed();
     }
 
     /// Insert a single value into the matrix (add to existing value)
@@ -117,6 +125,8 @@ namespace AMDiS
     ColBasis const*	colBasis_;
 
     Backend backend_;
+
+    double timeInsertion_ = 0;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/DOFVectorBase.hpp

@@ -59,7 +59,7 @@ namespace AMDiS
       std::enable_if_t<Concepts::Arithmetic<Scalar>, int> = 0>
     Self& operator=(Scalar value)
     {
-      vector() = value;
+      backend_.set(value);
       return *this;
     }
 
@@ -95,9 +95,9 @@ namespace AMDiS
     /// Resize the \ref vector to the size of the \ref basis and set to zero
     virtual void compress() override
     {
-      if (backend_.size() != size()) {
+      if (size_type(backend_.size()) != size()) {
         backend_.resize(size());
-        vector() = 0;
+        backend_.set(0);
       }
     }
 

src/amdis/linear_algebra/eigen/CMakeLists.txt

+install(FILES
+    DirectRunner.hpp
+    DOFMatrix.hpp
+    DOFVector.hpp
+    IterativeRunner.hpp
+    PreconConfig.hpp
+    SolverConfig.hpp
+    SolverCreator.hpp
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/eigen)

src/amdis/linear_algebra/eigen/DOFMatrix.hpp

+#pragma once
+
+#include <list>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <Eigen/SparseCore>
+
+#include <dune/common/timer.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+
+namespace AMDiS
+{
+  /// \brief The basic container that stores a base matrix and a corresponding
+  /// row/column feSpace.
+  template <class ValueType, int Orientation = Eigen::RowMajor>
+  class EigenMatrix
+  {
+  public:
+    /// The matrix type of the underlying base matrix
+    using BaseMatrix = Eigen::SparseMatrix<ValueType, Orientation>;
+
+    /// The type of the elements of the DOFMatrix
+    using value_type = ValueType;
+
+    /// The index/size - type
+    using size_type = typename BaseMatrix::Index;
+
+  public:
+    /// Constructor. Constructs new BaseMatrix.
+    EigenMatrix() = default;
+
+    /// Return a reference to the data-matrix \ref matrix
+    BaseMatrix& matrix()
+    {
+      return matrix_;
+    }
+
+    /// Return a reference to the data-matrix \ref matrix
+    BaseMatrix const& matrix() const
+    {
+      return matrix_;
+    }
+
+
+    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
+    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
+    /// be created by \ref init and must be closed by \ref finish after insertion.
+    void insert(size_type r, size_type c, value_type const& value)
+    {
+      test_exit_dbg(r < matrix_.rows() && c < matrix_.cols(),
+          "Indices out of range [0,{})x[0,{})", matrix_.rows(), matrix_.cols());
+      tripletList_.emplace_back(r, c, value);
+    }
+
+
+    /// Create inserter. Assumes that no inserter is currently active on this matrix.
+    template <class RowBasis, class ColBasis>
+    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
+      }
+    }
+
+    /// Delete inserter -> finish insertion. Must be called in order to fill the
+    /// 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
+    }
+
+    /// Return the number of nonzeros in the matrix
+    std::size_t nnz() const
+    {
+      return matrix_.nonZeros();
+    }
+
+    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
+    /// 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);
+      }
+      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)
+            it.valueRef() = (it.row() == it.col() ? value_type(1) : value_type(0));
+        }
+      }
+    }
+
+    void clearDirichletRow(size_type idx, int_t<Eigen::RowMajor>)
+    {
+      Eigen::SparseVector<ValueType> unit_vector(matrix_.cols());
+      unit_vector.coeffRef(idx) = value_type(1);
+      matrix_.row(idx) = unit_vector;
+    }
+
+  private:
+    /// The data-matrix
+    BaseMatrix matrix_;
+
+    /// A list of row,col indices and values
+    std::vector<Eigen::Triplet<ValueType, Eigen::Index>> tripletList_;
+  };
+
+  template <class RowBasisType, class ColBasisType,
+            class T = double, int O = Eigen::RowMajor>
+  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, EigenMatrix<T, O>>;
+
+} // end namespace AMDiS

src/amdis/linear_algebra/eigen/DOFVector.hpp

+#pragma once
+
+#include <Eigen/Dense>
+
+#include <dune/common/ftraits.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
+
+namespace AMDiS
+{
+  /// The basic container that stores a base vector and a corresponding basis
+  template <class ValueType>
+  class EigenVector
+  {
+  public:
+    /// The type of the elements of the DOFVector
+    using value_type = ValueType;
+
+    /// The type of the elements of the DOFVector
+    using block_type = ValueType;
+
+    /// The underlying field type
+    using field_type = typename Dune::FieldTraits<ValueType>::field_type;
+
+    /// The type of the base vector
+    using BaseVector = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>;
+
+    /// The index/size - type
+    using size_type  = typename BaseVector::Index;
+
+  public:
+    /// Constructor. Constructs new BaseVector.
+    EigenVector() = default;
+
+    /// Return the data-vector \ref vector_
+    BaseVector const& vector() const
+    {
+      return vector_;
+    }
+
+    /// Return the data-vector \ref vector_
+    BaseVector& vector()
+    {
+      return vector_;
+    }
+
+    /// Return the current size of the \ref vector_
+    size_type size() const
+    {
+      return vector_.size();
+    }
+
+    /// Resize the \ref vector_ to the size \p s
+    void resize(size_type s)
+    {
+      vector_.resize(s);
+    }
+
+
+    /// Access the entry \p i of the \ref vector with read-access.
+    value_type const& operator[](size_type i) const
+    {
+      test_exit_dbg(i < size(),
+        "Index {} out of range [0,{})", i, size());
+      return vector_.coeff(i);
+    }
+
+    /// Access the entry \p i of the \ref vector with write-access.
+    value_type& operator[](size_type i)
+    {
+      test_exit_dbg(i < size(),
+        "Index {} out of range [0,{})", i, size());
+      return vector_.coeffRef(i);
+    }
+
+    void set(field_type value)
+    {
+      vector_.setConstant(value);
+    }
+
+  private:
+    /// The data-vector (can hold a new BaseVector or a pointer to external data
+    BaseVector vector_;
+  };
+
+
+  template <class BasisType, class ValueType = double>
+  struct DOFVector : public DOFVectorBase<BasisType, EigenVector<ValueType>>
+  {
+    using Super = DOFVectorBase<BasisType, EigenVector<ValueType>>;
+    using Super::operator=;
+    DOFVector(BasisType const& basis)
+      : Super(basis)
+    {}
+  };
+
+
+  /// Constructor a dofvector from given basis and name
+  template <class ValueType = double, class Basis>
+  DOFVector<Basis, ValueType>
+  makeDOFVector(Basis const& basis)
+  {
+    return {basis};
+  }
+
+} // end namespace AMDiS

src/amdis/linear_algebra/eigen/DirectRunner.hpp

+#pragma once
+
+#include <algorithm>
+#include <string>
+
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+#include <amdis/linear_algebra/eigen/SolverConfig.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \ingroup Solver
+   * \class AMDiS::DirectRunner
+   * \brief \implements RunnerInterface for the (external) direct solvers
+   */
+  template <class LU, class Matrix, class VectorX, class VectorB>
+  class DirectRunner
+      : public RunnerInterface<Matrix, VectorX, VectorB>
+  {
+  protected:
+    using Super = RunnerInterface<Matrix, VectorX, VectorB>;
+
+  public:
+    /// Constructor.
+    DirectRunner(std::string prefix)
+      : solver_{}
+    {
+      SolverConfig<LU>::init(prefix, solver_);
+    }
+
+    /// Implementes \ref RunnerInterface::init()
+    virtual void init(Matrix const& A) override
+    {
+      solver_.compute(A);
+      test_exit(solver_.info() == Eigen::Success, "Error in solver.compute(matrix)");
+    }
+
+    /// Implementes \ref RunnerInterface::exit()
+    virtual void exit() override {}
+
+    /// Implementes \ref RunnerInterface::solve()
+    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
+                      SolverInfo& solverInfo) override
+    {
+      x = solver_.solve(b);
+
+      auto r = VectorB(b);
+      if (x.norm() != 0)
+        r -= A * x;
+
+      solverInfo.setAbsResidual(r.norm());
+      solverInfo.setError(solver_.info());
+
+      return solver_.info() == Eigen::Success ? 0 : 1;
+    }
+
+  private:
+    LU solver_;
+  };
+}

src/amdis/linear_algebra/eigen/IterativeRunner.hpp

+#pragma once
+
+#include <dune/istl/preconditioner.hh>
+
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+#include <amdis/linear_algebra/eigen/PreconConfig.hpp>
+#include <amdis/linear_algebra/eigen/SolverConfig.hpp>
+
+namespace AMDiS
+{
+  template <class IterativeSolver, class Matrix, class VectorX, class VectorB>
+  class IterativeRunner
+      : public RunnerInterface<Matrix, VectorX, VectorB>
+  {
+    using SolverCfg = SolverConfig<IterativeSolver>;
+    using PreconCfg = PreconConfig<typename IterativeSolver::Preconditioner>;
+
+  public:
+    IterativeRunner(std::string prefix)
+      : solver_{}
+    {
+      SolverCfg::init(prefix, solver_);
+      PreconCfg::init(prefix + "->precon", solver_.preconditioner());
+    }
+
+
+    /// Implementes \ref RunnerInterface::init()
+    virtual void init(Matrix const& A) override
+    {
+      solver_.compute(A);
+      test_exit(solver_.info() == Eigen::Success, "Error in solver.compute(matrix)");
+    }
+
+    /// Implementes \ref RunnerInterface::exit()
+    virtual void exit() override {}
+
+    /// Implementes \ref RunnerInterface::solve()
+    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
+                      SolverInfo& solverInfo) override
+    {
+      solver_.setTolerance(solverInfo.getRelTolerance());
+      x = solver_.solveWithGuess(b, x);
+
+      auto r = VectorB(b);
+      if (x.norm() != 0)
+        r -= A * x;
+
+      solverInfo.setAbsResidual(r.norm());
+      solverInfo.setRelResidual(solver_.error());
+      solverInfo.setError(solver_.info());
+      msg("number of iteration: {}", solver_.iterations());
+
+      return solver_.info() == Eigen::Success ? 0 : int(solver_.info());
+    }
+
+  private:
+    IterativeSolver solver_;
+  };