From 67437f72bf82814e86cda1533a29cbfedbeb58ee Mon Sep 17 00:00:00 2001 From: Norihiro Watanabe <norihiro.watanabe@ufz.de> Date: Tue, 4 Dec 2012 20:17:14 +0100 Subject: [PATCH] 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 | 6 + .../LisLinearSystem.cpp | 193 +++++++++++++++++ .../SystemOfLinearEquations/LisLinearSystem.h | 194 ++++++++++++++++++ .../SystemOfLinearEquations/LisOption.cpp | 99 +++++++++ .../SystemOfLinearEquations/LisOption.h | 143 +++++++++++++ 5 files changed, 635 insertions(+) create mode 100644 MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp create mode 100644 MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h create mode 100644 MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp create mode 100644 MathLib/LinAlg/SystemOfLinearEquations/LisOption.h diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt index ee8412af6b9..4d2707d271f 100644 --- a/MathLib/CMakeLists.txt +++ b/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) diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp new file mode 100644 index 00000000000..fa51527a947 --- /dev/null +++ b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp @@ -0,0 +1,193 @@ +/** + * 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 diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h new file mode 100644 index 00000000000..e6493d8bc9d --- /dev/null +++ b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h @@ -0,0 +1,194 @@ +/** + * 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_ + diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp b/MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp new file mode 100644 index 00000000000..68ba72317a7 --- /dev/null +++ b/MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp @@ -0,0 +1,99 @@ +/** + * 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 diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisOption.h b/MathLib/LinAlg/SystemOfLinearEquations/LisOption.h new file mode 100644 index 00000000000..698bc5f13a5 --- /dev/null +++ b/MathLib/LinAlg/SystemOfLinearEquations/LisOption.h @@ -0,0 +1,143 @@ +/** + * 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_ -- GitLab