diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp index 0c77a338ff5209fbbf2da1d868f59077fd74bf0d..109f49c0ab9f6250638624fda8bd750fdfce94d5 100644 --- a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp +++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp @@ -30,32 +30,7 @@ EigenLisLinearSolver::EigenLisLinearSolver(EigenMatrix &A, : _A(A) { if (option) - setOption(*option); -} - -void EigenLisLinearSolver::setOption(BaseLib::ConfigTree const& option) -{ - boost::optional<BaseLib::ConfigTree> 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; - } + _option.addOptions(*option); } void EigenLisLinearSolver::solve(EigenVector &b_, EigenVector &x_) diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h index 09eae598691394dad20def466d7866b67e4c7eff..1e397702dc2bba0d5fb6cf12da2068fcf0811f58 100644 --- a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h +++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h @@ -40,21 +40,11 @@ public: EigenLisLinearSolver(EigenMatrix &A, const std::string solver_name = "", BaseLib::ConfigTree const*const option = nullptr); - /** - * parse linear solvers configuration - */ - void setOption(BaseLib::ConfigTree const& option); - /** * copy linear solvers options */ void setOption(const LisOption &option) { _option = option; } - /** - * get linear solver options - */ - LisOption &getOption() { return _option; } - /** * solve a given linear equations * diff --git a/MathLib/LinAlg/Lis/LisLinearSolver.cpp b/MathLib/LinAlg/Lis/LisLinearSolver.cpp index 16af3b7d17e63a569f9a7ef131848b223be16945..11ea75bff755e6106304d0297bab18605448491d 100644 --- a/MathLib/LinAlg/Lis/LisLinearSolver.cpp +++ b/MathLib/LinAlg/Lis/LisLinearSolver.cpp @@ -12,14 +12,13 @@ * */ -#include "LisLinearSolver.h" - -#include <string> -#include <sstream> #ifdef _OPENMP #include <omp.h> #endif -#include "logog/include/logog.hpp" + +#include "LisLinearSolver.h" + +#include <logog/include/logog.hpp> #include "LisCheck.h" @@ -32,35 +31,12 @@ LisLinearSolver::LisLinearSolver(LisMatrix &A, : _A(A) { if (option) - setOption(*option); -} + { + _option.addOptions(*option); -void LisLinearSolver::setOption(BaseLib::ConfigTree const& option) -{ - boost::optional<BaseLib::ConfigTree> 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; +#ifndef NDEBUG + _option.printInfo(); +#endif } } @@ -70,62 +46,67 @@ void LisLinearSolver::solve(LisVector &b, LisVector &x) 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 - - // 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; int ierr = lis_solver_create(&solver); checkLisError(ierr); - ierr = lis_solver_set_option(const_cast<char*>(solver_options.c_str()), solver); - checkLisError(ierr); - ierr = lis_solver_set_option(const_cast<char*>(tol_option.c_str()), solver); - checkLisError(ierr); - ierr = lis_solver_set_option(const_cast<char*>("-print mem"), solver); - checkLisError(ierr); - ierr = lis_solver_set_optionC(solver); - checkLisError(ierr); + + { + std::string opt; + for (auto const& it : _option.settings) + { + opt = it.first + " " + it.second; + + ierr = lis_solver_set_option(const_cast<char*>(opt.c_str()), solver); + checkLisError(ierr); + } + } +#ifdef _OPENMP + INFO("-> number of threads: %i", (int) omp_get_max_threads()); +#endif + { + int precon; + ierr = lis_solver_get_precon(solver, &precon); + INFO("-> precon: %i", precon); + } + { + int slv; + ierr = lis_solver_get_solver(solver, &slv); + INFO("-> solver: %i", slv); + } // solve INFO("-> solve"); ierr = lis_solve(_A.getRawMatrix(), b.getRawVector(), x.getRawVector(), solver); checkLisError(ierr); - int iter = 0; - double resid = 0.0; - ierr = lis_solver_get_iter(solver, &iter); - checkLisError(ierr); - ierr = lis_solver_get_residualnorm(solver, &resid); - checkLisError(ierr); - INFO("\t iteration: %d/%ld\n", iter, _option.max_iterations); - INFO("\t residual: %e\n", resid); + { + int iter = 0; + ierr = lis_solver_get_iter(solver, &iter); + checkLisError(ierr); + + std::string max_iter = _option.settings["-maxiter"]; + if (max_iter.empty()) max_iter = "--"; + INFO("-> iteration: %d/%s", iter, max_iter.c_str()); + } + { + double resid = 0.0; + ierr = lis_solver_get_residualnorm(solver, &resid); + checkLisError(ierr); + INFO("-> residual: %g", resid); + } + { + double time, itime, ptime, p_ctime, p_itime; + ierr = lis_solver_get_timeex(solver, &time, &itime, + &ptime, &p_ctime, &p_itime); + checkLisError(ierr); + INFO("-> time total (s): %g", time); + INFO("-> time iterations (s): %g", itime); + INFO("-> time preconditioning (s): %g", ptime); + INFO("-> time precond. create (s): %g", p_ctime); + INFO("-> time precond. iter (s): %g", p_itime); + } // Clear solver ierr = lis_solver_destroy(solver); diff --git a/MathLib/LinAlg/Lis/LisLinearSolver.h b/MathLib/LinAlg/Lis/LisLinearSolver.h index e1fd67c68201c618e2b8fd421f457c407d816913..9e3d0113b12a74ae1fef744d30bd245cb51f041c 100644 --- a/MathLib/LinAlg/Lis/LisLinearSolver.h +++ b/MathLib/LinAlg/Lis/LisLinearSolver.h @@ -18,7 +18,8 @@ #include <vector> #include <string> -#include "lis.h" +#include <lis.h> + #include "BaseLib/ConfigTree.h" #include "LisOption.h" @@ -32,7 +33,7 @@ namespace MathLib * \brief Linear solver using Lis (http://www.ssisc.org/lis/) * */ -class LisLinearSolver +class LisLinearSolver final { public: /** @@ -47,26 +48,12 @@ public: LisLinearSolver(LisMatrix &A, const std::string solver_name = "", BaseLib::ConfigTree const*const option = nullptr); - virtual ~LisLinearSolver() {} - - /** - * configure linear solvers - * @param option - */ - void setOption(BaseLib::ConfigTree const& 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 * diff --git a/MathLib/LinAlg/Lis/LisOption.cpp b/MathLib/LinAlg/Lis/LisOption.cpp index abaaa53d4b3f546cc788c25099feafd31d975f1e..3c59fd2e1c25c6ad8c68deef103b301a9be98aa6 100644 --- a/MathLib/LinAlg/Lis/LisOption.cpp +++ b/MathLib/LinAlg/Lis/LisOption.cpp @@ -14,89 +14,54 @@ #include "LisOption.h" -namespace MathLib -{ - -LisOption::LisOption() -{ - solver_type = SolverType::CG; - precon_type = PreconType::NONE; - matrix_type = MatrixType::CRS; - max_iterations = 1e6; - error_tolerance = 1.e-16; -} +#include <logog/include/logog.hpp> -LisOption::SolverType LisOption::getSolverType(const std::string &solver_name) +namespace MathLib { -#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) +void LisOption::addOptions(const BaseLib::ConfigTree & options) { -#define RETURN_PRECOM_ENUM_IF_SAME_STRING(str, TypeName) \ - if (#TypeName==str) return PreconType::TypeName; + auto const range = options.equal_range("lis_option"); + for (auto it=range.first; it!=range.second; ++it) + { + auto key = it->second.get_optional<Key>("option"); + auto val = it->second.get_optional<Value>("value"); - 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); + if (key && val) { + settings[*key] = *val; + } + } - return PreconType::NONE; -#undef RETURN_PRECOM_ENUM_IF_SAME_STRING + auto solver_type = options.get_optional<std::string>("solver_type"); + if (solver_type) { + settings["-i"] = *solver_type; + } + auto precon_type = options.get_optional<std::string>("precon_type"); + if (precon_type) { + settings["-p"] = *precon_type; + } + auto error_tolerance = options.get_optional<std::string>("error_tolerance"); + if (error_tolerance) { + settings["-tol"] = *error_tolerance; + } + auto max_iteration_step = options.get_optional<std::string>("max_iteration_step"); + if (max_iteration_step) { + settings["-maxiter"] = *max_iteration_step; + } } -LisOption::MatrixType LisOption::getMatrixType(const std::string &matrix_name) +void LisOption::printInfo() const { -#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); + if (settings.empty()) { + INFO("Lis options: none set."); + } else { + INFO("Lis options:"); - return MatrixType::CRS; -#undef RETURN_MATRIX_ENUM_IF_SAME_STRING + for (auto const& it : settings) { + INFO("-> %s %s ", it.first.c_str(), it.second.c_str()); + } + INFO(""); + } } } //MathLib diff --git a/MathLib/LinAlg/Lis/LisOption.h b/MathLib/LinAlg/Lis/LisOption.h index f880ae95afde8084612dfb2960e01f7e3d466545..6ae0eea5f8ef85380545b02db835d6d31e030df9 100644 --- a/MathLib/LinAlg/Lis/LisOption.h +++ b/MathLib/LinAlg/Lis/LisOption.h @@ -16,6 +16,9 @@ #define LIS_OPTION_H_ #include <string> +#include <map> + +#include "BaseLib/ConfigTree.h" namespace MathLib { @@ -25,48 +28,60 @@ namespace MathLib */ 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 - }; + using Key = std::string; + using Value = std::string; - /// Preconditioner type - enum class PreconType : int + std::map<Key, Value> settings; + + LisOption() { - NONE = 0, - JACOBI = 1, - ILU = 2, - SSOR = 3, - Hybrid = 4, - IplusS = 5, - SAINV = 6, - SAAMG = 7, - CroutILU = 8, - ILUT = 9 - }; + // by default do not zero out the solution vector + settings["-initxzeros"] = "0"; + } + + /** + * Adds options from a \c ConfigTree. + * + * Options are stored as strings as obtained from the ConfigTree. + * In particular it will not be checked by this method if options are valid. + * + * It is only guaranteed that each given option will be passed only once to Lis. + * + * The \c ConfigTree section must have the layout as shown in the + * following example: + * + * \code{.xml} + * <linear_solver> + * <solver_type> gmres </solver_type> + * <precon_type> ilu </precon_type> + * <error_tolerance> 1.0e-10 </error_tolerance> + * <max_iteration_step> 200 </max_iteration_step> + * + * <lis_option><option> -ilu_fill </option><value> 3 </value></lis_option> + * <lis_option><option> -print </option><value> mem </value></lis_option> + * </linear_solver> + * \endcode + * + * The first four tags are special in the sense that those options + * will take precedence even if there is an equivalent Lis option specified. + * + * Any possible Lis option can be supplied by a line like + * + * \code{.xml} + * <lis_option><option> -ilu_fill </option><value> 3 </value></lis_option> + * \endcode + * + * For possible values, please refer to the Lis User Guide: + * http://www.ssisc.org/lis/lis-ug-en.pdf + * + * Note: Option -omp_num_threads cannot be used with this class since Lis + * currently (version 1.5.57) only sets the number of threads in + * \c lis_initialize(). Refer to the Lis source code for details. + */ + void addOptions(BaseLib::ConfigTree const& options); + + void printInfo() const; + /// Matrix type enum class MatrixType : int @@ -83,62 +98,6 @@ struct LisOption 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); }; } diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h index cdad7033f53bf70b69f2eccf7a0a2e772300e5b8..39623c1b88c78a2a994c53e9e861c12345a1f09b 100644 --- a/ProcessLib/GroundwaterFlowProcess.h +++ b/ProcessLib/GroundwaterFlowProcess.h @@ -129,6 +129,16 @@ public: std::abort(); } } + + // Linear solver options + { + auto const par = config.get_child_optional("linear_solver"); + + if (par) + { + _linear_solver_options.reset(new BaseLib::ConfigTree(*par)); + } + } } template <unsigned GlobalDim> @@ -237,7 +247,8 @@ public: _x.reset(_global_setup.createVector(num_unknowns)); _rhs.reset(_global_setup.createVector(num_unknowns)); - _linearSolver.reset(new typename GlobalSetup::LinearSolver(*_A, "gw_")); + _linear_solver.reset(new typename GlobalSetup::LinearSolver( + *_A, "gw_", _linear_solver_options.get())); setInitialConditions(*_hydraulic_head); @@ -304,7 +315,7 @@ public: MathLib::applyKnownSolution(*_A, *_rhs, _dirichlet_bc.global_ids, _dirichlet_bc.values); #endif - _linearSolver->solve(*_rhs, *_x); + _linear_solver->solve(*_rhs, *_x); return true; } @@ -389,7 +400,8 @@ private: std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets; GlobalSetup _global_setup; - std::unique_ptr<typename GlobalSetup::LinearSolver> _linearSolver; + std::unique_ptr<BaseLib::ConfigTree> _linear_solver_options; + std::unique_ptr<typename GlobalSetup::LinearSolver> _linear_solver; std::unique_ptr<typename GlobalSetup::MatrixType> _A; std::unique_ptr<typename GlobalSetup::VectorType> _rhs; std::unique_ptr<typename GlobalSetup::VectorType> _x;