Merge pull request #889 from chleh/lis-options

Rewrite of LisOption

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_)

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
      *

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);

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
      *

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

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);
 };
 
 }

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;