add LisLinearSystem class

updating LisOption
- documentation
- add more solver and preconditioner type
- add matrix_type

update LisLinearSysetem
- impose prescribed values by penalty coefficient
- support matrix_type, solver_precon_arg in LisOption

MathLib/CMakeLists.txt

@@ -21,6 +21,12 @@ 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/LisLinearSystem.h
+        LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp
+    )
+ENDIF ()
 SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_LEQS})
 
 IF (METIS_FOUND)

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"
+
+#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, RowMajorSparsity* /*sp*/)
+: _dim(dimension), _is_initialized(false), _max_diag_coeff(0)
+{
+    setZero();
+}
+
+LisLinearSystem::~LisLinearSystem()
+{
+    lis_matrix_destroy(_AA);
+    lis_vector_destroy(_bb);
+    lis_vector_destroy(_xx);
+}
+
+void LisLinearSystem::setOption(const boost::property_tree::ptree &option)
+{
+    using boost::property_tree::ptree;
+
+    boost::optional<ptree> ptSolver = option.get_child("LinearSolver");
+    if (!ptSolver.is_initialized())
+        return;
+
+    ptree op = ptSolver.get();
+
+    boost::optional<std::string> solver_type = op.get_optional<std::string>("solver_type");
+    if (solver_type.is_initialized()) {
+        _option.solver_type = _option.getSolverType(solver_type.get());
+    }
+    boost::optional<std::string> precon_type = op.get_optional<std::string>("precon_type");
+    if (precon_type.is_initialized()) {
+        _option.precon_type = _option.getPreconType(precon_type.get());
+    }
+    boost::optional<std::string> matrix_type = op.get_optional<std::string>("matrix_type");
+    if (matrix_type) {
+        _option.matrix_type = _option.getMatrixType(matrix_type.get());
+    }
+    boost::optional<double> error_tolerance = op.get_optional<double>("error_tolerance");
+    if (error_tolerance.is_initialized()) {
+        _option.error_tolerance = error_tolerance.get();
+    }
+    boost::optional<int> max_iteration_step = op.get_optional<int>("max_iteration_step");
+    if (max_iteration_step.is_initialized()) {
+        _option.max_iterations = max_iteration_step.get();
+    }
+}
+
+void LisLinearSystem::setZero()
+{
+    int ierr = 0;
+    // create a matrix
+    if (_is_initialized) {
+        // 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);
+    }
+#ifndef USE_MPI
+    ierr = lis_matrix_create(0, &_AA);
+    ierr = lis_matrix_set_size(_AA, 0, getDimension());
+#else
+    lis_matrix_create(MPI_COMM_WORLD, &_AA);
+    lis_matrix_set_size(_AA, dimension, 0);
+//    lis_matrix_get_range(A,&is,&ie);
+#endif
+    checkError(ierr);
+
+    // crate or zero RHS, x
+    if (!_is_initialized) {
+        ierr = lis_vector_duplicate(_AA, &_bb);
+        ierr = lis_vector_duplicate(_AA, &_xx);
+        _is_initialized = true;
+    } else {
+        ierr = lis_vector_set_all(0.0, _bb);
+        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) = epsilon
+        setMatEntry(rowId, rowId, penalty);
+        //b(k) = x*epsilon
+        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());
+    const int nthreads = omp_get_max_threads();
+#else
+    const int nthreads = 1;
+#endif
+
+    int ierr = 0;
+
+    applyKnownSolution();
+    // assemble a matrix
+    ierr = lis_matrix_set_type(_AA, static_cast<int>(_option.matrix_type));
+    ierr = lis_matrix_assemble(_AA);
+    checkError(ierr);
+
+    // configure option
+    const std::size_t MAX_ZEILE = 512;
+    char solver_options[MAX_ZEILE], tol_option[MAX_ZEILE];
+    if (_option.solver_precon_arg.length()==0) {
+        sprintf(solver_options, "-i %d -p %d %s ", static_cast<int>(_option.solver_type), static_cast<int>(_option.precon_type), _option.extra_arg.c_str());
+    } else {
+        sprintf(solver_options, "%s ", _option.solver_precon_arg.c_str());
+    }
+    sprintf(tol_option, "-tol %e -maxiter %ld -omp_num_threads %d -initx_zeros 0", _option.error_tolerance, _option.max_iterations, nthreads);
+
+    // Create solver
+    LIS_SOLVER solver;
+    ierr = lis_solver_create(&solver);
+    ierr = lis_solver_set_option(solver_options, solver);
+    ierr = lis_solver_set_option(tol_option, solver);
+    ierr = lis_solver_set_option((char*)"-print mem", solver);
+    
+    // solve
+    INFO("-> solve");
+    ierr = lis_solve(_AA, _bb, _xx, solver);
+    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);
+    ierr = lis_solver_get_residualnorm(solver, &resid);
+    printf("\t iteration: %d/%ld\n", iter, _option.max_iterations);
+    printf("\t residuals: %e\n", resid);
+
+    // Clear solver
+    lis_solver_destroy(solver);
+    std::cout << "------------------------------------------------------------------" << std::endl;
+}
+
+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 "MathLib/LinAlg/Sparse/Sparsity.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
+     * @param sp        Sparsity
+     */
+    LisLinearSystem(std::size_t length, RowMajorSparsity* sp);
+
+    /**
+     *
+     */
+    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();
+
+    /// get entry in A
+    virtual double getMatEntry(std::size_t rowId, std::size_t colId) const
+    {
+        //TODO
+        return .0;
+    }
+
+    /// 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:
+    std::size_t _dim;
+    bool _is_initialized;
+    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 (str.compare(#TypeName)==0) 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 (str.compare(#TypeName)==0) 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 (str.compare(#TypeName)==0) 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, VBR);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, DNS);
+    RETURN_MATRIX_ENUM_IF_SAME_STRING(matrix_name, COO);
+
+    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>
+#include "lis.h"
+
+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 = LIS_MATRIX_CRS,
+        CCS = LIS_MATRIX_CCS,
+        MSR = LIS_MATRIX_MSR,
+        DIA = LIS_MATRIX_DIA,
+        ELL = LIS_MATRIX_ELL,
+        JDS = LIS_MATRIX_JDS,
+        BSR = LIS_MATRIX_BSR,
+        VBR = LIS_MATRIX_VBR,
+        DNS = LIS_MATRIX_DNS,
+        COO = LIS_MATRIX_COO
+    };
+
+    /// 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_