Merge pull request #64 from norihiro-w/merge-LinearSystem2

System of linear equations inteface and Lis solver interface.

BaseLib/TemplateLogogFormatterSuppressedGCC.h

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file TemplateLogogFormatterSuppressedGCC.h
+ *
+ * Created on 2012-12-06 by Norihiro Watanabe
+ */
+
+#ifndef TEMPLATELOGOGFORMATTERSUPPRESSEDGCC_H_
+#define TEMPLATELOGOGFORMATTERSUPPRESSEDGCC_H_
+
+// ** INCLUDES **
+#include "logog/include/logog.hpp"
+
+namespace BaseLib {
+
+/**
+ * \brief TemplateLogogFormatterSuppressedGCC strips topics given as a template
+ * parameter from logog::FormatterGCC.
+ * See http://johnwbyrd.github.com/logog/customformatting.html for details.
+ **/
+template <int T_SUPPPRESS_TOPIC_FLAG>
+class TemplateLogogFormatterSuppressedGCC : public logog::FormatterGCC
+{
+
+	virtual TOPIC_FLAGS GetTopicFlags( const logog::Topic &topic )
+	{
+	return ( logog::Formatter::GetTopicFlags( topic ) &
+		~( T_SUPPPRESS_TOPIC_FLAG ));
+	}
+
+};
+
+#endif // TEMPLATELOGOGFORMATTERSUPPRESSEDGCC_H_
+
+} // namespace BaseLib

CMakeLists.txt

@@ -55,6 +55,9 @@ OPTION(OGS_BUILD_UTILS "Should the utilities programms be built?" OFF)
 
 OPTION(OGS_NO_EXTERNAL_LIBS "Builds OGS without any external dependencies." OFF)
 
+# Linear solvers
+OPTION(OGS_USE_LIS "Use Lis" OFF)
+
 # Logging
 OPTION(OGS_DISABLE_LOGGING "Disables all logog messages." OFF)
 
@@ -98,6 +101,10 @@ INCLUDE_DIRECTORIES( ${CMAKE_SOURCE_DIR}/ThirdParty )
 INCLUDE_DIRECTORIES( ${CMAKE_SOURCE_DIR}/ThirdParty/quickcheck )
 INCLUDE_DIRECTORIES( ${CMAKE_BINARY_DIR}/ThirdParty/zlib )
 
+IF(OGS_USE_LIS)
+    ADD_DEFINITIONS(-DUSE_LIS)
+ENDIF()
+
 ADD_SUBDIRECTORY( BaseLib )
 # TODO This is a hack but we have to make sure that Boost is built first
 ADD_DEPENDENCIES(BaseLib Boost)

MathLib/CMakeLists.txt

@@ -20,6 +20,16 @@ SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
 GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
 SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
 
+GET_SOURCE_FILES(SOURCES_LINALG_LEQS LinAlg/SystemOfLinearEquations)
+IF (NOT OGS_USE_LIS)
+    LIST(REMOVE_ITEM SOURCES_LINALG_LEQS 
+        LinAlg/SystemOfLinearEquations/LisOption.h
+        LinAlg/SystemOfLinearEquations/LisOption.cpp
+        LinAlg/SystemOfLinearEquations/LisLinearSystem.h
+        LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp
+    )
+ENDIF ()
+SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_LEQS})
 
 IF (METIS_FOUND)
 	GET_SOURCE_FILES(SOURCES_LINALG_SPARSE_NESTEDDISSECTION LinAlg/Sparse/NestedDissectionPermutation)
@@ -38,8 +48,17 @@ IF(METIS_FOUND)
 	)
 ENDIF()
 
+IF (LIS_FOUND)
+    INCLUDE_DIRECTORIES(${LIS_INCLUDE_DIR})
+ENDIF()
+
+
 # Create the library
 ADD_LIBRARY( MathLib STATIC ${SOURCES} )
 
 SET_TARGET_PROPERTIES(MathLib PROPERTIES LINKER_LANGUAGE CXX)
 
+IF (LIS_FOUND)
+    TARGET_LINK_LIBRARIES( MathLib ${LIS_LIBRARIES} )
+ENDIF()
+

MathLib/LinAlg/Sparse/Sparsity.h

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file Sparsity.h
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#ifndef SPARSITY_H_
+#define SPARSITY_H_
+
+#include <vector>
+#include <set>
+
+namespace MathLib 
+{
+
+/**
+ * \brief Row-major sparse pattern
+ */
+typedef std::vector<std::set<std::size_t> > RowMajorSparsity;
+
+} //MathLib
+
+#endif // SPARSITY_H_

MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.h

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file ISystemOfLinearEquations.h
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#ifndef ISYSTEMOFLINEAREQUATIONS_H_
+#define ISYSTEMOFLINEAREQUATIONS_H_
+
+#include <iostream>
+#include <vector>
+#include <boost/property_tree/ptree.hpp>
+#include "MathLib/LinAlg/Sparse/Sparsity.h"
+
+namespace MathLib
+{
+
+/**
+ * \brief Interface to any system of linear equations Ax=b
+ *
+ * This class defines unified interface to a system of linear equations
+ * including linear solvers.
+ *
+ */
+class ISystemOfLinearEquations
+{
+public:
+    /// Index representing invalid
+    const static std::size_t index_npos = -1;
+
+    /// 
+    virtual ~ISystemOfLinearEquations()
+    {
+    }
+
+    /**
+     * set properties
+     *
+     * @param option
+     */
+    virtual void setOption(const boost::property_tree::ptree &option) = 0;
+
+    /**
+     * initialize this system with zero
+     *
+     * This function initialize the system by setting zero to all entries in
+     * a matrix, a RHS vector and solution vector. Dimension of the system is
+     * not changed.
+     */
+    virtual void setZero() = 0;
+
+    /// return dimension of this equation
+    virtual std::size_t getDimension() const = 0;
+
+    /**
+     * set an entry in a matrix
+     *
+     * @param rowId
+     * @param colId
+     * @param v
+     */
+    virtual void setMatEntry(std::size_t rowId, std::size_t colId, double v) = 0;
+
+    /**
+     * add a value into a matrix
+     *
+     * @param rowId
+     * @param colId
+     * @param v
+     */
+    virtual void addMatEntry(std::size_t rowId, std::size_t colId, double v) = 0;
+
+    /**
+     * add a sub-matrix into a matrix at given row and column positions
+     *
+     * @param vec_row_pos
+     *  A vector of global row index. The number of entries of vec_row_pos have
+     *  to be identical to the number of rows of the local matrix.
+     * @param vec_col_pos
+     *  A vector of global column index. The number of vec_col_pos have to be
+     *  identical to the number of columns of the local matrix.
+     * @param sub_matrix    A sub-matrix
+     * @param fac           scaling parameter
+     */
+    template <class T_DENSE_MATRIX>
+    void addSubMat(const std::vector<std::size_t> &vec_row_pos,
+            const std::vector<std::size_t> &vec_col_pos,
+            const T_DENSE_MATRIX &sub_matrix, double fac = 1.0);
+
+    /**
+     * add a sub-matrix into a matrix
+     *
+     * @param vec_row_col_pos
+     *  A vector of global row and column index. The number of entries of
+     *  vec_row_col_pos have to be identical to the number of rows and columns
+     *  of the local matrix.
+     * @param sub_matrix        A sub-matrix
+     * @param fac               scaling parameter
+     */
+    template <class T_DENSE_MATRIX>
+    void addSubMat(const std::vector<std::size_t> &vec_row_col_pos,
+            const T_DENSE_MATRIX &sub_matrix, double fac = 1.0);
+
+    /**
+     * get RHS entry
+     *
+     * @param rowId
+     * @return
+     */
+    virtual double getRHSVec(std::size_t rowId) const = 0;
+
+    /**
+     * set RHS entry
+     *
+     * @param rowId
+     * @param v
+     */
+    virtual void setRHSVec(std::size_t rowId, double v) = 0;
+
+    /**
+     * add RHS entry
+     *
+     * @param rowId
+     * @param v
+     */
+    virtual void addRHSVec(std::size_t rowId, double v) = 0;
+
+    /**
+     * add a sub vector to RHS
+     *
+     * @param vec_row_pos
+     *  A vector of global row index. The number of entries of vec_row_pos
+     *  have to be identical to the number of rows of the sub vector.
+     * @param sub_vector
+     *  Pointer to a sub-vector. Its length must be the same as the length of vec_row_pos.
+     * @param fac
+     *  Scaling parameter. Default value is 1.
+     */
+    inline virtual void addSubRHS(const std::vector<std::size_t> &vec_row_pos,
+            const double* sub_vector, double fac = 1.0);
+
+    /**
+     * add a sub vector to RHS
+     *
+     * @param vec_row_pos
+     *  A vector of global row index. The number of entries of vec_row_pos
+     *  have to be identical to the number of rows of the sub vector.
+     * @param sub_vector    A sub-vector
+     * @param fac           Scaling parameter
+     */
+    template <class T_DENSE_VECTOR>
+    void addSubRHS(const std::vector<std::size_t> &vec_row_pos,
+            const T_DENSE_VECTOR &sub_vector, double fac = 1.0);
+
+    /**
+     * get an entry in a solution vector
+     * @param rowId
+     * @return
+     */
+    virtual double getSolVec(std::size_t rowId) = 0;
+
+    /**
+     * set an entry in a solution vector
+     * @param rowId
+     * @param v
+     */
+    virtual void setSolVec(std::size_t rowId, double v) = 0;
+
+    /**
+     * set prescribed value
+     * @param row_id
+     * @param x
+     */
+    virtual void setKnownSolution(std::size_t row_id, double x) = 0;
+
+    /**
+     * set prescribed value
+     * @param vec_id    A vector of global row index
+     * @param vec_x     A vector of prescribed values
+     */
+    virtual void setKnownSolution(const std::vector<std::size_t> &vec_id,
+            const std::vector<double> &vec_x) = 0;
+
+    /// solve the linear system
+    virtual void solve() = 0;
+
+    /// print out the equation for debugging
+    virtual void printout(std::ostream &os = std::cout) const = 0;
+
+};
+
+}
+
+#include "ISystemOfLinearEquations.tpp"
+
+#endif //ISYSTEMOFLINEAREQUATIONS_H_

MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.tpp

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file ISystemOfLinearEquations.tpp
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#ifndef ISYSTEMOFLINEAREQUATIONS_TPP_
+#define ISYSTEMOFLINEAREQUATIONS_TPP_
+
+namespace MathLib
+{
+
+template <class T_DENSE_MATRIX>
+void ISystemOfLinearEquations::addSubMat(const std::vector<std::size_t> &vec_row_pos, const std::vector<std::size_t> &vec_col_pos, const T_DENSE_MATRIX &sub_matrix, double fkt)
+{
+    const std::size_t n_rows = vec_row_pos.size();
+    const std::size_t n_cols = vec_col_pos.size();
+    for (std::size_t i=0; i<n_rows; i++) {
+        const std::size_t rowId = vec_row_pos[i];
+        if (rowId==index_npos) continue;
+        for (std::size_t j=0; j<n_cols; j++) {
+            const std::size_t colId = vec_col_pos[j];
+            if (colId==index_npos) continue;
+            addMatEntry(rowId, colId, fkt*sub_matrix(i,j));
+        }
+    }
+}
+
+template <class T_DENSE_MATRIX>
+void ISystemOfLinearEquations::addSubMat(const std::vector<std::size_t> &vec_pos, const T_DENSE_MATRIX &sub_matrix, double fkt)
+{
+    addSubMat(vec_pos, vec_pos, sub_matrix, fkt);
+}
+
+inline void ISystemOfLinearEquations::addSubRHS(const std::vector<std::size_t> &vec_row_pos, const double *sub_vector, double fkt)
+{
+    for (std::size_t i=0; i<vec_row_pos.size(); i++) {
+        const std::size_t rowId = vec_row_pos[i];
+        if (rowId==index_npos) continue;
+        addRHSVec(rowId, sub_vector[i]*fkt);
+    }
+}
+
+template <class T_DENSE_VECTOR>
+void ISystemOfLinearEquations::addSubRHS(const std::vector<std::size_t> &vec_row_pos, const T_DENSE_VECTOR &sub_vector, double fkt)
+{
+    for (std::size_t i=0; i<vec_row_pos.size(); i++) {
+        const std::size_t rowId = vec_row_pos[i];
+        if (rowId==index_npos) continue;
+        addRHSVec(rowId, sub_vector(i)*fkt);
+    }
+}
+
+}
+
+#endif //ISYSTEMOFLINEAREQUATIONS_TPP_
+

MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file LisLinearSystem.cpp
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#include "LisLinearSystem.h"
+
+#include <string>
+#include <sstream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+#include "logog/include/logog.hpp"
+
+namespace MathLib
+{
+
+bool LisLinearSystem::checkError(int err)
+{
+    bool ok = (err == LIS_SUCCESS);
+    if (!ok) {
+        ERR("***ERROR: Lis error code = %d", err);
+    }
+    return ok;
+}
+
+LisLinearSystem::LisLinearSystem(std::size_t dimension)
+: _dim(dimension), _max_diag_coeff(.0)
+{
+    int ierr = 0;
+    ierr = lis_matrix_create(0, &_AA); checkError(ierr);
+    ierr = lis_matrix_set_size(_AA, 0, dimension); checkError(ierr);
+    ierr = lis_vector_duplicate(_AA, &_bb); checkError(ierr);
+    ierr = lis_vector_duplicate(_AA, &_xx); checkError(ierr);
+}
+
+LisLinearSystem::~LisLinearSystem()
+{
+    int ierr = 0;
+    ierr = lis_matrix_destroy(_AA); checkError(ierr);
+    ierr = lis_vector_destroy(_bb); checkError(ierr);
+    ierr = lis_vector_destroy(_xx); checkError(ierr);
+}
+
+void LisLinearSystem::setOption(const boost::property_tree::ptree &option)
+{
+    using boost::property_tree::ptree;
+
+    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<std::string> matrix_type = ptSolver->get_optional<std::string>("matrix_type");
+    if (matrix_type) {
+        _option.matrix_type = _option.getMatrixType(*matrix_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 LisLinearSystem::setZero()
+{
+    int ierr = 0;
+    // A matrix has to be removed and created every time because Lis doesn't provide a
+    // function to set matrix entries to zero
+    ierr = lis_matrix_destroy(_AA); checkError(ierr);
+    ierr = lis_matrix_create(0, &_AA); checkError(ierr);
+    ierr = lis_matrix_set_size(_AA, 0, getDimension()); checkError(ierr);
+    // set zero RHS, x
+    ierr = lis_vector_set_all(0.0, _bb); checkError(ierr);
+    ierr = lis_vector_set_all(0.0, _xx); checkError(ierr);
+
+    _max_diag_coeff = .0;
+}
+
+void LisLinearSystem::applyKnownSolution()
+{
+    //Use penalty parameter
+    const double penalty_scaling = 1e+10;
+    const double penalty = _max_diag_coeff * penalty_scaling;
+    INFO("-> max. absolute value of diagonal entries = %e", _max_diag_coeff);
+    INFO("-> penalty scaling = %e", penalty_scaling);
+    const std::size_t n_bc = _vec_knownX_id.size();
+    for (std::size_t i_bc=0; i_bc<n_bc; i_bc++) {
+        const std::size_t rowId = _vec_knownX_id[i_bc];
+        const double x = _vec_knownX_x[i_bc];
+        //A(k, k) = penalty
+        setMatEntry(rowId, rowId, penalty);
+        //b(k) = x*penalty
+        setRHSVec(rowId, x*penalty);
+    }
+}
+
+void LisLinearSystem::solve()
+{
+    INFO("------------------------------------------------------------------");
+    INFO("*** LIS solver computation");
+#ifdef _OPENMP
+    INFO("-> max number of threads = %d", omp_get_num_procs());
+    INFO("-> number of threads = %d", omp_get_max_threads());
+#endif
+
+    applyKnownSolution();
+    // assemble a matrix
+    int ierr = 0;
+    ierr = lis_matrix_set_type(_AA, static_cast<int>(_option.matrix_type)); checkError(ierr);
+    ierr = lis_matrix_assemble(_AA); checkError(ierr);
+
+    // configure option
+    std::string solver_options;
+    if (_option.solver_precon_arg.empty()) {
+        std::stringstream ss;
+        ss << "-i " << static_cast<int>(_option.solver_type);
+        ss << " -p " << static_cast<int>(_option.precon_type);
+        if (!_option.extra_arg.empty())
+            ss << " " << _option.extra_arg;
+        solver_options = ss.str();
+    } else {
+        solver_options = _option.solver_precon_arg;
+    }
+    std::string tol_option;
+    {
+        std::stringstream ss;
+        ss << "-tol " << _option.error_tolerance;
+        ss << " -maxiter " << _option.max_iterations;
+        ss << " -initx_zeros 0"; //0: use given x as initial guess, 1: x0=0
+#ifdef _OPENMP
+        const int nthreads = omp_get_max_threads();
+        ss << " -omp_num_threads " << nthreads;
+#endif
+        tol_option = ss.str();
+    }
+
+    // Create solver
+    LIS_SOLVER solver;
+    ierr = lis_solver_create(&solver); checkError(ierr);
+    ierr = lis_solver_set_option(const_cast<char*>(solver_options.c_str()), solver); checkError(ierr);
+    ierr = lis_solver_set_option(const_cast<char*>(tol_option.c_str()), solver); checkError(ierr);
+    ierr = lis_solver_set_option(const_cast<char*>("-print mem"), solver); checkError(ierr);
+    
+    // solve
+    INFO("-> solve");
+    ierr = lis_solve(_AA, _bb, _xx, solver); checkError(ierr);
+    checkError(ierr);
+    //lis_output(_AA, _bb, _xx, LIS_FMT_MM, "/home/localadmin/tasks/20120814_ogs6test/matrix1.txt");
+
+    int iter = 0;
+    double resid = 0.0;
+    ierr = lis_solver_get_iters(solver, &iter); checkError(ierr);
+    ierr = lis_solver_get_residualnorm(solver, &resid); checkError(ierr);
+    INFO("\t iteration: %d/%ld\n", iter, _option.max_iterations);
+    INFO("\t residuals: %e\n", resid);
+
+    // Clear solver
+    ierr = lis_solver_destroy(solver); checkError(ierr);
+    INFO("------------------------------------------------------------------");
+}
+
+void LisLinearSystem::printout(std::ostream &os) const
+{
+    os << "#A=" << std::endl;
+    os << "#x=" << std::endl;
+    lis_vector_print(_xx);
+    os << "#b=" << std::endl;
+    lis_vector_print(_bb);
+}
+
+} //MathLib

MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file LisLinearSystem.h
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#ifndef LISLINEARSYSTEM_H_
+#define LISLINEARSYSTEM_H_
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include <boost/property_tree/ptree.hpp>
+
+#include "lis.h"
+
+#include "LisOption.h"
+
+namespace MathLib
+{
+
+/**
+ * \brief Linear system using Lis solver (http://www.ssisc.org/lis/)
+ *
+ * This class utilizes Lis solver.
+ */
+class LisLinearSystem
+{
+public:
+    //---------------------------------------------------------------
+    // realization of ISystemOfLinearEquations
+    //---------------------------------------------------------------
+
+    /**
+     * Create a linear system using Lis
+     *
+     * @param length    System dimension
+     */
+    LisLinearSystem(std::size_t length);
+
+    /**
+     *
+     */
+    virtual ~LisLinearSystem();
+
+    /**
+     * configure linear solvers
+     * @param option
+     */
+    virtual void setOption(const boost::property_tree::ptree &option);
+
+    /// return the system dimension
+    virtual std::size_t getDimension() const { return _dim; };
+
+    /// reset this equation
+    virtual void setZero();
+
+    /// set entry in A
+    virtual void setMatEntry(std::size_t rowId, std::size_t colId, double v)
+    {
+        if (rowId==colId)
+            _max_diag_coeff = std::max(_max_diag_coeff, std::abs(v));
+        lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
+    }
+
+    /// add value into A
+    virtual void addMatEntry(std::size_t rowId, std::size_t colId, double v)
+    {
+        if (rowId==colId)
+            _max_diag_coeff = std::max(_max_diag_coeff, std::abs(v));
+        lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
+    }
+
+    /// get RHS entry
+    virtual double getRHSVec(std::size_t rowId) const
+    {
+        double v;
+        lis_vector_get_value(_bb, rowId, &v);
+        return v;
+    }
+
+    /// set RHS entry
+    virtual void setRHSVec(std::size_t rowId, double v)
+    {
+        lis_vector_set_value(LIS_INS_VALUE, rowId, v, _bb);
+    }
+
+    /// add RHS entry
+    virtual void addRHSVec(std::size_t rowId, double v)
+    {
+        lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _bb);
+    }
+
+    /// get an entry in a solution vector
+    virtual double getSolVec(std::size_t rowId)
+    {
+        double v;
+        lis_vector_get_value(_xx, rowId, &v);
+        return v;
+    }
+
+    /// set a solution vector
+    virtual void setSolVec(std::size_t rowId, double v)
+    {
+        lis_vector_set_value(LIS_INS_VALUE, rowId, v, _xx);
+    }
+
+    /// set prescribed value
+    virtual void setKnownSolution(std::size_t rowId, double x)
+    {
+        _vec_knownX_id.push_back(rowId);
+        _vec_knownX_x.push_back(x);
+    }
+
+    /// set prescribed values
+    virtual void setKnownSolution(const std::vector<std::size_t> &vec_id, const std::vector<double> &vec_x)
+    {
+        _vec_knownX_id.insert(_vec_knownX_id.end(), vec_id.begin(), vec_id.end());
+        _vec_knownX_x.insert(_vec_knownX_x.end(), vec_x.begin(), vec_x.end());
+    }
+
+    /// solve this equation
+    virtual void solve();
+
+    /// printout this equation for debugging
+    virtual void printout(std::ostream &os=std::cout) const;
+
+    //---------------------------------------------------------------
+    // specific to this class
+    //---------------------------------------------------------------
+    /**
+     * configure linear solvers
+     * @param option
+     */
+    void setOption(const LisOption &option)
+    {
+        _option = option;
+    }
+
+    /**
+     * get linear solver options
+     * @return
+     */
+    LisOption &getOption()
+    {
+        return _option;
+    }
+
+    /**
+     * get a solution vector
+     * @param x
+     */
+    void getSolVec(double* x)
+    {
+        lis_vector_gather(_xx, x);
+    }
+
+private:
+    bool checkError(int err);
+    void applyKnownSolution();
+
+private:
+    const std::size_t _dim;
+    double _max_diag_coeff;
+    LisOption _option;
+    LIS_MATRIX _AA;
+    LIS_VECTOR _bb;
+    LIS_VECTOR _xx;
+    std::vector<std::size_t> _vec_knownX_id;
+    std::vector<double> _vec_knownX_x;
+};
+
+
+} // MathLib
+
+#endif //LISLINEARSYSTEM_H_
+

MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file LisOption.cpp
+ *
+ * Created on 2012-12-05 by Norihiro Watanabe
+ */
+
+#include "LisOption.h"
+
+namespace MathLib
+{
+
+LisOption::LisOption()
+{
+    solver_type = SolverType::CG;
+    precon_type = PreconType::NONE;
+    matrix_type = MatrixType::CRS;
+    max_iterations = 500;
+    error_tolerance = 1.e-10;
+}
+
+LisOption::SolverType LisOption::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, BiCG);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, CGS);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCGSTAB);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCGSTABl);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, GPBiCG);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, TFQMR);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, Orthomin);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, GMRES);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, Jacobi);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, GaussSeidel);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, SOR);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCGSafe);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, CR);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCR);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, CRS);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCRSTAB);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, GPBiCR);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCRSafe);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, FGMRESm);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, IDRs);
+    RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, MINRES);
+
+    return SolverType::INVALID;
+#undef RETURN_SOLVER_ENUM_IF_SAME_STRING
+}
+
+LisOption::PreconType LisOption::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, JACOBI);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, ILU);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, SSOR);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, Hybrid);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, IplusS);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, SAINV);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, SAAMG);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, CroutILU);
+    RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, ILUT);
+
+    return PreconType::NONE;
+#undef RETURN_PRECOM_ENUM_IF_SAME_STRING
+}
+
+LisOption::MatrixType LisOption::getMatrixType(const std::string &matrix_name)
+{
+#define RETURN_MATRIX_ENUM_IF_SAME_STRING(str, TypeName) \
+    if (#TypeName==str) return MatrixType::TypeName;
+
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, CRS);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, CCS);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, MSR);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, DIA);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, ELL);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, JDS);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, BSR);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, BSC);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, VBR);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, COO);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, DNS);
+
+    return MatrixType::CRS;
+#undef RETURN_MATRIX_ENUM_IF_SAME_STRING
+}
+
+} //MathLib

MathLib/LinAlg/SystemOfLinearEquations/LisOption.h

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file LisOption.h
+ *
+ * Created on 2012-06-25 by Norihiro Watanabe
+ */
+
+#ifndef LIS_OPTION_H_
+#define LIS_OPTION_H_
+
+#include <string>
+
+namespace MathLib
+{
+
+/**
+ * \brief Option for Lis solver
+ */
+struct LisOption
+{
+    /// Solver type
+    enum class SolverType : int
+    {
+        INVALID = 0,
+        CG = 1,
+        BiCG = 2,
+        CGS = 3,
+        BiCGSTAB = 4,
+        BiCGSTABl = 5,
+        GPBiCG = 6,
+        TFQMR = 7,
+        Orthomin = 8,
+        GMRES = 9,
+        Jacobi = 10,
+        GaussSeidel = 11,
+        SOR = 12,
+        BiCGSafe = 13,
+        CR = 14,
+        BiCR = 15,
+        CRS = 16,
+        BiCRSTAB = 17,
+        GPBiCR = 18,
+        BiCRSafe = 19,
+        FGMRESm = 20,
+        IDRs = 21,
+        MINRES = 22
+    };
+
+    /// Preconditioner type
+    enum class PreconType : int
+    {
+        NONE = 0,
+        JACOBI = 1,
+        ILU = 2,
+        SSOR = 3,
+        Hybrid = 4,
+        IplusS = 5,
+        SAINV = 6,
+        SAAMG = 7,
+        CroutILU = 8,
+        ILUT = 9
+    };
+
+    /// Matrix type
+    enum class MatrixType : int
+    {
+        CRS = 1,
+        CCS = 2,
+        MSR = 3,
+        DIA = 4,
+        ELL = 5,
+        JDS = 6,
+        BSR = 7,
+        BSC = 8,
+        VBR = 9,
+        COO = 10,
+        DNS = 11
+    };
+
+    /// Linear solver type
+    SolverType solver_type;
+    /// Preconditioner type
+    PreconType precon_type;
+    /// Matrix type
+    MatrixType matrix_type;
+    /// Maximum iteration count
+    long max_iterations;
+    /// Error tolerance
+    double error_tolerance;
+    /// Extra option
+    std::string extra_arg;
+    /// Arguments for solver and preconditioner. This variable is always preferred
+    /// to other variables.
+    std::string solver_precon_arg;
+
+    /**
+     * Constructor
+     *
+     * Default options are CG, no preconditioner, iteration count 500 and
+     * tolerance 1e-10. Default matrix storage type is CRS.
+     */
+    LisOption();
+
+    /// Destructor
+    ~LisOption() {};
+
+    /**
+     * 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);
+
+    /**
+     * return a matrix type
+     * @param matrix_name
+     * @return
+     *      If there is no matrix type matched with the given name, CRS
+     *      is returned.
+     */
+    static MatrixType getMatrixType(const std::string &matrix_name);
+};
+
+}
+#endif //LIS_OPTION_H_

Tests/CMakeLists.txt

@@ -7,13 +7,19 @@ ENDIF()
 APPEND_SOURCE_FILES(TEST_SOURCES)
 APPEND_SOURCE_FILES(TEST_SOURCES BaseLib)
 APPEND_SOURCE_FILES(TEST_SOURCES GeoLib)
+APPEND_SOURCE_FILES(TEST_SOURCES MathLib)
 
 INCLUDE_DIRECTORIES(
 	${CMAKE_SOURCE_DIR}/Tests/gtest
 	${CMAKE_SOURCE_DIR}/BaseLib
 	${CMAKE_SOURCE_DIR}/GeoLib
+	${CMAKE_SOURCE_DIR}/MathLib
 )
 
+IF (LIS_FOUND)
+    INCLUDE_DIRECTORIES(${LIS_INCLUDE_DIR})
+ENDIF()
+
 ADD_EXECUTABLE (testrunner testrunner.cpp ${TEST_SOURCES})
 SET_TARGET_PROPERTIES(testrunner PROPERTIES FOLDER Testing)
 
@@ -21,6 +27,8 @@ TARGET_LINK_LIBRARIES(testrunner
 	GTest
 	BaseLib
 	GeoLib
+	MathLib
+	logog
 	${Boost_LIBRARIES}
 	${CMAKE_THREAD_LIBS_INIT}
 )

Tests/MathLib/TestLinearSystem.cpp

+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.com)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.com/LICENSE.txt
+ *
+ *
+ * \file TestLinearSystem.cpp
+ *
+ * Created on 2012-08-03 by Norihiro Watanabe
+ */
+
+#include <gtest/gtest.h>
+#include <boost/property_tree/ptree.hpp>
+
+#include "MathLib/LinAlg/Sparse/Sparsity.h"
+#ifdef USE_LIS
+#include "MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h"
+#endif
+
+namespace
+{
+
+inline void ASSERT_DOUBLE_ARRAY_EQ(const double* Expected, const double* Actual, size_t N, double epsilon=1.0e-8) {
+    for (size_t i=0; i<N; i++) \
+        ASSERT_NEAR(Expected[i], Actual[i], epsilon);
+}
+
+struct Example1
+{
+    std::vector<double> mat;
+    std::vector<size_t> list_dirichlet_bc_id;
+    std::vector<double> list_dirichlet_bc_value;
+    static const size_t dim_eqs = 9;
+    std::vector<double> exH;
+    MathLib::RowMajorSparsity sparse;
+
+    Example1()
+    {
+        double d_mat[] = {
+            6.66667e-012, -1.66667e-012, 0, -1.66667e-012, -3.33333e-012, 0, 0, 0, 0, 
+            -1.66667e-012, 1.33333e-011, -1.66667e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 0, 0, 0, 
+            0, -1.66667e-012, 6.66667e-012, 0, -3.33333e-012, -1.66667e-012, 0, 0, 0, 
+            -1.66667e-012, -3.33333e-012, 0, 1.33333e-011, -3.33333e-012, 0, -1.66667e-012, -3.33333e-012, 0, 
+            -3.33333e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 2.66667e-011, -3.33333e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 
+            0, -3.33333e-012, -1.66667e-012, 0, -3.33333e-012, 1.33333e-011, 0, -3.33333e-012, -1.66667e-012, 
+            0, 0, 0, -1.66667e-012, -3.33333e-012, 0, 6.66667e-012, -1.66667e-012, 0, 
+            0, 0, 0, -3.33333e-012, -3.33333e-012, -3.33333e-012, -1.66667e-012, 1.33333e-011, -1.66667e-012, 
+            0, 0, 0, 0, -3.33333e-012, -1.66667e-012, 0, -1.66667e-012, 6.66667e-012
+        };
+        mat.assign(d_mat, d_mat+dim_eqs*dim_eqs);
+        size_t int_dirichlet_bc_id[] = {2,5,8,0,3,6};
+        list_dirichlet_bc_id.assign(int_dirichlet_bc_id, int_dirichlet_bc_id+6);
+        list_dirichlet_bc_value.resize(6);
+        fill(list_dirichlet_bc_value.begin(), list_dirichlet_bc_value.begin()+3, .0);
+        fill(list_dirichlet_bc_value.begin()+3, list_dirichlet_bc_value.end(), 1.0);
+        exH.resize(9);
+        for (size_t i=0; i<9; i++) {
+            if (i%3==0) exH[i] = 1.0;
+            if (i%3==1) exH[i] = 0.5;
+            if (i%3==2) exH[i] = 0.;
+        }
+        sparse.resize(dim_eqs);
+        for (size_t i=0; i<dim_eqs; i++) {
+            for (size_t j=0; j<dim_eqs; j++) {
+                if (mat[i*dim_eqs+j]!=.0)
+                    sparse[i].insert(j);
+            }
+        }
+
+    }
+};
+}
+
+#ifdef USE_LIS
+TEST(Math, LinearSystemLis)
+{
+    // set a problem
+    Example1 ex1;
+
+    // create a linear system
+    MathLib::LisLinearSystem eqs(ex1.dim_eqs);
+
+    // construct
+    for (size_t i=0; i<ex1.dim_eqs; i++) {
+        for (size_t j=0; j<ex1.dim_eqs; j++) {
+            double v = ex1.mat[i*ex1.dim_eqs+j];
+            if (v!=.0)
+                eqs.addMatEntry(i, j, v);
+        }
+    }
+
+    // apply BC
+    eqs.setKnownSolution(ex1.list_dirichlet_bc_id, ex1.list_dirichlet_bc_value);
+
+    // 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("matrix_type", "CCS");
+    t_solver.put("error_tolerance", 1e-15);
+    t_solver.put("max_iteration_step", 1000);
+    t_root.put_child("LinearSolver", t_solver);
+    eqs.setOption(t_root);
+
+    // check if the option was correctly parsed
+    MathLib::LisOption &lisOption = eqs.getOption();
+    ASSERT_EQ(MathLib::LisOption::SolverType::CG, lisOption.solver_type);
+    ASSERT_EQ(MathLib::LisOption::PreconType::NONE, lisOption.precon_type);
+    ASSERT_EQ(MathLib::LisOption::MatrixType::CCS, lisOption.matrix_type);
+    ASSERT_EQ(1e-15, lisOption.error_tolerance);
+    ASSERT_EQ(1000, lisOption.max_iterations);
+
+    // solve
+    eqs.solve();
+
+//    eqs.printout();
+
+    // check solution
+    std::vector<double> vec_x(eqs.getDimension());
+    eqs.getSolVec(&vec_x[0]);
+    ASSERT_DOUBLE_ARRAY_EQ(&ex1.exH[0], &vec_x[0], ex1.dim_eqs, 1.e-5);
+}
+
+#endif
+

Tests/testrunner.cpp

@@ -12,10 +12,43 @@
 
 // ** INCLUDES **
 #include "gtest/gtest.h"
+#include "logog/include/logog.hpp"
+#ifdef USE_LIS
+#include "lis.h"
+#endif
+#include "BaseLib/TemplateLogogFormatterSuppressedGCC.h"
 
 /// Implementation of the googletest testrunner
 int main(int argc, char* argv[])
 {
-	testing::InitGoogleTest ( &argc, argv );
-	return RUN_ALL_TESTS();
+    int ret = 0;
+    LOGOG_INITIALIZE();
+    {
+        logog::Cout out;
+        BaseLib::TemplateLogogFormatterSuppressedGCC<TOPIC_LEVEL_FLAG | TOPIC_FILE_NAME_FLAG | TOPIC_LINE_NUMBER_FLAG> custom_format;
+        out.SetFormatter(custom_format);
+
+        try {
+            // initialize libraries which will be used while testing
+#ifdef USE_LIS
+            lis_initialize(&argc, &argv);
+#endif
+            // start google test
+            testing::InitGoogleTest ( &argc, argv );
+            ret = RUN_ALL_TESTS();
+        } catch (char* e) {
+            ERR(e);
+        } catch (std::exception& e) {
+            ERR(e.what());
+        } catch (...) {
+            ERR("Unknown exception occurred!");
+        }
+        // finalize libraries
+#ifdef USE_LIS
+        lis_finalize();
+#endif
+    } // make sure no logog objects exist when LOGOG_SHUTDOWN() is called.
+    LOGOG_SHUTDOWN();
+
+    return ret;
 }

scripts/cmake/Find.cmake

@@ -138,6 +138,11 @@ IF(Shapelib_FOUND)
 	ADD_DEFINITIONS(-DShapelib_FOUND)
 ENDIF() # Shapelib_FOUND
 
+## lis ##
+IF(OGS_USE_LIS)
+    FIND_PACKAGE( LIS REQUIRED )
+ENDIF()
+
 ########################
 ### Find other stuff ###
 ########################

scripts/cmake/Functions.cmake

@@ -15,8 +15,9 @@ MACRO(GET_SOURCE_FILES SOURCE_FILES)
 	ENDIF()
 
 	# Get all files in the directory
-	FILE(GLOB GET_SOURCE_FILES_HEADERS ${DIR}/*.h)
-	FILE(GLOB GET_SOURCE_FILES_SOURCES ${DIR}/*.cpp)
+	FILE(GLOB GET_SOURCE_FILES_HEADERS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} ${DIR}/*.h)
+	FILE(GLOB GET_SOURCE_FILES_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} ${DIR}/*.tpp)
+	FILE(GLOB GET_SOURCE_FILES_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} ${DIR}/*.cpp)
 
 	SET(${SOURCE_FILES} ${GET_SOURCE_FILES_HEADERS} ${GET_SOURCE_FILES_SOURCES})
 
@@ -54,4 +55,4 @@ MACRO(ADD_GOOGLE_TESTS executable)
 			# message ("Adding test: ${test_name}")
 		ENDFOREACH(hit)
 	ENDFOREACH()
-ENDMACRO()
+ENDMACRO()
\ No newline at end of file

scripts/cmake/cmake/FindLIS.cmake

+# - Try to find LIS
+# Once done, this will define
+#
+#  LIS_FOUND
+#  LIS_INCLUDE_DIRS
+#  LIS_LIBRARIES
+
+set(LIS_ROOT_DIR
+    "${LIS_ROOT_DIR}"
+    CACHE
+    PATH
+    "Directory to search for Lis library")
+	
+find_path( LIS_INCLUDE_DIR
+	NAMES lis.h
+	HINTS
+	${LIS_ROOT_DIR}/include
+    /usr/include/lis
+    $ENV{HOME}/include/
+	)
+
+find_library(LIS_LIBRARY
+    NAMES lis
+    HINTS 
+    ${LIS_ROOT_DIR}/lib
+    /usr/lib
+    $ENV{HOME}/lib/
+    ) 
+
+SET(LIS_LIBRARIES ${LIS_LIBRARY})
+
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(LIS DEFAULT_MSG LIS_LIBRARY LIS_INCLUDE_DIR)
+