From 4bcf59c2c5c8b12b13c93537abccd5c447670f4c Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Thu, 9 Jul 2015 11:09:51 +0200
Subject: [PATCH] [Mat] add LinAlg/Eigen and LinAlg/EigenLis
---
MathLib/CMakeLists.txt | 10 ++
MathLib/LinAlg/ApplyKnownSolution.h | 4 +
MathLib/LinAlg/Eigen/EigenLinearSolver.cpp | 145 +++++++++++++++
MathLib/LinAlg/Eigen/EigenLinearSolver.h | 80 +++++++++
MathLib/LinAlg/Eigen/EigenMatrix.h | 170 ++++++++++++++++++
MathLib/LinAlg/Eigen/EigenOption.cpp | 48 +++++
MathLib/LinAlg/Eigen/EigenOption.h | 82 +++++++++
MathLib/LinAlg/Eigen/EigenTools.cpp | 58 ++++++
MathLib/LinAlg/Eigen/EigenTools.h | 36 ++++
MathLib/LinAlg/Eigen/EigenVector.h | 114 ++++++++++++
.../LinAlg/EigenLis/EigenLisLinearSolver.cpp | 82 +++++++++
.../LinAlg/EigenLis/EigenLisLinearSolver.h | 76 ++++++++
Tests/MathLib/TestGlobalMatrixInterface.cpp | 12 ++
Tests/MathLib/TestGlobalVectorInterface.cpp | 11 ++
Tests/MathLib/TestLinearSolver.cpp | 44 +++++
15 files changed, 972 insertions(+)
create mode 100644 MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
create mode 100644 MathLib/LinAlg/Eigen/EigenLinearSolver.h
create mode 100644 MathLib/LinAlg/Eigen/EigenMatrix.h
create mode 100644 MathLib/LinAlg/Eigen/EigenOption.cpp
create mode 100644 MathLib/LinAlg/Eigen/EigenOption.h
create mode 100644 MathLib/LinAlg/Eigen/EigenTools.cpp
create mode 100644 MathLib/LinAlg/Eigen/EigenTools.h
create mode 100644 MathLib/LinAlg/Eigen/EigenVector.h
create mode 100644 MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
create mode 100644 MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 4887c539178..09c09926eac 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -23,11 +23,21 @@ set(SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
+if(OGS_USE_EIGEN)
+ GET_SOURCE_FILES(SOURCES_LINALG_EIGEN LinAlg/Eigen)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGEN})
+endif()
+
if(OGS_USE_LIS)
GET_SOURCE_FILES(SOURCES_LINALG_LIS LinAlg/Lis)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_LIS})
endif()
+if(OGS_USE_EIGEN AND OGS_USE_LIS)
+ GET_SOURCE_FILES(SOURCES_LINALG_EIGENLIS LinAlg/EigenLis)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGENLIS})
+endif()
+
if(OGS_USE_PETSC)
GET_SOURCE_FILES(SOURCES_LINALG_PETSC LinAlg/PETSc)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC})
diff --git a/MathLib/LinAlg/ApplyKnownSolution.h b/MathLib/LinAlg/ApplyKnownSolution.h
index 31c99e56cc2..72c28fe3481 100644
--- a/MathLib/LinAlg/ApplyKnownSolution.h
+++ b/MathLib/LinAlg/ApplyKnownSolution.h
@@ -13,6 +13,10 @@
#include "MathLib/LinAlg/Dense/DenseTools.h"
+#ifdef OGS_USE_EIGEN
+#include "MathLib/LinAlg/Eigen/EigenTools.h"
+#endif // USE_LIS
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisTools.h"
#endif // USE_LIS
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
new file mode 100644
index 00000000000..8071c650e30
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -0,0 +1,145 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EigenLinearSolver.h"
+
+#include <logog/include/logog.hpp>
+#include <unsupported/Eigen/IterativeSolvers>
+
+#include "EigenVector.h"
+#include "EigenMatrix.h"
+#include "EigenTools.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+namespace details
+{
+
+template <class T_SOLVER>
+class EigenDirectSolver : public EigenLinearSolver::IEigenSolver
+{
+public:
+ EigenDirectSolver(EigenMatrix::RawMatrixType &A) : _A(A)
+ {
+ // fill A and b;
+ // Compute the ordering permutation vector from the structural pattern of A
+ _solver.analyzePattern(A);
+ // Compute the numerical factorization
+ _solver.factorize(A);
+ if(_solver.info()!=Eigen::Success) {
+ INFO("\t failed");
+ return;
+ }
+ }
+
+ virtual ~EigenDirectSolver() {}
+
+ void solve(EigenVector::RawVectorType &b, EigenVector::RawVectorType &x, EigenOption &) override
+ {
+ //Use the factors to solve the linear system
+ INFO("-> solve");
+ x = _solver.solve(b);
+ }
+
+private:
+ T_SOLVER _solver;
+ EigenMatrix::RawMatrixType& _A;
+};
+
+template <class T_SOLVER>
+class EigenIterativeSolver : public EigenLinearSolver::IEigenSolver
+{
+public:
+ EigenIterativeSolver(EigenMatrix::RawMatrixType &A) : _A(A)
+ {
+ INFO("-> initialize with the coefficient matrix");
+ _solver.compute(A);
+ if(_solver.info()!=Eigen::Success) {
+ INFO("\t failed");
+ return;
+ }
+ }
+
+ virtual ~EigenIterativeSolver() {}
+
+ void solve(EigenVector::RawVectorType &b, EigenVector::RawVectorType &x, EigenOption &opt) override
+ {
+ INFO("-> solve");
+ _solver.setTolerance(opt.error_tolerance);
+ _solver.setMaxIterations(opt.max_iterations);
+ x = _solver.solveWithGuess(b, x);
+ if(_solver.info()!=Eigen::Success) {
+ INFO("\t solving failed");
+ return;
+ }
+ INFO("\t iteration: %d/%ld", _solver.iterations(), opt.max_iterations);
+ INFO("\t residual: %e\n", _solver.error());
+ }
+
+private:
+ T_SOLVER _solver;
+ EigenMatrix::RawMatrixType& _A;
+};
+
+} // details
+
+EigenLinearSolver::EigenLinearSolver(EigenMatrix &A, ptree const*const option)
+{
+ if (option)
+ setOption(*option);
+
+ A.getRawMatrix().makeCompressed();
+ if (_option.solver_type==EigenOption::SolverType::SparseLU) {
+ typedef Eigen::SparseLU<EigenMatrix::RawMatrixType, Eigen::COLAMDOrdering<int> > SolverType;
+ _solver = new details::EigenDirectSolver<SolverType>(A.getRawMatrix());
+ } else if (_option.solver_type==EigenOption::SolverType::BiCGSTAB) {
+ typedef Eigen::BiCGSTAB<EigenMatrix::RawMatrixType, Eigen::DiagonalPreconditioner<double>> SolverType;
+ _solver = new details::EigenIterativeSolver<SolverType>(A.getRawMatrix());
+ } else if (_option.solver_type==EigenOption::SolverType::CG) {
+ typedef Eigen::ConjugateGradient<EigenMatrix::RawMatrixType, Eigen::Lower, Eigen::DiagonalPreconditioner<double>> SolverType;
+ _solver = new details::EigenIterativeSolver<SolverType>(A.getRawMatrix());
+ }
+}
+
+void EigenLinearSolver::setOption(const ptree &option)
+{
+ boost::optional<ptree> ptSolver = option.get_child("LinearSolver");
+ if (!ptSolver)
+ return;
+
+ boost::optional<std::string> solver_type = ptSolver->get_optional<std::string>("solver_type");
+ if (solver_type) {
+ _option.solver_type = _option.getSolverType(*solver_type);
+ }
+ boost::optional<std::string> precon_type = ptSolver->get_optional<std::string>("precon_type");
+ if (precon_type) {
+ _option.precon_type = _option.getPreconType(*precon_type);
+ }
+ boost::optional<double> error_tolerance = ptSolver->get_optional<double>("error_tolerance");
+ if (error_tolerance) {
+ _option.error_tolerance = *error_tolerance;
+ }
+ boost::optional<int> max_iteration_step = ptSolver->get_optional<int>("max_iteration_step");
+ if (max_iteration_step) {
+ _option.max_iterations = *max_iteration_step;
+ }
+}
+
+void EigenLinearSolver::solve(EigenVector &b, EigenVector &x)
+{
+ INFO("------------------------------------------------------------------");
+ INFO("*** Eigen solver computation");
+ _solver->solve(b.getRawVector(), x.getRawVector(), _option);
+ INFO("------------------------------------------------------------------");
+}
+
+} //MathLib
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
new file mode 100644
index 00000000000..d536b157c9a
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
@@ -0,0 +1,80 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENLINEARSOLVER_H_
+#define EIGENLINEARSOLVER_H_
+
+#include <vector>
+
+#include <boost/any.hpp>
+#include <boost/property_tree/ptree.hpp>
+
+#include "EigenVector.h"
+#include "EigenOption.h"
+
+namespace MathLib
+{
+
+class EigenMatrix;
+
+/**
+ */
+class EigenLinearSolver
+{
+public:
+ EigenLinearSolver(EigenMatrix &A, boost::property_tree::ptree const*const option = nullptr);
+
+ virtual ~EigenLinearSolver()
+ {
+ delete _solver;
+ }
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const boost::property_tree::ptree &option);
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const EigenOption &option) { _option = option; }
+
+ /**
+ * get linear solver options
+ * @return
+ */
+ EigenOption &getOption() { return _option; }
+
+ /**
+ * solve a given linear equations
+ *
+ * @param b RHS vector
+ * @param x Solution vector
+ */
+ void solve(EigenVector &b, EigenVector &x);
+
+protected:
+ class IEigenSolver
+ {
+ public:
+ virtual ~IEigenSolver() {}
+ /// execute a linear solver
+ virtual void solve(EigenVector::RawVectorType &b, EigenVector::RawVectorType &x, EigenOption &) = 0;
+ };
+
+ EigenOption _option;
+ IEigenSolver* _solver;
+};
+
+} // MathLib
+
+#endif //EIGENLINEARSOLVER_H_
+
diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
new file mode 100644
index 00000000000..fb254a58e0f
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenMatrix.h
@@ -0,0 +1,170 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENMATRIX_H_
+#define EIGENMATRIX_H_
+
+#include <cassert>
+#include <fstream>
+#include <iosfwd>
+#include <string>
+
+#include <Eigen/Sparse>
+
+#include "MathLib/LinAlg/RowColumnIndices.h"
+#include "EigenVector.h"
+
+namespace MathLib
+{
+
+/**
+ * Global matrix based on Eigen sparse matrix
+ */
+class EigenMatrix
+{
+public:
+ typedef Eigen::SparseMatrix<double, Eigen::RowMajor> RawMatrixType;
+
+ /**
+ * constructor
+ * @param n_rows the number of rows (that is equal to the number of columns)
+ */
+ explicit EigenMatrix(std::size_t n) :_mat(n, n) {}
+
+ virtual ~EigenMatrix() {}
+
+ /// return the number of rows
+ virtual std::size_t getNRows() const { return _mat.rows(); }
+
+ /// return the number of columns
+ virtual std::size_t getNCols() const { return _mat.cols(); }
+
+ /// return a start index of the active data range
+ virtual std::size_t getRangeBegin() const { return 0; }
+
+ /// return an end index of the active data range
+ virtual std::size_t getRangeEnd() const { return getNRows(); }
+
+ /// Reset data entries to zero.
+ virtual void setZero()
+ {
+ _mat.setZero();
+ }
+
+ /// set entry
+ int setValue(std::size_t row, std::size_t col, double val)
+ {
+ _mat.coeffRef(row, col) = val;
+ return 0;
+ }
+
+ /// add value
+ int add(std::size_t row, std::size_t col, double val)
+ {
+ assert(row < getNRows() && col < getNCols());
+ _mat.coeffRef(row, col) += val;
+ return 0;
+ }
+
+ /// Add sub-matrix at positions \c row_pos and same column positions as the
+ /// given row positions.
+ template<class T_DENSE_MATRIX>
+ void add(std::vector<std::size_t> const& row_pos,
+ const T_DENSE_MATRIX &sub_matrix,
+ double fkt = 1.0)
+ {
+ this->add(row_pos, row_pos, sub_matrix, fkt);
+ }
+
+ /// Add sub-matrix at positions given by \c indices.
+ template<class T_DENSE_MATRIX>
+ void add(RowColumnIndices<std::size_t> const& indices,
+ const T_DENSE_MATRIX &sub_matrix,
+ double fkt = 1.0)
+ {
+ this->add(indices.rows, indices.columns, sub_matrix, fkt);
+ }
+
+ ///
+ template <class T_DENSE_MATRIX>
+ void add(std::vector<std::size_t> const& row_pos,
+ std::vector<std::size_t> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
+ double fkt = 1.0);
+
+ /// get value
+ double get(std::size_t row, std::size_t col)
+ {
+ assert(row < getNRows() && col < getNCols());
+ return _mat.coeff(row, col);
+ }
+
+ /// get value
+ double operator() (std::size_t row, std::size_t col) const
+ {
+ return _mat.coeff(row, col);
+ }
+
+ /// get a maximum value in diagonal entries
+ virtual double getMaxDiagCoeff() {return 1;}
+
+ /// y = mat * x
+ virtual void multiply(const EigenVector &x, EigenVector &y) const
+ {
+ y.getRawVector() = _mat * x.getRawVector();
+ }
+
+ /// return if this matrix is already assembled or not
+ bool isAssembled() const { return true; }
+
+ /// printout this matrix for debugging
+ void write(const std::string &filename) const
+ {
+ std::ofstream of(filename);
+ if (of)
+ write(of);
+ }
+
+ /// printout this matrix for debugging
+ void write(std::ostream &os) const
+ {
+ for (int k=0; k<_mat.outerSize(); ++k)
+ for (Eigen::SparseMatrix<double>::InnerIterator it(_mat,k); it; ++it)
+ os << it.row() << " " << it.col() << ": " << it.value() << "\n";
+ os << std::endl;
+ }
+
+ RawMatrixType& getRawMatrix() { return _mat; }
+ const RawMatrixType& getRawMatrix() const { return _mat; }
+
+protected:
+ RawMatrixType _mat;
+};
+
+
+template<class T_DENSE_MATRIX>
+void
+EigenMatrix::add(std::vector<std::size_t> const& row_pos, std::vector<std::size_t> const& col_pos,
+ const T_DENSE_MATRIX &sub_matrix, double fkt)
+{
+ const std::size_t n_rows = row_pos.size();
+ const std::size_t n_cols = col_pos.size();
+ for (std::size_t i = 0; i < n_rows; i++) {
+ const std::size_t row = row_pos[i];
+ for (std::size_t j = 0; j < n_cols; j++) {
+ const std::size_t col = col_pos[j];
+ add(row, col, fkt * sub_matrix(i, j));
+ }
+ }
+};
+
+
+} // end namespace MathLib
+
+#endif
+
diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
new file mode 100644
index 00000000000..03533f688a3
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -0,0 +1,48 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EigenOption.h"
+
+namespace MathLib
+{
+
+EigenOption::EigenOption()
+{
+ solver_type = SolverType::SparseLU;
+ precon_type = PreconType::NONE;
+ max_iterations = 1e6;
+ error_tolerance = 1.e-16;
+}
+
+EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_name)
+{
+#define RETURN_SOLVER_ENUM_IF_SAME_STRING(str, TypeName) \
+ if (#TypeName==str) return SolverType::TypeName;
+
+ RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, CG);
+ RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCGSTAB);
+ RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, SparseLU);
+
+ return SolverType::INVALID;
+#undef RETURN_SOLVER_ENUM_IF_SAME_STRING
+}
+
+EigenOption::PreconType EigenOption::getPreconType(const std::string &precon_name)
+{
+#define RETURN_PRECOM_ENUM_IF_SAME_STRING(str, TypeName) \
+ if (#TypeName==str) return PreconType::TypeName;
+
+ RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, NONE);
+ RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, DIAGONAL);
+
+ return PreconType::NONE;
+#undef RETURN_PRECOM_ENUM_IF_SAME_STRING
+}
+
+} //MathLib
diff --git a/MathLib/LinAlg/Eigen/EigenOption.h b/MathLib/LinAlg/Eigen/EigenOption.h
new file mode 100644
index 00000000000..ba911c49a9d
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenOption.h
@@ -0,0 +1,82 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENOPTION_H_
+#define EIGENOPTION_H_
+
+#include <string>
+
+namespace MathLib
+{
+
+/**
+ * \brief Option for Eigen sparse solver
+ */
+struct EigenOption
+{
+ /// Solver type
+ enum class SolverType : int
+ {
+ INVALID,
+ CG,
+ BiCGSTAB,
+ SparseLU
+ };
+
+ /// Preconditioner type
+ enum class PreconType : int
+ {
+ NONE,
+ DIAGONAL
+ };
+
+ /// Linear solver type
+ SolverType solver_type;
+ /// Preconditioner type
+ PreconType precon_type;
+ /// Maximum iteration count
+ long max_iterations;
+ /// Error tolerance
+ double error_tolerance;
+
+ /**
+ * Constructor
+ *
+ * Default options are CG, no preconditioner, iteration count 500 and
+ * tolerance 1e-10. Default matrix storage type is CRS.
+ */
+ EigenOption();
+
+ /// Destructor
+ ~EigenOption() {}
+
+ /**
+ * return a linear solver type from the solver name
+ *
+ * @param solver_name
+ * @return a linear solver type
+ * If there is no solver type matched with the given name, INVALID
+ * is returned.
+ */
+ static SolverType getSolverType(const std::string &solver_name);
+
+ /**
+ * return a preconditioner type from the name
+ *
+ * @param precon_name
+ * @return a preconditioner type
+ * If there is no preconditioner type matched with the given name, NONE
+ * is returned.
+ */
+ static PreconType getPreconType(const std::string &precon_name);
+
+};
+
+}
+#endif //LIS_OPTION_H_
diff --git a/MathLib/LinAlg/Eigen/EigenTools.cpp b/MathLib/LinAlg/Eigen/EigenTools.cpp
new file mode 100644
index 00000000000..b820e768c84
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenTools.cpp
@@ -0,0 +1,58 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EigenTools.h"
+
+#include <logog/include/logog.hpp>
+
+#include "EigenMatrix.h"
+#include "EigenVector.h"
+
+namespace MathLib
+{
+
+void applyKnownSolution(EigenMatrix &A_, EigenVector &b_, const std::vector<std::size_t> &vec_knownX_id,
+ const std::vector<double> &vec_knownX_x, double /*penalty_scaling*/)
+{
+ typedef EigenMatrix::RawMatrixType SpMat;
+ auto &A = A_.getRawMatrix();
+ auto &b = b_.getRawVector();
+ const size_t n_rows = A.rows();
+ for (std::size_t ix=0; ix<vec_knownX_id.size(); ix++)
+ {
+ int row_id = vec_knownX_id[ix];
+ auto x = vec_knownX_x[ix];
+ //A(k, j) = 0.
+ for (SpMat::InnerIterator it(A,row_id); it; ++it)
+ it.valueRef() = .0;
+ //b_i -= A(i,k)*val, i!=k
+ for (size_t i=0; i<n_rows; i++)
+ for (SpMat::InnerIterator it(A,i); it; ++it)
+ {
+ if (it.col()!=row_id) continue;
+ b[i] -= it.value()*x;
+ }
+ //b_k = val
+ b[row_id] = x;
+ //A(i, k) = 0., i!=k
+ for (size_t i=0; i<n_rows; i++)
+ for (SpMat::InnerIterator it(A,i); it; ++it)
+ {
+ if (it.col()!=row_id) continue;
+ it.valueRef() = 0.0;
+ }
+ //A(k, k) = 1.0
+ A.coeffRef(row_id, row_id) = 1.0; //=x
+ }
+}
+
+} // MathLib
+
+
+
diff --git a/MathLib/LinAlg/Eigen/EigenTools.h b/MathLib/LinAlg/Eigen/EigenTools.h
new file mode 100644
index 00000000000..95995e0c259
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenTools.h
@@ -0,0 +1,36 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENTOOLS_H_
+#define EIGENTOOLS_H_
+
+#include <vector>
+
+namespace MathLib
+{
+class EigenMatrix;
+class EigenVector;
+
+/**
+ * apply known solutions to a system of linear equations
+ *
+ * @param A Coefficient matrix
+ * @param b RHS vector
+ * @param _vec_knownX_id a vector of known solution entry IDs
+ * @param _vec_knownX_x a vector of known solutions
+ * @param penalty_scaling value for scaling some matrix and right hand side
+ * entries to enforce some conditions
+ */
+void applyKnownSolution(EigenMatrix &A, EigenVector &b, const std::vector<std::size_t> &_vec_knownX_id,
+ const std::vector<double> &_vec_knownX_x, double penalty_scaling = 1e+10);
+
+} // MathLib
+
+#endif //EIGENTOOLS_H_
+
diff --git a/MathLib/LinAlg/Eigen/EigenVector.h b/MathLib/LinAlg/Eigen/EigenVector.h
new file mode 100644
index 00000000000..edc93f4af2e
--- /dev/null
+++ b/MathLib/LinAlg/Eigen/EigenVector.h
@@ -0,0 +1,114 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENVECTOR_H_
+#define EIGENVECTOR_H_
+
+#include <string>
+#include <vector>
+#include <fstream>
+
+#include <Eigen/Eigen>
+
+namespace MathLib
+{
+
+/**
+ * Global vector based on Eigen vector
+ */
+class EigenVector
+{
+public:
+ typedef Eigen::VectorXd RawVectorType;
+
+ /**
+ * Constructor for initialization of the number of rows
+ * @param length number of rows
+ */
+ explicit EigenVector(std::size_t length) : _vec(length) {}
+
+ /// copy constructor
+ EigenVector(EigenVector const &src) : _vec(src._vec) {}
+
+ /**
+ *
+ */
+ virtual ~EigenVector() {}
+
+ /// return a vector length
+ std::size_t size() const { return _vec.size(); }
+
+ /// return a start index of the active data range
+ std::size_t getRangeBegin() const { return 0;}
+
+ /// return an end index of the active data range
+ std::size_t getRangeEnd() const { return this->size(); }
+
+ /// set all values in this vector
+ EigenVector& operator= (double v) { _vec.setConstant(v); return *this; }
+
+ /// set all values in this vector
+ EigenVector& operator*= (double v) { _vec *= v; return *this; }
+
+ /// access entry
+ double operator[] (std::size_t rowId) const { return get(rowId); }
+
+ /// get entry
+ double get(std::size_t rowId) const
+ {
+ return _vec[rowId];
+ }
+
+ /// set entry
+ void set(std::size_t rowId, double v)
+ {
+ _vec[rowId] = v;
+ }
+
+ /// add entry
+ void add(std::size_t rowId, double v)
+ {
+ _vec[rowId] += v;
+ }
+
+ /// add entries
+ template<class T_SUBVEC>
+ void add(const std::vector<std::size_t> &pos, const T_SUBVEC &sub_vec)
+ {
+ for (std::size_t i=0; i<pos.size(); ++i) {
+ this->add(pos[i], sub_vec[i]);
+ }
+ }
+
+ /// printout this equation for debugging
+ void write (const std::string &filename) const { std::ofstream os(filename); os << _vec; }
+
+ /// return a raw Lis vector object
+ RawVectorType& getRawVector() {return _vec; }
+
+ /// return a raw Lis vector object
+ const RawVectorType& getRawVector() const {return _vec; }
+
+ /// vector operation: set data
+ EigenVector& operator= (const EigenVector &src) { _vec = static_cast<const EigenVector&>(src)._vec; return *this; }
+
+ /// vector operation: add
+ void operator+= (const EigenVector& v) { _vec += static_cast<const EigenVector&>(v)._vec; }
+
+ /// vector operation: subtract
+ void operator-= (const EigenVector& v) { _vec -= static_cast<const EigenVector&>(v)._vec; }
+
+private:
+ RawVectorType _vec;
+};
+
+} // MathLib
+
+#endif //EIGENVECTOR_H_
+
diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
new file mode 100644
index 00000000000..b2254372238
--- /dev/null
+++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
@@ -0,0 +1,82 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EigenLisLinearSolver.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+#include <logog/include/logog.hpp>
+
+#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
+#include "MathLib/LinAlg/Lis/LisMatrix.h"
+#include "MathLib/LinAlg/Lis/LisVector.h"
+#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+EigenLisLinearSolver::EigenLisLinearSolver(EigenMatrix &A, ptree const*const option)
+: _A(A)
+{
+ if (option)
+ setOption(*option);
+}
+
+void EigenLisLinearSolver::setOption(const ptree &option)
+{
+ boost::optional<ptree> ptSolver = option.get_child("LinearSolver");
+ if (!ptSolver)
+ return;
+
+ boost::optional<std::string> solver_type = ptSolver->get_optional<std::string>("solver_type");
+ if (solver_type) {
+ _option.solver_type = _option.getSolverType(*solver_type);
+ }
+ boost::optional<std::string> precon_type = ptSolver->get_optional<std::string>("precon_type");
+ if (precon_type) {
+ _option.precon_type = _option.getPreconType(*precon_type);
+ }
+ boost::optional<double> error_tolerance = ptSolver->get_optional<double>("error_tolerance");
+ if (error_tolerance) {
+ _option.error_tolerance = *error_tolerance;
+ }
+ boost::optional<int> max_iteration_step = ptSolver->get_optional<int>("max_iteration_step");
+ if (max_iteration_step) {
+ _option.max_iterations = *max_iteration_step;
+ }
+}
+
+void EigenLisLinearSolver::solve(EigenVector &b_, EigenVector &x_)
+{
+ auto &A = _A.getRawMatrix();
+ auto &b = b_.getRawVector();
+ auto &x = x_.getRawVector();
+
+ A.makeCompressed();
+ int nnz = A.nonZeros();
+ int* ptr = A.outerIndexPtr();
+ int* col = A.innerIndexPtr();
+ double* data = A.valuePtr();
+ LisMatrix lisA(_A.getNRows(), nnz, ptr, col, data);
+ LisVector lisb(b.rows(), b.data());
+ LisVector lisx(x.rows(), x.data());
+
+ LisLinearSolver lissol(lisA);
+ lissol.setOption(_option);
+ lissol.solve(lisb, lisx);
+
+ for (std::size_t i=0; i<lisx.size(); i++)
+ x[i] = lisx[i];
+}
+
+} //MathLib
diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
new file mode 100644
index 00000000000..796dd8b44bc
--- /dev/null
+++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
@@ -0,0 +1,76 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EIGENLISLINEARSOLVER_H_
+#define EIGENLISLINEARSOLVER_H_
+
+#include <vector>
+
+#include <boost/property_tree/ptree.hpp>
+#include <lis.h>
+
+#include "MathLib/LinAlg/Lis/LisOption.h"
+
+namespace MathLib
+{
+class EigenVector;
+class EigenMatrix;
+
+/**
+ * Linear solver using Lis library with Eigen matrix and vector objects
+ */
+class EigenLisLinearSolver
+{
+public:
+ /**
+ * Constructor
+ * @param A Coefficient matrix object
+ * @param option A pointer to a linear solver option. In case you omit
+ * this argument, default settings follow those of
+ * LisOption struct.
+ */
+ EigenLisLinearSolver(EigenMatrix &A, boost::property_tree::ptree const*const option = nullptr);
+
+ virtual ~EigenLisLinearSolver() {}
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const boost::property_tree::ptree &option);
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const LisOption &option) { _option = option; }
+
+ /**
+ * get linear solver options
+ * @return
+ */
+ LisOption &getOption() { return _option; }
+
+ /**
+ * solve a given linear equations
+ *
+ * @param b RHS vector
+ * @param x Solution vector
+ */
+ void solve(EigenVector &b, EigenVector &x);
+
+private:
+ EigenMatrix& _A;
+ LisOption _option;
+};
+
+} // MathLib
+
+#endif //EIGENLISLINEARSOLVER_H_
+
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index d9017896b55..430704887e4 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -19,6 +19,10 @@
#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
+#ifdef OGS_USE_EIGEN
+#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#endif
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisMatrix.h"
#endif
@@ -169,6 +173,14 @@ TEST(Math, CheckInterface_GlobalDenseMatrix)
checkGlobalMatrixInterface(m);
}
+#ifdef OGS_USE_EIGEN
+TEST(Math, CheckInterface_EigenMatrix)
+{
+ MathLib::EigenMatrix m(10);
+ checkGlobalMatrixInterface(m);
+}
+#endif
+
#ifdef USE_LIS
TEST(Math, CheckInterface_LisMatrix)
{
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 3ad962468ff..caa06479069 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -19,6 +19,10 @@
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
+#ifdef OGS_USE_EIGEN
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
+#endif
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisVector.h"
#endif
@@ -195,6 +199,13 @@ TEST(Math, CheckInterface_DenseVector)
checkGlobalVectorInterface<MathLib::DenseVector<double> >();
}
+#ifdef OGS_USE_EIGEN
+TEST(Math, CheckInterface_EigenVector)
+{
+ checkGlobalVectorInterface<MathLib::EigenVector >();
+}
+#endif
+
#ifdef USE_LIS
TEST(Math, CheckInterface_LisVector)
{
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 7c5415b3df0..95d9569ec0a 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -23,6 +23,16 @@
#include "MathLib/LinAlg/ApplyKnownSolution.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+#ifdef OGS_USE_EIGEN
+#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
+#include "MathLib/LinAlg/Eigen/EigenLinearSolver.h"
+#endif
+
+#if defined(OGS_USE_EIGEN) && defined(USE_LIS)
+#include "MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h"
+#endif
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
#endif
@@ -201,6 +211,40 @@ TEST(MathLib, CheckInterface_GaussAlgorithm)
checkLinearSolverInterface<MathLib::GlobalDenseMatrix<double>, MathLib::DenseVector<double>, LinearSolverType>(A, t_root);
}
+#ifdef OGS_USE_EIGEN
+TEST(Math, CheckInterface_Eigen)
+{
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+ boost::property_tree::ptree t_solver;
+ t_solver.put("solver_type", "CG");
+ t_solver.put("precon_type", "NONE");
+ t_solver.put("error_tolerance", 1e-15);
+ t_solver.put("max_iteration_step", 1000);
+ t_root.put_child("LinearSolver", t_solver);
+
+ MathLib::EigenMatrix A(Example1::dim_eqs);
+ checkLinearSolverInterface<MathLib::EigenMatrix, MathLib::EigenVector, MathLib::EigenLinearSolver>(A, t_root);
+}
+#endif
+
+#if defined(OGS_USE_EIGEN) && defined(USE_LIS)
+TEST(Math, CheckInterface_EigenLis)
+{
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+ boost::property_tree::ptree t_solver;
+ t_solver.put("solver_type", "CG");
+ t_solver.put("precon_type", "NONE");
+ t_solver.put("error_tolerance", 1e-15);
+ t_solver.put("max_iteration_step", 1000);
+ t_root.put_child("LinearSolver", t_solver);
+
+ MathLib::EigenMatrix A(Example1::dim_eqs);
+ checkLinearSolverInterface<MathLib::EigenMatrix, MathLib::EigenVector, MathLib::EigenLisLinearSolver>(A, t_root);
+}
+#endif
+
#ifdef USE_LIS
TEST(Math, CheckInterface_Lis)
{
--
GitLab