diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index c0d598e..ca362d3 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -69,7 +69,7 @@ jobs: cd build ctest --verbose --output-on-failure lcov --directory . --capture --output-file coverage.info - lcov --remove coverage.info '/usr/*' "${HOME}"'/.cache/*' '*tests/*.cpp' --output-file coverage.info + lcov --remove coverage.info '/usr/*' "${HOME}"'/.cache/*' '*tests/*.cpp' '*tests/*.h' --output-file coverage.info - name: Upload Coverage uses: codecov/codecov-action@v3 diff --git a/CMakeLists.txt b/CMakeLists.txt index f9c5bdd..3c9c53c 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -108,7 +108,7 @@ endif() # Set default minimum C++ standard if(POLYSOLVE_TOPLEVEL_PROJECT) - set(CMAKE_CXX_STANDARD 14) + set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) endif() @@ -145,14 +145,25 @@ endif() ################################################################################ # Add an empty library and fill in the list of sources in `src/CMakeLists.txt`. +add_library(polysolve_linear) +add_library(polysolve::linear ALIAS polysolve_linear) +add_subdirectory(src/polysolve) +add_subdirectory(src/polysolve/linear) + add_library(polysolve) add_library(polysolve::polysolve ALIAS polysolve) +add_subdirectory(src/polysolve/nonlinear) +target_link_libraries(polysolve_linear PUBLIC polysolve_coverage_config) target_link_libraries(polysolve PUBLIC polysolve_coverage_config) -add_subdirectory(src/polysolve) + + +target_compile_features(polysolve_linear PUBLIC cxx_std_17) +target_compile_features(polysolve PUBLIC cxx_std_17) # Public include directory for Polysolve +target_include_directories(polysolve_linear PUBLIC ${PROJECT_SOURCE_DIR}/src) target_include_directories(polysolve PUBLIC ${PROJECT_SOURCE_DIR}/src) ################################################################################ @@ -160,13 +171,39 @@ target_include_directories(polysolve PUBLIC ${PROJECT_SOURCE_DIR}/src) ################################################################################ if(POLYSOLVE_LARGE_INDEX) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_LARGE_INDEX) + target_compile_definitions(polysolve_linear PUBLIC -DPOLYSOLVE_LARGE_INDEX) endif() +target_compile_definitions(polysolve_linear PUBLIC -DPOLYSOLVE_LINEAR_SPEC="${CMAKE_SOURCE_DIR}/linear-solver-spec.json") +target_compile_definitions(polysolve_linear PUBLIC -DPOLYSOLVE_NON_LINEAR_SPEC="${CMAKE_SOURCE_DIR}/non-linear-solver-spec.json") + + ################################################################################ # Dependencies ################################################################################ +# polysolve_linear +target_link_libraries(polysolve PUBLIC polysolve::linear) + +# CppNumericalSolvers +include(cppoptlib) +target_link_libraries(polysolve PUBLIC cppoptlib) + +# LBFGSpp +include(LBFGSpp) +target_link_libraries(polysolve PUBLIC LBFGSpp::LBFGSpp) + + + + +# JSON Specification Engine library +include(jse) +target_link_libraries(polysolve_linear PUBLIC jse::jse) + +# spdlog +include(spdlog) +target_link_libraries(polysolve_linear PUBLIC spdlog::spdlog) + # Accelerate solver if(POLYSOLVE_WITH_ACCELERATE) include(CPM) @@ -179,61 +216,63 @@ if(POLYSOLVE_WITH_ACCELERATE) set(BLA_VENDOR Apple) find_package(BLAS REQUIRED) find_package(LAPACK REQUIRED) - target_link_libraries(polysolve PRIVATE BLAS::BLAS LAPACK::LAPACK) + target_link_libraries(polysolve_linear PRIVATE BLAS::BLAS LAPACK::LAPACK) endif() # Extra warnings include(polysolve_warnings) +target_link_libraries(polysolve_linear PRIVATE polysolve::warnings) target_link_libraries(polysolve PRIVATE polysolve::warnings) # Sanitizers if(POLYSOLVE_WITH_SANITIZERS) include(sanitizers) + add_sanitizers(polysolve_linear) add_sanitizers(polysolve) endif() include(eigen) -target_link_libraries(polysolve PUBLIC Eigen3::Eigen) +target_link_libraries(polysolve_linear PUBLIC Eigen3::Eigen) # Hypre (GNU Lesser General Public License) if(POLYSOLVE_WITH_HYPRE) include(hypre) - target_link_libraries(polysolve PUBLIC HYPRE::HYPRE) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_HYPRE) + target_link_libraries(polysolve_linear PUBLIC HYPRE::HYPRE) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_HYPRE) if(HYPRE_WITH_MPI) - target_compile_definitions(polysolve PUBLIC HYPRE_WITH_MPI) + target_compile_definitions(polysolve_linear PUBLIC HYPRE_WITH_MPI) endif() endif() # Json (MIT) include(json) -target_link_libraries(polysolve PUBLIC nlohmann_json::nlohmann_json) +target_link_libraries(polysolve_linear PUBLIC nlohmann_json::nlohmann_json) # Accelerate solvers if(POLYSOLVE_WITH_ACCELERATE) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_ACCELERATE) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_ACCELERATE) endif() # CHOLMOD solver if(POLYSOLVE_WITH_CHOLMOD) include(suitesparse) - target_link_libraries(polysolve PRIVATE SuiteSparse::CHOLMOD) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CHOLMOD) + target_link_libraries(polysolve_linear PRIVATE SuiteSparse::CHOLMOD) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_CHOLMOD) endif() # MKL library if(POLYSOLVE_WITH_MKL) include(mkl) - target_link_libraries(polysolve PRIVATE mkl::mkl) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_MKL) + target_link_libraries(polysolve_linear PRIVATE mkl::mkl) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_MKL) endif() # Pardiso solver if(POLYSOLVE_WITH_PARDISO) include(pardiso) if(TARGET Pardiso::Pardiso) - target_link_libraries(polysolve PRIVATE Pardiso::Pardiso) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PARDISO) + target_link_libraries(polysolve_linear PRIVATE Pardiso::Pardiso) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_PARDISO) else() message(WARNING "Pardiso not found, solver will not be available.") endif() @@ -242,16 +281,16 @@ endif() # UmfPack solver if(POLYSOLVE_WITH_UMFPACK) include(suitesparse) - target_link_libraries(polysolve PRIVATE SuiteSparse::UMFPACK) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_UMFPACK) + target_link_libraries(polysolve_linear PRIVATE SuiteSparse::UMFPACK) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_UMFPACK) endif() # SuperLU solver if(POLYSOLVE_WITH_SUPERLU) include(superlu) if(TARGET SuperLU::SuperLU) - target_link_libraries(polysolve PRIVATE SuperLU::SuperLU) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_SUPERLU) + target_link_libraries(polysolve_linear PRIVATE SuperLU::SuperLU) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_SUPERLU) else() message(WARNING "SuperLU not found, solver will not be available.") endif() @@ -260,23 +299,23 @@ endif() # AMGCL solver if(POLYSOLVE_WITH_AMGCL) include(amgcl) - target_link_libraries(polysolve PUBLIC amgcl::amgcl) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_AMGCL) + target_link_libraries(polysolve_linear PUBLIC amgcl::amgcl) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_AMGCL) endif() # Spectra (MPL 2.0) if(POLYSOLVE_WITH_SPECTRA) include(spectra) - target_link_libraries(polysolve PUBLIC Spectra::Spectra) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_SPECTRA) + target_link_libraries(polysolve_linear PUBLIC Spectra::Spectra) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_SPECTRA) endif() # cuSolver solvers if(POLYSOLVE_WITH_CUSOLVER) include(cusolverdn) if(TARGET CUDA::cusolver) - target_link_libraries(polysolve PUBLIC CUDA::cusolver) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUSOLVER) + target_link_libraries(polysolve_linear PUBLIC CUDA::cusolver) + target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_CUSOLVER) else() message(WARNING "cuSOLVER not found, solver will not be available.") endif() @@ -286,8 +325,6 @@ endif() # Compiler options ################################################################################ -# Use C++14 -target_compile_features(polysolve PUBLIC cxx_std_14) ################################################################################ # Tests diff --git a/README.md b/README.md index f1e17e0..231ca18 100644 --- a/README.md +++ b/README.md @@ -15,10 +15,10 @@ This library contains a cross-platform Eigen wrapper for many different external ```c++ const std::string solver_name = "Hypre" -auto solver = LinearSolver::create(solver_name, ""); +auto solver = Solver::create(solver_name, ""); // Configuration parameters like iteration or accuracy for iterative solvers -// solver->setParameters(params); +// solver->set_parameters(params); // System sparse matrix Eigen::SparseMatrix A; @@ -29,12 +29,12 @@ Eigen::VectorXd b; // Solution Eigen::VectorXd x(b.size()); -solver->analyzePattern(A, A.rows()); +solver->analyze_pattern(A, A.rows()); solver->factorize(A); solver->solve(b, x); ``` -You can use `LinearSolver::availableSolvers()` to obtain the list of available solvers. +You can use `Solver::available_solvers()` to obtain the list of available solvers. ## Parameters diff --git a/cmake/recipes/LBFGSpp.cmake b/cmake/recipes/LBFGSpp.cmake new file mode 100644 index 0000000..680a377 --- /dev/null +++ b/cmake/recipes/LBFGSpp.cmake @@ -0,0 +1,20 @@ +# LBFGSpp (https://github.com/yixuan/LBFGSpp) +# License: MIT + +if(TARGET LBGFSpp::LBGFSpp) + return() +endif() + +message(STATUS "Third-party: creating target 'LBGFSpp::LBGFSpp'") + +include(CPM) +CPMAddPackage( + NAME lbfgspp + GITHUB_REPOSITORY yixuan/LBFGSpp + GIT_TAG v0.2.0 + DOWNLOAD_ONLY TRUE +) + +add_library(LBFGSpp INTERFACE) +target_include_directories(LBFGSpp INTERFACE "${lbfgspp_SOURCE_DIR}/include") +add_library(LBFGSpp::LBFGSpp ALIAS LBFGSpp) diff --git a/cmake/recipes/cppoptlib.cmake b/cmake/recipes/cppoptlib.cmake new file mode 100644 index 0000000..e19709e --- /dev/null +++ b/cmake/recipes/cppoptlib.cmake @@ -0,0 +1,19 @@ +# CppNumericalSolvers (https://github.com/PatWie/CppNumericalSolvers.git) +# License: MIT + +if(TARGET cppoptlib) + return() +endif() + +message(STATUS "Third-party: creating target 'cppoptlib'") + +include(CPM) +CPMAddPackage( + NAME cppoptlib + GITHUB_REPOSITORY PatWie/CppNumericalSolvers + GIT_TAG 7eddf28fa5a8872a956d3c8666055cac2f5a535d + DOWNLOAD_ONLY TRUE +) + +add_library(cppoptlib INTERFACE) +target_include_directories(cppoptlib SYSTEM INTERFACE ${cppoptlib_SOURCE_DIR}/include) \ No newline at end of file diff --git a/cmake/recipes/jse.cmake b/cmake/recipes/jse.cmake new file mode 100644 index 0000000..4b67c02 --- /dev/null +++ b/cmake/recipes/jse.cmake @@ -0,0 +1,11 @@ +# json spec engine (https://github.com/geometryprocessing/json-spec-engine) +# License: MIT + +if(TARGET jse::jse) + return() +endif() + +message(STATUS "Third-party: creating target 'jse::jse'") + +include(CPM) +CPMAddPackage("gh:geometryprocessing/json-spec-engine#1261dc89478c7646ff99cbed8bc5357c2813565d") \ No newline at end of file diff --git a/cmake/recipes/spdlog.cmake b/cmake/recipes/spdlog.cmake new file mode 100644 index 0000000..a61cd2a --- /dev/null +++ b/cmake/recipes/spdlog.cmake @@ -0,0 +1,37 @@ +# +# Copyright 2020 Adobe. All rights reserved. +# This file is licensed to you under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. You may obtain a copy +# of the License at http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under +# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS +# OF ANY KIND, either express or implied. See the License for the specific language +# governing permissions and limitations under the License. +# + +# spdlog (https://github.com/gabime/spdlog) +# License: MIT + +if(TARGET spdlog::spdlog) + return() +endif() + +message(STATUS "Third-party: creating target 'spdlog::spdlog'") + +option(SPDLOG_INSTALL "Generate the install target" ON) +set(CMAKE_INSTALL_DEFAULT_COMPONENT_NAME "spdlog") + +include(CPM) +CPMAddPackage("gh:gabime/spdlog@1.9.2") + +set_target_properties(spdlog PROPERTIES POSITION_INDEPENDENT_CODE ON) + +set_target_properties(spdlog PROPERTIES FOLDER external) + +if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang" OR + "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + target_compile_options(spdlog PRIVATE + "-Wno-sign-conversion" + ) +endif() \ No newline at end of file diff --git a/linear-solver-spec.json b/linear-solver-spec.json new file mode 100644 index 0000000..1e59528 --- /dev/null +++ b/linear-solver-spec.json @@ -0,0 +1,425 @@ +[ + { + "pointer": "/", + "default": null, + "type": "object", + "optional": [ + "enable_overwrite_solver", + "solver", + "precond", + "Eigen::LeastSquaresConjugateGradient", + "Eigen::DGMRES", + "Eigen::ConjugateGradient", + "Eigen::BiCGSTAB", + "Eigen::GMRES", + "Eigen::MINRES", + "Pardiso", + "Hypre", + "AMGCL" + ], + "doc": "Settings for the linear solver." + }, + { + "pointer": "/enable_overwrite_solver", + "default": false, + "type": "bool", + "doc": "If solver name is not present, falls back to default" + }, + { + "pointer": "/solver", + "default": "", + "type": "string", + "doc": "Linear solver type.", + "options": [ + "Eigen::SimplicialLDLT", + "Eigen::SparseLU", + "Eigen::CholmodSupernodalLLT", + "Eigen::UmfPackLU", + "Eigen::SuperLU", + "Eigen::PardisoLDLT", + "Eigen::PardisoLLT", + "Eigen::PardisoLU", + "Pardiso", + "Hypre", + "AMGCL", + "Eigen::LeastSquaresConjugateGradient", + "Eigen::DGMRES", + "Eigen::ConjugateGradient", + "Eigen::BiCGSTAB", + "Eigen::GMRES", + "Eigen::MINRES" + ] + }, + { + "pointer": "/precond", + "default": "", + "type": "string", + "doc": "Preconditioner used if using an iterative linear solver.", + "options": [ + "Eigen::IdentityPreconditioner", + "Eigen::DiagonalPreconditioner", + "Eigen::IncompleteCholesky", + "Eigen::LeastSquareDiagonalPreconditioner", + "Eigen::IncompleteLUT" + ] + }, + { + "pointer": "/Eigen::LeastSquaresConjugateGradient", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's Least Squares Conjugate Gradient solver." + }, + { + "pointer": "/Eigen::DGMRES", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's DGMRES solver." + }, + { + "pointer": "/Eigen::ConjugateGradient", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's Conjugate Gradient solver." + }, + { + "pointer": "/Eigen::BiCGSTAB", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's BiCGSTAB solver." + }, + { + "pointer": "/Eigen::GMRES", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's GMRES solver." + }, + { + "pointer": "/Eigen::MINRES", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "tolerance" + ], + "doc": "Settings for the Eigen's MINRES solver." + }, + { + "pointer": "/Pardiso", + "default": null, + "type": "object", + "optional": [ + "mtype" + ], + "doc": "Settings for the Pardiso solver." + }, + { + "pointer": "/Hypre", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "pre_max_iter", + "tolerance" + ], + "doc": "Settings for the Hypre solver." + }, + { + "pointer": "/AMGCL", + "default": null, + "type": "object", + "optional": [ + "solver", + "precond" + ], + "doc": "Settings for the AMGCL solver." + }, + { + "pointer": "/Eigen::LeastSquaresConjugateGradient/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::LeastSquaresConjugateGradient/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Eigen::DGMRES/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::DGMRES/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Eigen::ConjugateGradient/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::ConjugateGradient/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Eigen::BiCGSTAB/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::BiCGSTAB/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Eigen::GMRES/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::GMRES/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Eigen::MINRES/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Eigen::MINRES/tolerance", + "default": 1e-12, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/Pardiso/mtype", + "default": 11, + "type": "int", + "options": [ + 1, + 2, + -2, + 3, + 4, + -4, + 6, + 11, + 13 + ], + "doc": "Matrix type." + }, + { + "pointer": "/Hypre/max_iter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/Hypre/pre_max_iter", + "default": 1, + "type": "int", + "doc": "Maximum number of pre iterations." + }, + { + "pointer": "/Hypre/tolerance", + "default": 1e-10, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/AMGCL/solver", + "default": null, + "type": "object", + "optional": [ + "tol", + "maxiter", + "type" + ], + "doc": "Solver settings for the AMGCL." + }, + { + "pointer": "/AMGCL/precond", + "default": null, + "type": "object", + "optional": [ + "relax", + "class", + "max_levels", + "direct_coarse", + "ncycle", + "coarsening" + ], + "doc": "Preconditioner settings for the AMGCL." + }, + { + "pointer": "/AMGCL/solver/maxiter", + "default": 1000, + "type": "int", + "doc": "Maximum number of iterations." + }, + { + "pointer": "/AMGCL/solver/tol", + "default": 1e-10, + "type": "float", + "doc": "Convergence tolerance." + }, + { + "pointer": "/AMGCL/solver/type", + "default": "cg", + "type": "string", + "doc": "Type of solver to use." + }, + { + "pointer": "/AMGCL/precond/relax", + "default": null, + "type": "object", + "optional": [ + "degree", + "type", + "power_iters", + "higher", + "lower", + "scale" + ], + "doc": "Preconditioner settings for the AMGCL." + }, + { + "pointer": "/AMGCL/precond/class", + "default": "amg", + "type": "string", + "doc": "Type of preconditioner to use." + }, + { + "pointer": "/AMGCL/precond/max_levels", + "default": 6, + "type": "int", + "doc": "Maximum number of levels." + }, + { + "pointer": "/AMGCL/precond/direct_coarse", + "default": false, + "type": "bool", + "doc": "Use direct solver for the coarsest level." + }, + { + "pointer": "/AMGCL/precond/ncycle", + "default": 2, + "type": "int", + "doc": "Number of cycles." + }, + { + "pointer": "/AMGCL/precond/coarsening", + "default": null, + "type": "object", + "optional": [ + "type", + "estimate_spectral_radius", + "relax", + "aggr" + ], + "doc": "Coarsening parameters." + }, + { + "pointer": "/AMGCL/precond/relax/degree", + "default": 16, + "type": "int", + "doc": "Degree of the polynomial." + }, + { + "pointer": "/AMGCL/precond/relax/type", + "default": "chebyshev", + "type": "string", + "doc": "Type of relaxation to use." + }, + { + "pointer": "/AMGCL/precond/relax/power_iters", + "default": 100, + "type": "int", + "doc": "Number of power iterations." + }, + { + "pointer": "/AMGCL/precond/relax/higher", + "default": 2, + "type": "float", + "doc": "Higher level relaxation." + }, + { + "pointer": "/AMGCL/precond/relax/lower", + "default": 0.008333333333, + "type": "float", + "doc": "Lower level relaxation." + }, + { + "pointer": "/AMGCL/precond/relax/scale", + "default": true, + "type": "bool", + "doc": "Scale." + }, + { + "pointer": "/AMGCL/precond/coarsening/type", + "default": "smoothed_aggregation", + "type": "string", + "doc": "Coarsening type." + }, + { + "pointer": "/AMGCL/precond/coarsening/estimate_spectral_radius", + "default": true, + "type": "bool", + "doc": "Should the spectral radius be estimated." + }, + { + "pointer": "/AMGCL/precond/coarsening/relax", + "default": 1, + "type": "float", + "doc": "Coarsening relaxation." + }, + { + "pointer": "/AMGCL/precond/coarsening/aggr", + "default": null, + "type": "object", + "optional": [ + "eps_strong" + ], + "doc": "Aggregation settings." + }, + { + "pointer": "/AMGCL/precond/coarsening/aggr/eps_strong", + "default": 0, + "type": "float", + "doc": "Aggregation epsilon strong." + } +] \ No newline at end of file diff --git a/non-linear-solver-spec.json b/non-linear-solver-spec.json new file mode 100644 index 0000000..8c577f8 --- /dev/null +++ b/non-linear-solver-spec.json @@ -0,0 +1,256 @@ +[ + { + "pointer": "/", + "default": null, + "type": "object", + "optional": [ + "solver", + "x_delta", + "grad_norm", + "first_grad_norm_tol", + "max_iterations", + "iterations_per_strategy", + "line_search", + "allow_out_of_iterations", + "LBFGS", + "Newton", + "box_constraints" + ], + "doc": "Settings for nonlinear solver. Interior-loop linear solver settings are defined in the solver/linear section." + }, + { + "pointer": "/solver", + "default": "newton", + "type": "string", + "options": [ + "Newton", + "DenseNewton", + "GradientDescent", + "L-BFGS", + "BFGS", + "L-BFGS-B", + "MMA" + ], + "doc": "Nonlinear solver type" + }, + { + "pointer": "/x_delta", + "default": 0, + "type": "float", + "min": 0, + "doc": "Stopping criterion: minimal change of the variables x for the iterations to continue. Computed as the L2 norm of x divide by the time step." + }, + { + "pointer": "/grad_norm", + "default": 1e-08, + "type": "float", + "min": 0, + "doc": "Stopping criterion: Minimal gradient norm for the iterations to continue." + }, + { + "pointer": "/first_grad_norm_tol", + "default": 1e-10, + "type": "float", + "doc": "Minimal gradient norm for the iterations to not start, assume we already are at a minimum." + }, + { + "pointer": "/max_iterations", + "default": 500, + "type": "int", + "doc": "Maximum number of iterations for a nonlinear solve." + }, + { + "pointer": "/iterations_per_strategy", + "default": 5, + "type": "int", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/iterations_per_strategy", + "type": "list", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/iterations_per_strategy/*", + "default": 5, + "type": "int", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/allow_out_of_iterations", + "default": false, + "type": "bool", + "doc": "If false (default), an exception will be thrown when the nonlinear solver reaches the maximum number of iterations." + }, + { + "pointer": "/LBFGS", + "default": null, + "type": "object", + "optional": [ + "history_size" + ], + "doc": "Options for LBFGS." + }, + { + "pointer": "/LBFGS/history_size", + "default": 6, + "type": "int", + "doc": "The number of corrections to approximate the inverse Hessian matrix." + }, + { + "pointer": "/Newton", + "default": null, + "type": "object", + "optional": [ + "reg_weight_min", + "reg_weight_max", + "reg_weight_inc", + "reg_weight_dec", + "force_psd_projection" + ], + "doc": "Options for Newton." + }, + { + "pointer": "/Newton/reg_weight_min", + "default": 1e-8, + "type": "float", + "doc": "Minimum regulariztion weight." + }, + { + "pointer": "/Newton/reg_weight_max", + "default": 1e8, + "type": "float", + "doc": "Maximum regulariztion weight." + }, + { + "pointer": "/Newton/reg_weight_inc", + "default": 10, + "type": "float", + "doc": "Regulariztion weight increment." + }, + { + "pointer": "/Newton/reg_weight_dec", + "default": 2, + "type": "float", + "doc": "Regulariztion weight decrement." + }, + { + "pointer": "/Newton/force_psd_projection", + "default": false, + "type": "bool", + "doc": "Force the Hessian to be PSD when using second order solvers (i.e., Newton's method)." + }, + { + "pointer": "/line_search", + "default": null, + "type": "object", + "optional": [ + "method", + "use_grad_norm_tol", + "min_step_size", + "max_step_size_iter", + "min_step_size_final", + "max_step_size_iter_final", + "default_init_step_size", + "step_ratio", + "Armijo" + ], + "doc": "Settings for line-search in the nonlinear solver" + }, + { + "pointer": "/line_search/method", + "default": "Backtracking", + "type": "string", + "options": [ + "Armijo", + "ArmijoAlt", + "Backtracking", + "MoreThuente", + "None" + ], + "doc": "Line-search type" + }, + { + "pointer": "/line_search/use_grad_norm_tol", + "default": 1e-6, + "type": "float", + "doc": "When the energy is smaller than use_grad_norm_tol, line-search uses norm of gradient instead of energy" + }, + { + "pointer": "/line_search/min_step_size", + "default": 1e-10, + "type": "float", + "doc": "Mimimum step size" + }, + { + "pointer": "/line_search/max_step_size_iter", + "default": 30, + "type": "int", + "doc": "Number of iterations" + }, + { + "pointer": "/line_search/min_step_size_final", + "default": 1e-20, + "type": "float", + "doc": "Mimimum step size for last descent strategy" + }, + { + "pointer": "/line_search/max_step_size_iter_final", + "default": 100, + "type": "int", + "doc": "Number of iterations for last descent strategy" + }, + { + "pointer": "/line_search/default_init_step_size", + "default": 1, + "type": "float", + "doc": "Initial step size" + }, + { + "pointer": "/line_search/step_ratio", + "default": 0.5, + "type": "float", + "doc": "Ratio used to decrease the step" + }, + { + "pointer": "/line_search/Armijo", + "default": null, + "type": "object", + "optional": [ + "c" + ], + "doc": "Options for Armijo." + }, + { + "pointer": "/line_search/Armijo/c", + "default": 0.5, + "type": "float", + "doc": "Armijo c parameter." + }, + { + "pointer": "/box_constraints", + "type": "object", + "optional": [ + "bounds", + "max_change" + ], + "default": null + }, + { + "pointer": "/box_constraints/bounds", + "default": [], + "type": "list", + "doc": "Box constraints on optimization variables." + }, + { + "pointer": "/box_constraints/bounds/*", + "type": "float", + "doc": "Box constraint values on optimization variables." + }, + { + "pointer": "/box_constraints/max_change", + "default": -1, + "type": "float", + "doc": "Maximum change of optimization variables in one iteration, only for solvers with box constraints. Negative value to disable this constraint." + } +] \ No newline at end of file diff --git a/src/polysolve/CMakeLists.txt b/src/polysolve/CMakeLists.txt index 831fab2..5477687 100644 --- a/src/polysolve/CMakeLists.txt +++ b/src/polysolve/CMakeLists.txt @@ -1,22 +1,8 @@ set(SOURCES - FEMSolver.cpp - FEMSolver.hpp - LbfgsSolver.hpp - LinearSolver.cpp - LinearSolver.hpp - LinearSolverAMGCL.cpp - LinearSolverAMGCL.hpp - LinearSolverCuSolverDN.cu - LinearSolverCuSolverDN.cuh - LinearSolverEigen.hpp - LinearSolverEigen.tpp - LinearSolverHypre.cpp - LinearSolverHypre.hpp - LinearSolverPardiso.cpp - LinearSolverPardiso.hpp - SaddlePointSolver.cpp - SaddlePointSolver.hpp + Utils.hpp + Utils.cpp ) source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES}) -target_sources(polysolve PRIVATE ${SOURCES}) +target_sources(polysolve_linear PRIVATE ${SOURCES}) + diff --git a/src/polysolve/FEMSolver.cpp b/src/polysolve/FEMSolver.cpp deleted file mode 100644 index 12a856e..0000000 --- a/src/polysolve/FEMSolver.cpp +++ /dev/null @@ -1,343 +0,0 @@ -//////////////////////////////////////////////////////////////////////////////// -#include - -#ifdef POLYSOLVE_WITH_SPECTRA -#include -#include -#include -#include -#endif - -#include - -//////////////////////////////////////////////////////////////////////////////// - -namespace polysolve -{ - namespace - { - Eigen::Vector4d compute_spectrum(const StiffnessMatrix &mat) - { -#ifdef POLYSOLVE_WITH_SPECTRA - typedef Spectra::SparseSymMatProd MatOp; - typedef Spectra::SparseSymShiftSolve InvMatOp; - Eigen::Vector4d res; - res.setConstant(NAN); - - InvMatOp invOpt(mat); - Spectra::SymEigsShiftSolver small_eig(&invOpt, 2, 4, 0); - - small_eig.init(); - const int n_small = small_eig.compute(100000); //, 1e-8, Spectra::SMALLEST_MAGN); - if (small_eig.info() == Spectra::SUCCESSFUL) - { - res(0) = small_eig.eigenvalues()(1); - res(1) = small_eig.eigenvalues()(0); - } - - MatOp op(mat); - Spectra::SymEigsSolver large_eig(&op, 2, 4); - - large_eig.init(); - const int n_large = large_eig.compute(100000); //, 1e-8, Spectra::LARGEST_MAGN); - // std::cout< &S, StiffnessMatrix &out) - { - const int am = A.rows(); - const int sm = S.size(); - - // assert(S.minCoeff() >= 0); - // assert(S.maxCoeff() < sm); - - // Build reindexing maps for columns and rows - std::vector> SI(am); - for (int i = 0; i < sm; i++) - { - SI[S[i]].push_back(i); - } - - // Take a guess at the number of nonzeros (this assumes uniform distribution - // not banded or heavily diagonal) - std::vector> entries; - entries.reserve((A.nonZeros() / (am * am)) * (sm * sm)); - - // Iterate over outside - for (int k = 0; k < A.outerSize(); ++k) - { - // Iterate over inside - for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) - { - for (auto rit = SI[it.row()].begin(); rit != SI[it.row()].end(); rit++) - { - for (auto cit = SI[it.col()].begin(); cit != SI[it.col()].end(); cit++) - { - entries.emplace_back(*rit, *cit, it.value()); - } - } - } - } - out.resize(sm, sm); - out.setFromTriplets(entries.begin(), entries.end()); - out.makeCompressed(); - } - } // namespace -} // namespace polysolve - -Eigen::Vector4d polysolve::dirichlet_solve( - LinearSolver &solver, StiffnessMatrix &A, Eigen::VectorXd &f, - const std::vector &dirichlet_nodes, Eigen::VectorXd &u, - const int precond_num, - const std::string &save_path, - bool compute_spectrum, - const bool remove_zero_cols, - const bool skip_last_cols) -{ - // Let Γ be the set of Dirichlet dofs. - // To implement nonzero Dirichlet boundary conditions, we seek to replace - // the linear system Au = f with a new system Ãx = g, where - // - Ã is the matrix A with rows and cols of i ∈ Γ set to identity - // - g[i] = f[i] for i ∈ Γ - // - g[i] = f[i] - Σ_{j ∈ Γ} a_ij f[j] for i ∉ Γ - // In matrix terms, if we call N = diag({1 iff i ∈ Γ}), then we have that - // g = f - (I-N)*A*N*f - - int n = A.outerSize(); - Eigen::VectorXd N(n); - N.setZero(); - for (int i : dirichlet_nodes) - { - N(i) = 1; - } - - Eigen::VectorXd g = f - ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); - - // if (0) { - // Eigen::MatrixXd rhs(g.size(), 6); - // rhs.col(0) = N; - // rhs.col(1) = f; - // rhs.col(2) = N.asDiagonal() * f; - // rhs.col(3) = A * (N.asDiagonal() * f); - // rhs.col(4) = ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); - // rhs.col(5) = g; - // std::cout << rhs << std::endl; - // } - - std::vector> coeffs; - coeffs.reserve(A.nonZeros()); - assert(A.rows() == A.cols()); - for (int k = 0; k < A.outerSize(); ++k) - { - for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) - { - // it.value(); - // it.row(); // row index - // it.col(); // col index (here it is equal to k) - // it.index(); // inner index, here it is equal to it.row() - if (N(it.row()) != 1 && N(it.col()) != 1) - { - coeffs.emplace_back(it.row(), it.col(), it.value()); - } - } - } - // TODO: For numerical stability, we should set diagonal values of the same - // magnitude than the other entries in the matrix - for (int k = 0; k < n; ++k) - { - coeffs.emplace_back(k, k, N(k)); - } - // Eigen::saveMarket(A, "A_before.mat"); - A.setFromTriplets(coeffs.begin(), coeffs.end()); - A.makeCompressed(); - - // std::cout << A << std::endl; - - //remove zero cols - if (remove_zero_cols) - { - std::vector zero_col(A.cols(), true); - zero_col.back() = !skip_last_cols; - for (int k = 0; k < A.outerSize(); ++k) - { - for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) - { - if (skip_last_cols) - { - if (it.row() != A.rows() - 1 && it.col() != A.cols() - 1 && fabs(it.value()) > 1e-12) - zero_col[it.col()] = false; - } - else - { - if (fabs(it.value()) > 1e-12) - zero_col[it.col()] = false; - } - } - } - - std::vector valid; - for (int i = 0; i < A.rows(); ++i) - { - if (!zero_col[i]) - valid.push_back(i); - else if (skip_last_cols) - { - A.coeffRef(A.rows() - 1, i) = 0; - A.coeffRef(i, A.cols() - 1) = 0; - } - } - - StiffnessMatrix As; - slice(A, valid, As); - - Eigen::VectorXd gs(As.rows()); - int index = 0; - for (int i = 0; i < zero_col.size(); ++i) - { - if (!zero_col[i]) - { - gs[index++] = g[i]; - } - } - - if (u.size() != gs.size()) - { - u.resize(gs.size()); - u.setZero(); - } - - Eigen::VectorXd us = u; - solver.analyzePattern(As, precond_num); - solver.factorize(As); - solver.solve(gs, us); - f = g; - - u.resize(A.rows()); - index = 0; - for (int i = 0; i < zero_col.size(); ++i) - { - if (zero_col[i]) - { - u[i] = 0; - f[i] = 0; - } - else - u[i] = us[index++]; - } - } - else - { - // Eigen::saveMarket(A, "A.mat"); - // Eigen::saveMarketVector(g, "b.mat"); - - if (u.size() != n) - { - u.resize(n); - u.setZero(); - } - - solver.analyzePattern(A, precond_num); - solver.factorize(A); - solver.solve(g, u); - f = g; - } - - if (!save_path.empty()) - { - Eigen::saveMarket(A, save_path); - } - - if (compute_spectrum) - { - return polysolve::compute_spectrum(A); - } - else - { - return Eigen::Vector4d::Zero(); - } -} - -void polysolve::prefactorize( - LinearSolver &solver, StiffnessMatrix &A, - const std::vector &dirichlet_nodes, const int precond_num, - const std::string &save_path) -{ - int n = A.outerSize(); - Eigen::VectorXd N(n); - N.setZero(); - for (int i : dirichlet_nodes) - { - N(i) = 1; - } - - std::vector> coeffs; - coeffs.reserve(A.nonZeros()); - assert(A.rows() == A.cols()); - for (int k = 0; k < A.outerSize(); ++k) - { - for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) - { - // it.value(); - // it.row(); // row index - // it.col(); // col index (here it is equal to k) - // it.index(); // inner index, here it is equal to it.row() - if (N(it.row()) != 1 && N(it.col()) != 1) - { - coeffs.emplace_back(it.row(), it.col(), it.value()); - } - } - } - // TODO: For numerical stability, we should set diagonal values of the same - // magnitude than the other entries in the matrix - for (int k = 0; k < n; ++k) - { - coeffs.emplace_back(k, k, N(k)); - } - // Eigen::saveMarket(A, "A_before.mat"); - A.setFromTriplets(coeffs.begin(), coeffs.end()); - A.makeCompressed(); - - solver.analyzePattern(A, precond_num); - solver.factorize(A); - - if (!save_path.empty()) - { - Eigen::saveMarket(A, save_path); - } -} - -void polysolve::dirichlet_solve_prefactorized( - LinearSolver &solver, const StiffnessMatrix &A, Eigen::VectorXd &f, - const std::vector &dirichlet_nodes, Eigen::VectorXd &u) -{ - // pre-factorized version of dirichlet_solve - - int n = A.outerSize(); - Eigen::VectorXd N(n); - N.setZero(); - for (int i : dirichlet_nodes) - { - N(i) = 1; - } - - Eigen::VectorXd g = f - ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); - - if (u.size() != n) - { - u.resize(n); - u.setZero(); - } - - solver.solve(g, u); - f = g; -} diff --git a/src/polysolve/LbfgsSolver.hpp b/src/polysolve/LbfgsSolver.hpp deleted file mode 100644 index 9bc5ccc..0000000 --- a/src/polysolve/LbfgsSolver.hpp +++ /dev/null @@ -1,166 +0,0 @@ -#pragma once - -#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include - -namespace cppoptlib -{ - - template - class LbfgsSolverL2 : public ISolver - { - public: - using Superclass = ISolver; - using typename Superclass::Scalar; - using typename Superclass::THessian; - using typename Superclass::TVector; - using MatrixType = Eigen::Matrix; - - enum class LineSearch - { - Armijo, - MoreThuente, - }; - - LineSearch line_search = LineSearch::Armijo; - - void setLineSearch(const std::string &name) - { - if (name == "armijo") - { - line_search = LineSearch::Armijo; - } - else if (name == "more_thuente") - { - line_search = LineSearch::MoreThuente; - } - else - { - throw std::invalid_argument("[SparseNewtonDescentSolver] Unknown line search."); - } - polysolve::logger().debug("\tline search {}", name); - } - - void minimize(ProblemType &objFunc, TVector &x0) - { - polysolve::AssemblerUtils::instance().clear_cache(); - - const size_t m = 10; - const size_t DIM = x0.rows(); - MatrixType sVector = MatrixType::Zero(DIM, m); - MatrixType yVector = MatrixType::Zero(DIM, m); - Eigen::Matrix alpha = Eigen::Matrix::Zero(m); - TVector grad(DIM), q(DIM), grad_old(DIM), s(DIM), y(DIM); - objFunc.gradient(x0, grad); - TVector x_old = x0; - TVector x_old2 = x0; - - size_t iter = 0, globIter = 0; - Scalar H0k = 1; - this->m_current.reset(); - do - { - const Scalar relativeEpsilon = static_cast(0.0001) * std::max(static_cast(1.0), x0.norm()); - - if (grad.norm() < relativeEpsilon) - break; - - // Algorithm 7.4 (L-BFGS two-loop recursion) - q = grad; - const int k = std::min(m, iter); - - // for i = k − 1, k − 2, . . . , k − m§ - for (int i = k - 1; i >= 0; i--) - { - // alpha_i <- rho_i*s_i^T*q - const double rho = 1.0 / static_cast(sVector.col(i)).dot(static_cast(yVector.col(i))); - alpha(i) = rho * static_cast(sVector.col(i)).dot(q); - // q <- q - alpha_i*y_i - q = q - alpha(i) * yVector.col(i); - } - // r <- H_k^0*q - q = H0k * q; - // for i k − m, k − m + 1, . . . , k − 1 - for (int i = 0; i < k; i++) - { - // beta <- rho_i * y_i^T * r - const Scalar rho = 1.0 / static_cast(sVector.col(i)).dot(static_cast(yVector.col(i))); - const Scalar beta = rho * static_cast(yVector.col(i)).dot(q); - // r <- r + s_i * ( alpha_i - beta) - q = q + sVector.col(i) * (alpha(i) - beta); - } - // stop with result "H_k*f_f'=q" - - // any issues with the descent direction ? - Scalar descent = -grad.dot(q); - Scalar alpha_init = 1.0 / grad.norm(); - if (descent > -0.0001 * relativeEpsilon) - { - q = -1 * grad; - iter = 0; - alpha_init = 1.0; - } - - // find steplength - // const Scalar rate = MoreThuente::linesearch(x0, -q, objFunc, alpha_init) ; - Scalar rate; - switch (line_search) - { - case LineSearch::Armijo: - rate = Armijo::linesearch(x0, -q, objFunc); - break; - case LineSearch::MoreThuente: - rate = MoreThuente::linesearch(x0, -q, objFunc); - break; - } - // update guess - x0 = x0 - rate * q; - - grad_old = grad; - objFunc.gradient(x0, grad); - - s = x0 - x_old; - y = grad - grad_old; - - // update the history - if (iter < m) - { - sVector.col(iter) = s; - yVector.col(iter) = y; - } - else - { - - sVector.leftCols(m - 1) = sVector.rightCols(m - 1).eval(); - sVector.rightCols(1) = s; - yVector.leftCols(m - 1) = yVector.rightCols(m - 1).eval(); - yVector.rightCols(1) = y; - } - // update the scaling factor - H0k = y.dot(s) / static_cast(y.dot(y)); - - x_old = x0; - polysolve::logger().debug("\titer: {}, f = {}, ‖g‖_2 = {}", globIter, objFunc.value(x0), grad.norm()); - - iter++; - globIter++; - ++this->m_current.iterations; - this->m_current.gradNorm = grad.norm(); // template lpNorm(); - this->m_status = checkConvergence(this->m_stop, this->m_current); - } while ((objFunc.callback(this->m_current, x0)) && (this->m_status == Status::Continue)); - } - }; - -} // namespace cppoptlib diff --git a/src/polysolve/LinearSolverEigen.tpp b/src/polysolve/LinearSolverEigen.tpp deleted file mode 100644 index 574d8ba..0000000 --- a/src/polysolve/LinearSolverEigen.tpp +++ /dev/null @@ -1,144 +0,0 @@ -#pragma once - -//////////////////////////////////////////////////////////////////////////////// -#include -#include -//////////////////////////////////////////////////////////////////////////////// - -//////////////////////////////////////////////////////////////////////////////// -// Direct solvers -//////////////////////////////////////////////////////////////////////////////// - -// Get info on the last solve step -template -void polysolve::LinearSolverEigenDirect::getInfo(json ¶ms) const -{ - switch (m_Solver.info()) - { - case Eigen::Success: - params["solver_info"] = "Success"; - break; - case Eigen::NumericalIssue: - params["solver_info"] = "NumericalIssue"; - break; - case Eigen::NoConvergence: - params["solver_info"] = "NoConvergence"; - break; - case Eigen::InvalidInput: - params["solver_info"] = "InvalidInput"; - break; - default: - assert(false); - } -} - -// Analyze sparsity pattern -template -void polysolve::LinearSolverEigenDirect::analyzePattern(const StiffnessMatrix &A, const int precond_num) -{ - m_Solver.analyzePattern(A); -} - -// Factorize system matrix -template -void polysolve::LinearSolverEigenDirect::factorize(const StiffnessMatrix &A) -{ - m_Solver.factorize(A); - if (m_Solver.info() == Eigen::NumericalIssue) - { - throw std::runtime_error("[LinearSolverEigenDirect] NumericalIssue encountered."); - } -} - -// Solve the linear system -template -void polysolve::LinearSolverEigenDirect::solve( - const Ref b, Ref x) -{ - x = m_Solver.solve(b); -} - -//////////////////////////////////////////////////////////////////////////////// -// Iterative solvers -//////////////////////////////////////////////////////////////////////////////// - -// Set solver parameters -template -void polysolve::LinearSolverEigenIterative::setParameters(const json ¶ms) -{ - const std::string solver_name = name(); - if (params.contains(solver_name)) - { - if (params[solver_name].contains("max_iter")) - { - m_Solver.setMaxIterations(params[solver_name]["max_iter"]); - } - if (params[solver_name].contains("tolerance")) - { - m_Solver.setTolerance(params[solver_name]["tolerance"]); - } - } -} - -// Get info on the last solve step -template -void polysolve::LinearSolverEigenIterative::getInfo(json ¶ms) const -{ - params["solver_iter"] = m_Solver.iterations(); - params["solver_error"] = m_Solver.error(); -} - -// Analyze sparsity pattern -template -void polysolve::LinearSolverEigenIterative::analyzePattern(const StiffnessMatrix &A, const int precond_num) -{ - m_Solver.analyzePattern(A); -} - -// Factorize system matrix -template -void polysolve::LinearSolverEigenIterative::factorize(const StiffnessMatrix &A) -{ - m_Solver.factorize(A); -} - -// Solve the linear system -template -void polysolve::LinearSolverEigenIterative::solve( - const Ref b, Ref x) -{ - assert(x.size() == b.size()); - x = m_Solver.solveWithGuess(b, x); -} - -//////////////////////////////////////////////////////////////////////////////// -// Dense solvers -//////////////////////////////////////////////////////////////////////////////// - -// Get info on the last solve step -template -void polysolve::LinearSolverEigenDense::getInfo(json ¶ms) const -{ - params["solver_info"] = "Success"; -} - -template -void polysolve::LinearSolverEigenDense::factorize(const StiffnessMatrix &A) -{ - factorize_dense(Eigen::MatrixXd(A)); -} - -// Factorize system matrix -template -void polysolve::LinearSolverEigenDense::factorize_dense(const Eigen::MatrixXd &A) -{ - m_Solver.compute(A); -} - -// Solve the linear system -template -void polysolve::LinearSolverEigenDense::solve( - const Ref b, Ref x) -{ - x = m_Solver.solve(b); -} diff --git a/src/polysolve/Types.hpp b/src/polysolve/Types.hpp new file mode 100644 index 0000000..f3a44f0 --- /dev/null +++ b/src/polysolve/Types.hpp @@ -0,0 +1,19 @@ +#pragma once + +#include +#include + +#include + +namespace polysolve +{ + +#ifdef POLYSOLVE_LARGE_INDEX + typedef Eigen::SparseMatrix StiffnessMatrix; +#else + typedef Eigen::SparseMatrix StiffnessMatrix; +#endif + + using json = nlohmann::json; + +} // namespace polysolve \ No newline at end of file diff --git a/src/polysolve/Utils.cpp b/src/polysolve/Utils.cpp new file mode 100644 index 0000000..ae433dd --- /dev/null +++ b/src/polysolve/Utils.cpp @@ -0,0 +1,74 @@ +#include "Utils.hpp" + +#include + +namespace polysolve +{ + + StopWatch::StopWatch(const std::string &name, spdlog::logger &logger) + : m_name(name), m_logger(logger) + { + start(); + } + + StopWatch::StopWatch(const std::string &name, double &total_time, spdlog::logger &logger) + : m_name(name), m_total_time(&total_time), m_logger(logger) + { + start(); + } + + StopWatch::~StopWatch() + { + stop(); + } + + void StopWatch::start() + { + is_running = true; + m_start = clock::now(); + } + + void StopWatch::stop() + { + if (!is_running) + return; + m_stop = clock::now(); + + is_running = false; + log_msg(); + if (m_total_time) + *m_total_time += getElapsedTimeInSec(); + if (m_count) + ++(*m_count); + } + + double StopWatch::getElapsedTimeInSec() + { + return std::chrono::duration(m_stop - m_start).count(); + } + + void StopWatch::log_msg() + { + const static std::string log_fmt_text = + fmt::format("[{}] {{}} {{:.3g}}s", fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing")); + + if (!m_name.empty()) + { + m_logger.trace(log_fmt_text, m_name, getElapsedTimeInSec()); + } + } + + void log_and_throw_error(spdlog::logger &logger, const std::string &msg) + { + logger.error(msg); + throw std::runtime_error(msg); + } + + Eigen::SparseMatrix sparse_identity(int rows, int cols) + { + Eigen::SparseMatrix I(rows, cols); + I.setIdentity(); + return I; + } + +} // namespace polysolve \ No newline at end of file diff --git a/src/polysolve/Utils.hpp b/src/polysolve/Utils.hpp new file mode 100644 index 0000000..f9eaf10 --- /dev/null +++ b/src/polysolve/Utils.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include "Types.hpp" + +#include + +#define POLYSOLVE_SCOPED_STOPWATCH(...) polysolve::StopWatch __polysolve_stopwatch(__VA_ARGS__) + +namespace polysolve +{ + + struct Timing + { + operator double() const { return time; } + + void operator+=(const double t) + { + time += t; + ++count; + } + + double time = 0; + size_t count = 0; + }; + + class StopWatch + { + private: + using clock = std::chrono::steady_clock; + + public: + StopWatch(const std::string &name, spdlog::logger &logger); + StopWatch(const std::string &name, double &total_time, spdlog::logger &logger); + + virtual ~StopWatch(); + + void start(); + void stop(); + + double getElapsedTimeInSec(); + void log_msg(); + + private: + std::string m_name; + std::chrono::time_point m_start, m_stop; + double *m_total_time = nullptr; + size_t *m_count = nullptr; + bool is_running = false; + + spdlog::logger &m_logger; + }; + + [[noreturn]] void log_and_throw_error(spdlog::logger &logger, const std::string &msg); + + template + [[noreturn]] void log_and_throw_error(spdlog::logger &logger, const std::string &msg, const Args &...args) + { + log_and_throw_error(logger, fmt::format(msg, args...)); + } + + Eigen::SparseMatrix sparse_identity(int rows, int cols); + +} // namespace polysolve \ No newline at end of file diff --git a/src/polysolve/LinearSolverAMGCL.cpp b/src/polysolve/linear/AMGCL.cpp similarity index 88% rename from src/polysolve/LinearSolverAMGCL.cpp rename to src/polysolve/linear/AMGCL.cpp index a5cef5a..cc63864 100644 --- a/src/polysolve/LinearSolverAMGCL.cpp +++ b/src/polysolve/linear/AMGCL.cpp @@ -1,13 +1,13 @@ #ifdef POLYSOLVE_WITH_AMGCL //////////////////////////////////////////////////////////////////////////////// -#include +#include "AMGCL.hpp" #include #include //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { namespace @@ -94,7 +94,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - LinearSolverAMGCL::LinearSolverAMGCL() + AMGCL::AMGCL() { params_ = default_params(); // NOTE: usolver and psolver parameters are only used if the @@ -103,7 +103,7 @@ namespace polysolve } // Set solver parameters - void LinearSolverAMGCL::setParameters(const json ¶ms) + void AMGCL::set_parameters(const json ¶ms) { if (params.contains("AMGCL")) { @@ -114,12 +114,12 @@ namespace polysolve } if (block_size_ == 2) { - block2_solver_.setParameters(params); + block2_solver_.set_parameters(params); return; } else if (block_size_ == 3) { - block3_solver_.setParameters(params); + block3_solver_.set_parameters(params); return; } @@ -127,16 +127,16 @@ namespace polysolve } } - void LinearSolverAMGCL::getInfo(json ¶ms) const + void AMGCL::get_info(json ¶ms) const { if (block_size_ == 2) { - block2_solver_.getInfo(params); + block2_solver_.get_info(params); return; } else if (block_size_ == 3) { - block3_solver_.getInfo(params); + block3_solver_.get_info(params); return; } params["num_iterations"] = iterations_; @@ -145,7 +145,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverAMGCL::factorize(const StiffnessMatrix &Ain) + void AMGCL::factorize(const StiffnessMatrix &Ain) { if (block_size_ == 2) { @@ -187,7 +187,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverAMGCL::solve(const Eigen::Ref rhs, Eigen::Ref result) + void AMGCL::solve(const Eigen::Ref rhs, Eigen::Ref result) { if (block_size_ == 2) { @@ -213,12 +213,12 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - LinearSolverAMGCL::~LinearSolverAMGCL() + AMGCL::~AMGCL() { } template - LinearSolverAMGCL_Block::LinearSolverAMGCL_Block() + AMGCL_Block::AMGCL_Block() { params_ = default_params(); @@ -228,13 +228,13 @@ namespace polysolve } template - void LinearSolverAMGCL_Block::setParameters(const json ¶ms) + void AMGCL_Block::set_parameters(const json ¶ms) { set_params(params, params_); } template - void LinearSolverAMGCL_Block::getInfo(json ¶ms) const + void AMGCL_Block::get_info(json ¶ms) const { params["num_iterations"] = iterations_; params["final_res_norm"] = residual_error_; @@ -243,7 +243,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// template - void LinearSolverAMGCL_Block::factorize(const StiffnessMatrix &Ain) + void AMGCL_Block::factorize(const StiffnessMatrix &Ain) { assert(precond_num_ > 0); @@ -279,7 +279,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// template - void LinearSolverAMGCL_Block::solve(const Eigen::Ref rhs, Eigen::Ref result) + void AMGCL_Block::solve(const Eigen::Ref rhs, Eigen::Ref result) { assert(result.size() == rhs.size()); std::vector _rhs(rhs.data(), rhs.data() + rhs.size()); @@ -298,12 +298,12 @@ namespace polysolve } template - LinearSolverAMGCL_Block::~LinearSolverAMGCL_Block() + AMGCL_Block::~AMGCL_Block() { } - template class LinearSolverAMGCL_Block<2>; - template class LinearSolverAMGCL_Block<3>; -} // namespace polysolve + template class AMGCL_Block<2>; + template class AMGCL_Block<3>; +} // namespace polysolve::linear #endif diff --git a/src/polysolve/LinearSolverAMGCL.hpp b/src/polysolve/linear/AMGCL.hpp similarity index 76% rename from src/polysolve/LinearSolverAMGCL.hpp rename to src/polysolve/linear/AMGCL.hpp index cf78854..03a5e46 100644 --- a/src/polysolve/LinearSolverAMGCL.hpp +++ b/src/polysolve/linear/AMGCL.hpp @@ -1,9 +1,7 @@ #pragma once -#ifdef POLYSOLVE_WITH_AMGCL - //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" #include #include @@ -44,19 +42,19 @@ // column-major matrix, the solver will actually solve A^T x = b. // -namespace polysolve +namespace polysolve::linear { template - class LinearSolverAMGCL_Block : public LinearSolver + class AMGCL_Block : public Solver { public: - LinearSolverAMGCL_Block(); - ~LinearSolverAMGCL_Block(); + AMGCL_Block(); + ~AMGCL_Block(); private: - POLYSOLVE_DELETE_MOVE_COPY(LinearSolverAMGCL_Block) + POLYSOLVE_DELETE_MOVE_COPY(AMGCL_Block) public: ////////////////////// @@ -64,13 +62,13 @@ namespace polysolve ////////////////////// // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Retrieve information - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } // Factorize system matrix virtual void factorize(const StiffnessMatrix &A) override; @@ -97,15 +95,15 @@ namespace polysolve double residual_error_; }; - class LinearSolverAMGCL : public LinearSolver + class AMGCL : public Solver { public: - LinearSolverAMGCL(); - ~LinearSolverAMGCL(); + AMGCL(); + ~AMGCL(); private: - POLYSOLVE_DELETE_MOVE_COPY(LinearSolverAMGCL) + POLYSOLVE_DELETE_MOVE_COPY(AMGCL) public: ////////////////////// @@ -113,25 +111,25 @@ namespace polysolve ////////////////////// // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Retrieve information - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override { - if (block_size_==2) + if (block_size_ == 2) { - block2_solver_.analyzePattern(A, precond_num); + block2_solver_.analyze_pattern(A, precond_num); return; } else if (block_size_ == 3) { - block3_solver_.analyzePattern(A, precond_num); + block3_solver_.analyze_pattern(A, precond_num); return; } - precond_num_ = precond_num; + precond_num_ = precond_num; } // Factorize system matrix @@ -158,10 +156,8 @@ namespace polysolve size_t iterations_; double residual_error_; - LinearSolverAMGCL_Block<2> block2_solver_; - LinearSolverAMGCL_Block<3> block3_solver_; + AMGCL_Block<2> block2_solver_; + AMGCL_Block<3> block3_solver_; }; -} // namespace polysolve - -#endif +} // namespace polysolve::linear diff --git a/src/polysolve/linear/CMakeLists.txt b/src/polysolve/linear/CMakeLists.txt new file mode 100644 index 0000000..a742a35 --- /dev/null +++ b/src/polysolve/linear/CMakeLists.txt @@ -0,0 +1,21 @@ +set(SOURCES + FEMSolver.cpp + FEMSolver.hpp + Solver.cpp + Solver.hpp + AMGCL.cpp + AMGCL.hpp + CuSolverDN.cu + CuSolverDN.cuh + EigenSolver.hpp + EigenSolver.tpp + HypreSolver.cpp + HypreSolver.hpp + Pardiso.cpp + Pardiso.hpp + SaddlePointSolver.cpp + SaddlePointSolver.hpp +) + +source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES}) +target_sources(polysolve_linear PRIVATE ${SOURCES}) diff --git a/src/polysolve/LinearSolverCuSolverDN.cu b/src/polysolve/linear/CuSolverDN.cu similarity index 89% rename from src/polysolve/LinearSolverCuSolverDN.cu rename to src/polysolve/linear/CuSolverDN.cu index 4da356c..29cf757 100644 --- a/src/polysolve/LinearSolverCuSolverDN.cu +++ b/src/polysolve/linear/CuSolverDN.cu @@ -1,7 +1,7 @@ #ifdef POLYSOLVE_WITH_CUSOLVER //////////////////////////////////////////////////////////////////////////////// -#include "LinearSolverCuSolverDN.cuh" +#include "CuSolverDN.cuh" #include #include @@ -12,7 +12,7 @@ //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { namespace @@ -91,13 +91,13 @@ namespace polysolve } // namespace template - LinearSolverCuSolverDN::LinearSolverCuSolverDN() + CuSolverDN::CuSolverDN() { init(); } template - void LinearSolverCuSolverDN::init() + void CuSolverDN::init() { cusolverDnCreate(&cuHandle); cusolverDnCreateParams(&cuParams); @@ -106,18 +106,18 @@ namespace polysolve } template - void LinearSolverCuSolverDN::getInfo(json ¶ms) const + void CuSolverDN::get_info(json ¶ms) const { } template - void LinearSolverCuSolverDN::factorize(const StiffnessMatrix &A) + void CuSolverDN::factorize(const StiffnessMatrix &A) { factorize_dense(Eigen::MatrixXd(A)); } template - void LinearSolverCuSolverDN::factorize_dense(const Eigen::MatrixXd &A) + void CuSolverDN::factorize_dense(const Eigen::MatrixXd &A) { numrows = (int)A.rows(); @@ -153,7 +153,7 @@ namespace polysolve } template - void LinearSolverCuSolverDN::solve(const Ref b, Ref x) + void CuSolverDN::solve(const Ref b, Ref x) { // copy b to device if (!d_b_alloc) @@ -182,7 +182,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// template - LinearSolverCuSolverDN::~LinearSolverCuSolverDN() + CuSolverDN::~CuSolverDN() { if (d_A_alloc) { @@ -204,9 +204,9 @@ namespace polysolve cudaDeviceReset(); } -} // namespace polysolve +} // namespace polysolve::linear -template polysolve::LinearSolverCuSolverDN::LinearSolverCuSolverDN(); -template polysolve::LinearSolverCuSolverDN::LinearSolverCuSolverDN(); +template polysolve::CuSolverDN::CuSolverDN(); +template polysolve::CuSolverDN::CuSolverDN(); #endif \ No newline at end of file diff --git a/src/polysolve/LinearSolverCuSolverDN.cuh b/src/polysolve/linear/CuSolverDN.cuh similarity index 84% rename from src/polysolve/LinearSolverCuSolverDN.cuh rename to src/polysolve/linear/CuSolverDN.cuh index 586c936..8250c35 100644 --- a/src/polysolve/LinearSolverCuSolverDN.cuh +++ b/src/polysolve/linear/CuSolverDN.cuh @@ -3,7 +3,7 @@ #ifdef POLYSOLVE_WITH_CUSOLVER //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" #include #include @@ -15,18 +15,18 @@ // https://docs.nvidia.com/cuda/cusolver/index.html#cuSolverDN-function-reference // -namespace polysolve +namespace polysolve::linear { template - class LinearSolverCuSolverDN : public LinearSolver + class CuSolverDN : public Solver { public: - LinearSolverCuSolverDN(); - ~LinearSolverCuSolverDN(); + CuSolverDN(); + ~CuSolverDN(); private: - POLYSOLVE_DELETE_MOVE_COPY(LinearSolverCuSolverDN) + POLYSOLVE_DELETE_MOVE_COPY(CuSolverDN) public: ////////////////////// @@ -34,7 +34,7 @@ namespace polysolve ////////////////////// // Retrieve memory information from cuSolverDN - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Factorize system matrix (sparse) virtual void factorize(const StiffnessMatrix &A) override; @@ -42,6 +42,8 @@ namespace polysolve // Factorize system matrix (dense, preferred) virtual void factorize_dense(const Eigen::MatrixXd &A) override; + bool is_dense() const override { return true; } + // Solve the linear system Ax = b virtual void solve(const Ref b, Ref x) override; @@ -73,6 +75,6 @@ namespace polysolve int numrows; }; -} // namespace polysolve +} // namespace polysolve::linear #endif diff --git a/src/polysolve/LinearSolverEigen.hpp b/src/polysolve/linear/EigenSolver.hpp similarity index 72% rename from src/polysolve/LinearSolverEigen.hpp rename to src/polysolve/linear/EigenSolver.hpp index 0d59d30..048b706 100644 --- a/src/polysolve/LinearSolverEigen.hpp +++ b/src/polysolve/linear/EigenSolver.hpp @@ -1,16 +1,16 @@ #pragma once //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { // ----------------------------------------------------------------------------- template - class LinearSolverEigenDirect : public LinearSolver + class EigenDirect : public Solver { protected: // Solver class @@ -23,15 +23,15 @@ namespace polysolve // Name of the solver type (for debugging purposes) virtual std::string name() const override { return m_Name; } - // Constructor requires a solver name used for finding parameters in the json file passed to setParameters - LinearSolverEigenDirect(const std::string &name) { m_Name = name; } + // Constructor requires a solver name used for finding parameters in the json file passed to set_parameters + EigenDirect(const std::string &name) { m_Name = name; } public: // Get info on the last solve step - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &K, const int precond_num) override; + virtual void analyze_pattern(const StiffnessMatrix &K, const int precond_num) override; // Factorize system matrix virtual void factorize(const StiffnessMatrix &K) override; @@ -43,7 +43,7 @@ namespace polysolve // ----------------------------------------------------------------------------- template - class LinearSolverEigenIterative : public LinearSolver + class EigenIterative : public Solver { protected: // Solver class @@ -56,18 +56,18 @@ namespace polysolve // Name of the solver type (for debugging purposes) virtual std::string name() const override { return m_Name; } - // Constructor requires a solver name used for finding parameters in the json file passed to setParameters - LinearSolverEigenIterative(const std::string &name) { m_Name = name; } + // Constructor requires a solver name used for finding parameters in the json file passed to set_parameters + EigenIterative(const std::string &name) { m_Name = name; } public: // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Get info on the last solve step - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &K, const int precond_num) override; + virtual void analyze_pattern(const StiffnessMatrix &K, const int precond_num) override; // Factorize system matrix virtual void factorize(const StiffnessMatrix &K) override; @@ -81,7 +81,7 @@ namespace polysolve // ----------------------------------------------------------------------------- template - class LinearSolverEigenDense : public LinearSolver + class EigenDenseSolver : public Solver { protected: // Solver class @@ -94,12 +94,14 @@ namespace polysolve // Name of the solver type (for debugging purposes) virtual std::string name() const override { return m_Name; } - // Constructor requires a solver name used for finding parameters in the json file passed to setParameters - LinearSolverEigenDense(const std::string &name) { m_Name = name; } + // Constructor requires a solver name used for finding parameters in the json file passed to set_parameters + EigenDenseSolver(const std::string &name) { m_Name = name; } + + bool is_dense() const override { return true; } public: // Get info on the last solve step - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Factorize system matrix virtual void factorize(const StiffnessMatrix &K) override; @@ -113,8 +115,8 @@ namespace polysolve // ----------------------------------------------------------------------------- -} // namespace polysolve +} // namespace polysolve::linear //////////////////////////////////////////////////////////////////////////////// -#include +#include "EigenSolver.tpp" diff --git a/src/polysolve/linear/EigenSolver.tpp b/src/polysolve/linear/EigenSolver.tpp new file mode 100644 index 0000000..ab84874 --- /dev/null +++ b/src/polysolve/linear/EigenSolver.tpp @@ -0,0 +1,146 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "EigenSolver.hpp" +#include +//////////////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////////////// +// Direct solvers +//////////////////////////////////////////////////////////////////////////////// +namespace polysolve::linear +{ + // Get info on the last solve step + template + void EigenDirect::get_info(json ¶ms) const + { + switch (m_Solver.info()) + { + case Eigen::Success: + params["solver_info"] = "Success"; + break; + case Eigen::NumericalIssue: + params["solver_info"] = "NumericalIssue"; + break; + case Eigen::NoConvergence: + params["solver_info"] = "NoConvergence"; + break; + case Eigen::InvalidInput: + params["solver_info"] = "InvalidInput"; + break; + default: + assert(false); + } + } + + // Analyze sparsity pattern + template + void EigenDirect::analyze_pattern(const StiffnessMatrix &A, const int precond_num) + { + m_Solver.analyzePattern(A); + } + + // Factorize system matrix + template + void EigenDirect::factorize(const StiffnessMatrix &A) + { + m_Solver.factorize(A); + if (m_Solver.info() == Eigen::NumericalIssue) + { + throw std::runtime_error("[EigenDirect] NumericalIssue encountered."); + } + } + + // Solve the linear system + template + void EigenDirect::solve( + const Ref b, Ref x) + { + x = m_Solver.solve(b); + } + + //////////////////////////////////////////////////////////////////////////////// + // Iterative solvers + //////////////////////////////////////////////////////////////////////////////// + + // Set solver parameters + template + void EigenIterative::set_parameters(const json ¶ms) + { + const std::string solver_name = name(); + if (params.contains(solver_name)) + { + if (params[solver_name].contains("max_iter")) + { + m_Solver.setMaxIterations(params[solver_name]["max_iter"]); + } + if (params[solver_name].contains("tolerance")) + { + m_Solver.setTolerance(params[solver_name]["tolerance"]); + } + } + } + + // Get info on the last solve step + template + void EigenIterative::get_info(json ¶ms) const + { + params["solver_iter"] = m_Solver.iterations(); + params["solver_error"] = m_Solver.error(); + } + + // Analyze sparsity pattern + template + void EigenIterative::analyze_pattern(const StiffnessMatrix &A, const int precond_num) + { + m_Solver.analyzePattern(A); + } + + // Factorize system matrix + template + void EigenIterative::factorize(const StiffnessMatrix &A) + { + m_Solver.factorize(A); + } + + // Solve the linear system + template + void EigenIterative::solve( + const Ref b, Ref x) + { + assert(x.size() == b.size()); + x = m_Solver.solveWithGuess(b, x); + } + + //////////////////////////////////////////////////////////////////////////////// + // Dense solvers + //////////////////////////////////////////////////////////////////////////////// + + // Get info on the last solve step + template + void EigenDenseSolver::get_info(json ¶ms) const + { + params["solver_info"] = "Success"; + } + + template + void EigenDenseSolver::factorize(const StiffnessMatrix &A) + { + factorize_dense(Eigen::MatrixXd(A)); + } + + // Factorize system matrix + template + void EigenDenseSolver::factorize_dense(const Eigen::MatrixXd &A) + { + m_Solver.compute(A); + } + + // Solve the linear system + template + void EigenDenseSolver::solve( + const Ref b, Ref x) + { + x = m_Solver.solve(b); + } +} // namespace polysolve::linear diff --git a/src/polysolve/linear/FEMSolver.cpp b/src/polysolve/linear/FEMSolver.cpp new file mode 100644 index 0000000..5d1f8f8 --- /dev/null +++ b/src/polysolve/linear/FEMSolver.cpp @@ -0,0 +1,345 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "FEMSolver.hpp" + +#ifdef POLYSOLVE_WITH_SPECTRA +#include +#include +#include +#include +#endif + +#include + +//////////////////////////////////////////////////////////////////////////////// + +namespace polysolve::linear +{ + namespace + { + Eigen::Vector4d compute_spectrum(const StiffnessMatrix &mat) + { +#ifdef POLYSOLVE_WITH_SPECTRA + typedef Spectra::SparseSymMatProd MatOp; + typedef Spectra::SparseSymShiftSolve InvMatOp; + Eigen::Vector4d res; + res.setConstant(NAN); + + InvMatOp invOpt(mat); + Spectra::SymEigsShiftSolver small_eig(&invOpt, 2, 4, 0); + + small_eig.init(); + const int n_small = small_eig.compute(100000); //, 1e-8, Spectra::SMALLEST_MAGN); + if (small_eig.info() == Spectra::SUCCESSFUL) + { + res(0) = small_eig.eigenvalues()(1); + res(1) = small_eig.eigenvalues()(0); + } + + MatOp op(mat); + Spectra::SymEigsSolver large_eig(&op, 2, 4); + + large_eig.init(); + const int n_large = large_eig.compute(100000); //, 1e-8, Spectra::LARGEST_MAGN); + // std::cout< &S, StiffnessMatrix &out) + { + const int am = A.rows(); + const int sm = S.size(); + + // assert(S.minCoeff() >= 0); + // assert(S.maxCoeff() < sm); + + // Build reindexing maps for columns and rows + std::vector> SI(am); + for (int i = 0; i < sm; i++) + { + SI[S[i]].push_back(i); + } + + // Take a guess at the number of nonzeros (this assumes uniform distribution + // not banded or heavily diagonal) + std::vector> entries; + entries.reserve((A.nonZeros() / (am * am)) * (sm * sm)); + + // Iterate over outside + for (int k = 0; k < A.outerSize(); ++k) + { + // Iterate over inside + for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) + { + for (auto rit = SI[it.row()].begin(); rit != SI[it.row()].end(); rit++) + { + for (auto cit = SI[it.col()].begin(); cit != SI[it.col()].end(); cit++) + { + entries.emplace_back(*rit, *cit, it.value()); + } + } + } + } + out.resize(sm, sm); + out.setFromTriplets(entries.begin(), entries.end()); + out.makeCompressed(); + } + } // namespace + + Eigen::Vector4d dirichlet_solve( + Solver &solver, StiffnessMatrix &A, Eigen::VectorXd &f, + const std::vector &dirichlet_nodes, Eigen::VectorXd &u, + const int precond_num, + const std::string &save_path, + bool should_compute_spectrum, + const bool remove_zero_cols, + const bool skip_last_cols) + { + assert(!solver.is_dense()); + // Let Γ be the set of Dirichlet dofs. + // To implement nonzero Dirichlet boundary conditions, we seek to replace + // the linear system Au = f with a new system Ãx = g, where + // - Ã is the matrix A with rows and cols of i ∈ Γ set to identity + // - g[i] = f[i] for i ∈ Γ + // - g[i] = f[i] - Σ_{j ∈ Γ} a_ij f[j] for i ∉ Γ + // In matrix terms, if we call N = diag({1 iff i ∈ Γ}), then we have that + // g = f - (I-N)*A*N*f + + int n = A.outerSize(); + Eigen::VectorXd N(n); + N.setZero(); + for (int i : dirichlet_nodes) + { + N(i) = 1; + } + + Eigen::VectorXd g = f - ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); + + // if (0) { + // Eigen::MatrixXd rhs(g.size(), 6); + // rhs.col(0) = N; + // rhs.col(1) = f; + // rhs.col(2) = N.asDiagonal() * f; + // rhs.col(3) = A * (N.asDiagonal() * f); + // rhs.col(4) = ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); + // rhs.col(5) = g; + // std::cout << rhs << std::endl; + // } + + std::vector> coeffs; + coeffs.reserve(A.nonZeros()); + assert(A.rows() == A.cols()); + for (int k = 0; k < A.outerSize(); ++k) + { + for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) + { + // it.value(); + // it.row(); // row index + // it.col(); // col index (here it is equal to k) + // it.index(); // inner index, here it is equal to it.row() + if (N(it.row()) != 1 && N(it.col()) != 1) + { + coeffs.emplace_back(it.row(), it.col(), it.value()); + } + } + } + // TODO: For numerical stability, we should set diagonal values of the same + // magnitude than the other entries in the matrix + for (int k = 0; k < n; ++k) + { + coeffs.emplace_back(k, k, N(k)); + } + // Eigen::saveMarket(A, "A_before.mat"); + A.setFromTriplets(coeffs.begin(), coeffs.end()); + A.makeCompressed(); + + // std::cout << A << std::endl; + + // remove zero cols + if (remove_zero_cols) + { + std::vector zero_col(A.cols(), true); + zero_col.back() = !skip_last_cols; + for (int k = 0; k < A.outerSize(); ++k) + { + for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) + { + if (skip_last_cols) + { + if (it.row() != A.rows() - 1 && it.col() != A.cols() - 1 && fabs(it.value()) > 1e-12) + zero_col[it.col()] = false; + } + else + { + if (fabs(it.value()) > 1e-12) + zero_col[it.col()] = false; + } + } + } + + std::vector valid; + for (int i = 0; i < A.rows(); ++i) + { + if (!zero_col[i]) + valid.push_back(i); + else if (skip_last_cols) + { + A.coeffRef(A.rows() - 1, i) = 0; + A.coeffRef(i, A.cols() - 1) = 0; + } + } + + StiffnessMatrix As; + slice(A, valid, As); + + Eigen::VectorXd gs(As.rows()); + int index = 0; + for (int i = 0; i < zero_col.size(); ++i) + { + if (!zero_col[i]) + { + gs[index++] = g[i]; + } + } + + if (u.size() != gs.size()) + { + u.resize(gs.size()); + u.setZero(); + } + + Eigen::VectorXd us = u; + solver.analyze_pattern(As, precond_num); + solver.factorize(As); + solver.solve(gs, us); + f = g; + + u.resize(A.rows()); + index = 0; + for (int i = 0; i < zero_col.size(); ++i) + { + if (zero_col[i]) + { + u[i] = 0; + f[i] = 0; + } + else + u[i] = us[index++]; + } + } + else + { + // Eigen::saveMarket(A, "A.mat"); + // Eigen::saveMarketVector(g, "b.mat"); + + if (u.size() != n) + { + u.resize(n); + u.setZero(); + } + + solver.analyze_pattern(A, precond_num); + solver.factorize(A); + solver.solve(g, u); + f = g; + } + + if (!save_path.empty()) + { + Eigen::saveMarket(A, save_path); + } + + if (should_compute_spectrum) + { + return compute_spectrum(A); + } + else + { + return Eigen::Vector4d::Zero(); + } + } + + void prefactorize( + Solver &solver, StiffnessMatrix &A, + const std::vector &dirichlet_nodes, const int precond_num, + const std::string &save_path) + { + int n = A.outerSize(); + Eigen::VectorXd N(n); + N.setZero(); + for (int i : dirichlet_nodes) + { + N(i) = 1; + } + + std::vector> coeffs; + coeffs.reserve(A.nonZeros()); + assert(A.rows() == A.cols()); + for (int k = 0; k < A.outerSize(); ++k) + { + for (StiffnessMatrix::InnerIterator it(A, k); it; ++it) + { + // it.value(); + // it.row(); // row index + // it.col(); // col index (here it is equal to k) + // it.index(); // inner index, here it is equal to it.row() + if (N(it.row()) != 1 && N(it.col()) != 1) + { + coeffs.emplace_back(it.row(), it.col(), it.value()); + } + } + } + // TODO: For numerical stability, we should set diagonal values of the same + // magnitude than the other entries in the matrix + for (int k = 0; k < n; ++k) + { + coeffs.emplace_back(k, k, N(k)); + } + // Eigen::saveMarket(A, "A_before.mat"); + A.setFromTriplets(coeffs.begin(), coeffs.end()); + A.makeCompressed(); + + solver.analyze_pattern(A, precond_num); + solver.factorize(A); + + if (!save_path.empty()) + { + Eigen::saveMarket(A, save_path); + } + } + + void dirichlet_solve_prefactorized( + Solver &solver, const StiffnessMatrix &A, Eigen::VectorXd &f, + const std::vector &dirichlet_nodes, Eigen::VectorXd &u) + { + // pre-factorized version of dirichlet_solve + + int n = A.outerSize(); + Eigen::VectorXd N(n); + N.setZero(); + for (int i : dirichlet_nodes) + { + N(i) = 1; + } + + Eigen::VectorXd g = f - ((1.0 - N.array()).matrix()).asDiagonal() * (A * (N.asDiagonal() * f)); + + if (u.size() != n) + { + u.resize(n); + u.setZero(); + } + + solver.solve(g, u); + f = g; + } + +} // namespace polysolve::linear \ No newline at end of file diff --git a/src/polysolve/FEMSolver.hpp b/src/polysolve/linear/FEMSolver.hpp similarity index 73% rename from src/polysolve/FEMSolver.hpp rename to src/polysolve/linear/FEMSolver.hpp index 8e5d936..b4e7915 100644 --- a/src/polysolve/FEMSolver.hpp +++ b/src/polysolve/linear/FEMSolver.hpp @@ -1,17 +1,14 @@ #pragma once //////////////////////////////////////////////////////////////////////////////// -#include - -#include -using json = nlohmann::json; +#include "Solver.hpp" #include #include #include //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { /// @@ -35,18 +32,18 @@ namespace polysolve /// @param[in] dirichlet_nodes { List of ids of Dirichlet nodes } /// @param[in,out] x { Unknown vector } /// - Eigen::Vector4d dirichlet_solve(LinearSolver &solver, StiffnessMatrix &A, + Eigen::Vector4d dirichlet_solve(Solver &solver, StiffnessMatrix &A, Eigen::VectorXd &b, const std::vector &dirichlet_nodes, Eigen::VectorXd &x, const int precond_num, const std::string &save_path = "", bool compute_spectrum = false, const bool remove_zero_cols = false, const bool skip_last_cols = false); - void prefactorize(LinearSolver &solver, StiffnessMatrix &A, - const std::vector &dirichlet_nodes, const int precond_num, - const std::string &save_path = ""); + void prefactorize(Solver &solver, StiffnessMatrix &A, + const std::vector &dirichlet_nodes, const int precond_num, + const std::string &save_path = ""); - void dirichlet_solve_prefactorized(LinearSolver &solver, const StiffnessMatrix &A, Eigen::VectorXd &f, - const std::vector &dirichlet_nodes, Eigen::VectorXd &u); + void dirichlet_solve_prefactorized(Solver &solver, const StiffnessMatrix &A, Eigen::VectorXd &f, + const std::vector &dirichlet_nodes, Eigen::VectorXd &u); -} // namespace polysolve +} // namespace polysolve::linear diff --git a/src/polysolve/LinearSolverHypre.cpp b/src/polysolve/linear/HypreSolver.cpp similarity index 95% rename from src/polysolve/LinearSolverHypre.cpp rename to src/polysolve/linear/HypreSolver.cpp index 5d255b1..5a48cda 100644 --- a/src/polysolve/LinearSolverHypre.cpp +++ b/src/polysolve/linear/HypreSolver.cpp @@ -1,18 +1,18 @@ #ifdef POLYSOLVE_WITH_HYPRE //////////////////////////////////////////////////////////////////////////////// -#include +#include "HypreSolver.hpp" #include #include //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { //////////////////////////////////////////////////////////////////////////////// - LinearSolverHypre::LinearSolverHypre() + HypreSolver::HypreSolver() { precond_num_ = 0; #ifdef HYPRE_WITH_MPI @@ -35,7 +35,7 @@ namespace polysolve } // Set solver parameters - void LinearSolverHypre::setParameters(const json ¶ms) + void HypreSolver::set_parameters(const json ¶ms) { if (params.contains("Hypre")) { @@ -54,7 +54,7 @@ namespace polysolve } } - void LinearSolverHypre::getInfo(json ¶ms) const + void HypreSolver::get_info(json ¶ms) const { params["num_iterations"] = num_iterations; params["final_res_norm"] = final_res_norm; @@ -62,7 +62,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverHypre::factorize(const StiffnessMatrix &Ain) + void HypreSolver::factorize(const StiffnessMatrix &Ain) { assert(precond_num_ > 0); @@ -183,7 +183,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverHypre::solve(const Eigen::Ref rhs, Eigen::Ref result) + void HypreSolver::solve(const Eigen::Ref rhs, Eigen::Ref result) { HYPRE_IJVector b; HYPRE_ParVector par_b; @@ -295,7 +295,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - LinearSolverHypre::~LinearSolverHypre() + HypreSolver::~HypreSolver() { if (has_matrix_) { @@ -304,6 +304,6 @@ namespace polysolve } } -} // namespace polysolve +} // namespace polysolve::linear #endif diff --git a/src/polysolve/LinearSolverHypre.hpp b/src/polysolve/linear/HypreSolver.hpp similarity index 74% rename from src/polysolve/LinearSolverHypre.hpp rename to src/polysolve/linear/HypreSolver.hpp index 6458296..f0fb014 100644 --- a/src/polysolve/LinearSolverHypre.hpp +++ b/src/polysolve/linear/HypreSolver.hpp @@ -1,9 +1,7 @@ #pragma once -#ifdef POLYSOLVE_WITH_HYPRE - //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" #include #include #include @@ -19,18 +17,18 @@ // https://github.com/LLNL/hypre/blob/v2.14.0/docs/HYPRE_usr_manual.pdf // -namespace polysolve +namespace polysolve::linear { - class LinearSolverHypre : public LinearSolver + class HypreSolver : public Solver { public: - LinearSolverHypre(); - ~LinearSolverHypre(); + HypreSolver(); + ~HypreSolver(); private: - POLYSOLVE_DELETE_MOVE_COPY(LinearSolverHypre) + POLYSOLVE_DELETE_MOVE_COPY(HypreSolver) public: ////////////////////// @@ -38,13 +36,13 @@ namespace polysolve ////////////////////// // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Retrieve memory information from Pardiso - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } // Factorize system matrix virtual void factorize(const StiffnessMatrix &A) override; @@ -72,6 +70,4 @@ namespace polysolve HYPRE_ParCSRMatrix parcsr_A; }; -} // namespace polysolve - -#endif +} // namespace polysolve::linear diff --git a/src/polysolve/LinearSolverMumps.cpp b/src/polysolve/linear/Mumps.cpp similarity index 100% rename from src/polysolve/LinearSolverMumps.cpp rename to src/polysolve/linear/Mumps.cpp diff --git a/src/polysolve/LinearSolverMumps.h b/src/polysolve/linear/Mumps.hpp similarity index 100% rename from src/polysolve/LinearSolverMumps.h rename to src/polysolve/linear/Mumps.hpp diff --git a/src/polysolve/LinearSolverPardiso.cpp b/src/polysolve/linear/Pardiso.cpp similarity index 95% rename from src/polysolve/LinearSolverPardiso.cpp rename to src/polysolve/linear/Pardiso.cpp index f0b0ddf..1681a3f 100644 --- a/src/polysolve/LinearSolverPardiso.cpp +++ b/src/polysolve/linear/Pardiso.cpp @@ -1,7 +1,7 @@ #ifdef POLYSOLVE_WITH_PARDISO //////////////////////////////////////////////////////////////////////////////// -#include "LinearSolverPardiso.hpp" +#include "Pardiso.hpp" #include #ifdef POLYSOLVE_WITH_MKL #include @@ -23,26 +23,26 @@ extern "C" #endif //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { // #define PLOTS_PARDISO //////////////////////////////////////////////////////////////////////////////// - LinearSolverPardiso::LinearSolverPardiso() + Pardiso::Pardiso() : mtype(-1) { setType(11); } - void LinearSolverPardiso::setType(int _mtype) + void Pardiso::setType(int _mtype) { mtype = _mtype; init(); } // Set solver parameters - void LinearSolverPardiso::setParameters(const json ¶ms) + void Pardiso::set_parameters(const json ¶ms) { if (params.contains("Pardiso")) { @@ -53,7 +53,7 @@ namespace polysolve } } - void LinearSolverPardiso::getInfo(json ¶ms) const + void Pardiso::get_info(json ¶ms) const { params["mem_symbolic_peak"] = iparm[14]; params["mem_symbolic_perm"] = iparm[15]; @@ -64,7 +64,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverPardiso::init() + void Pardiso::init() { if (mtype == -1) { @@ -200,7 +200,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverPardiso::analyzePattern(const StiffnessMatrix &A, const int precond_num) + void Pardiso::analyze_pattern(const StiffnessMatrix &A, const int precond_num) { if (mtype == -1) { @@ -260,7 +260,7 @@ namespace polysolve // ----------------------------------------------------------------------------- - void LinearSolverPardiso::factorize(const StiffnessMatrix &A) + void Pardiso::factorize(const StiffnessMatrix &A) { if (mtype == -1) { @@ -294,7 +294,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverPardiso::solve(const Eigen::Ref rhs, Eigen::Ref result) + void Pardiso::solve(const Eigen::Ref rhs, Eigen::Ref result) { if (mtype == -1) { @@ -399,7 +399,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - void LinearSolverPardiso::freeNumericalFactorizationMemory() + void Pardiso::freeNumericalFactorizationMemory() { phase = 0; // Release internal memory #ifdef POLYSOLVE_WITH_MKL @@ -413,7 +413,7 @@ namespace polysolve //////////////////////////////////////////////////////////////////////////////// - LinearSolverPardiso::~LinearSolverPardiso() + Pardiso::~Pardiso() { if (mtype == -1) { @@ -433,6 +433,6 @@ namespace polysolve #endif } -} // namespace polysolve +} // namespace polysolve::linear #endif diff --git a/src/polysolve/LinearSolverPardiso.hpp b/src/polysolve/linear/Pardiso.hpp similarity index 85% rename from src/polysolve/LinearSolverPardiso.hpp rename to src/polysolve/linear/Pardiso.hpp index 05a75aa..00bd1ca 100644 --- a/src/polysolve/LinearSolverPardiso.hpp +++ b/src/polysolve/linear/Pardiso.hpp @@ -1,9 +1,7 @@ #pragma once -#ifdef POLYSOLVE_WITH_PARDISO - //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" #include #include #include @@ -16,18 +14,18 @@ // - OMP_NUM_THREADS: number of threads // - PARDISO_LIC_PATH: path to the folder containing the license file -namespace polysolve +namespace polysolve::linear { - class LinearSolverPardiso : public LinearSolver + class Pardiso : public Solver { public: - LinearSolverPardiso(); - ~LinearSolverPardiso(); + Pardiso(); + ~Pardiso(); private: - POLYSOLVE_DELETE_MOVE_COPY(LinearSolverPardiso) + POLYSOLVE_DELETE_MOVE_COPY(Pardiso) protected: void setType(int _mtype); @@ -40,13 +38,13 @@ namespace polysolve ////////////////////// // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Retrieve memory information from Pardiso - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override; + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override; // Factorize system matrix virtual void factorize(const StiffnessMatrix &A) override; @@ -117,6 +115,4 @@ namespace polysolve int numUniqueElements; }; -} // namespace polysolve - -#endif +} // namespace polysolve::linear diff --git a/src/polysolve/SaddlePointSolver.cpp b/src/polysolve/linear/SaddlePointSolver.cpp similarity index 88% rename from src/polysolve/SaddlePointSolver.cpp rename to src/polysolve/linear/SaddlePointSolver.cpp index ea57b5d..6d974b2 100644 --- a/src/polysolve/SaddlePointSolver.cpp +++ b/src/polysolve/linear/SaddlePointSolver.cpp @@ -1,10 +1,10 @@ -#include +#include "SaddlePointSolver.hpp" #include //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { namespace { @@ -65,7 +65,7 @@ namespace polysolve } // Set solver parameters - void SaddlePointSolver::setParameters(const json ¶ms) + void SaddlePointSolver::set_parameters(const json ¶ms) { if (params.contains("max_iter")) { @@ -102,7 +102,7 @@ namespace polysolve } } - void SaddlePointSolver::getInfo(json ¶ms) const + void SaddlePointSolver::get_info(json ¶ms) const { params["num_iterations"] = num_iterations_; params["final_res_norm"] = final_res_norm_; @@ -165,12 +165,12 @@ namespace polysolve Eigen::VectorXd alphau; Eigen::VectorXd alphap; - auto asymmetric_solver = LinearSolver::create(asymmetric_solver_name_, ""); - auto symmetric_solver = LinearSolver::create(symmetric_solver_name_, ""); - asymmetric_solver->setParameters(asymmetric_solver_params_); - symmetric_solver->setParameters(symmetric_solver_params_); + auto asymmetric_solver = Solver::create(asymmetric_solver_name_, ""); + auto symmetric_solver = Solver::create(symmetric_solver_name_, ""); + asymmetric_solver->set_parameters(asymmetric_solver_params_); + symmetric_solver->set_parameters(symmetric_solver_params_); - symmetric_solver->analyzePattern(Ss, Ss.rows()); + symmetric_solver->analyze_pattern(Ss, Ss.rows()); symmetric_solver->factorize(Ss); int i; @@ -182,33 +182,33 @@ namespace polysolve yu[i].setZero(); yp[i].setZero(); - //1 - // iters{i}.yu = gmres(As, iters{i}.Rms, iter_gmrs, eps_gm, outer_iter_gmrs); - asymmetric_solver->analyzePattern(As, As.rows()); + // 1 + // iters{i}.yu = gmres(As, iters{i}.Rms, iter_gmrs, eps_gm, outer_iter_gmrs); + asymmetric_solver->analyze_pattern(As, As.rows()); asymmetric_solver->factorize(As); assert(currentRms.size() == yu[i].size()); asymmetric_solver->solve(currentRms, yu[i]); - //2 - //Rcst = iters{i}.Rcs - Bs' * iters{i}.yu; + // 2 + // Rcst = iters{i}.Rcs - Bs' * iters{i}.yu; Rcst = currentRcs - BsT * yu[i]; - //3 - //iters{i}.yp = bicgstab(Ss, Rcst, eps_cg, 10000); + // 3 + // iters{i}.yp = bicgstab(Ss, Rcst, eps_cg, 10000); assert(Rcst.size() == yp[i].size()); symmetric_solver->solve(Rcst, yp[i]); - //4 - //Rmst = iters{i}.Rms - Bs*iters{i}.yp; + // 4 + // Rmst = iters{i}.Rms - Bs*iters{i}.yp; Rmst = currentRms - Bs * yp[i]; - //5 - // iters{i}.yu = gmres(As, Rmst, iter_gmrs, eps_gm, outer_iter_gmrs); + // 5 + // iters{i}.yu = gmres(As, Rmst, iter_gmrs, eps_gm, outer_iter_gmrs); assert(Rmst.size() == yu[i].size()); yu[i].setZero(); asymmetric_solver->solve(Rmst, yu[i]); - //update + // update Rmu.emplace_back(As * yu[i]); Rmp.emplace_back(Bs * yp[i]); Rcu.emplace_back(BsT * yu[i]); @@ -251,7 +251,7 @@ namespace polysolve // alpha = A\b; Eigen::VectorXd alpha = A.ldlt().solve(b); - // small_solver->analyzePattern(A, A.rows()); + // small_solver->analyze_pattern(A, A.rows()); // small_solver->solve(b, alpha); // alphau = alpha(1:i); @@ -259,7 +259,7 @@ namespace polysolve alphau = alpha.topRows(i + 1); alphap = alpha.bottomRows(i + 1); - //TODO stopping condition! + // TODO stopping condition! compute_solution(i + 1, alphau, alphap, yu, yp, Wm, Wc, result); final_res_norm_ = (Ain_ * result - rhs).norm(); @@ -292,4 +292,4 @@ namespace polysolve { } -} // namespace polysolve +} // namespace polysolve::linear diff --git a/src/polysolve/SaddlePointSolver.hpp b/src/polysolve/linear/SaddlePointSolver.hpp similarity index 79% rename from src/polysolve/SaddlePointSolver.hpp rename to src/polysolve/linear/SaddlePointSolver.hpp index 99c430d..0afabe1 100644 --- a/src/polysolve/SaddlePointSolver.hpp +++ b/src/polysolve/linear/SaddlePointSolver.hpp @@ -1,7 +1,7 @@ #pragma once //////////////////////////////////////////////////////////////////////////////// -#include +#include "Solver.hpp" #include #include #include @@ -9,10 +9,10 @@ //////////////////////////////////////////////////////////////////////////////// // -namespace polysolve +namespace polysolve::linear { - class SaddlePointSolver : public LinearSolver + class SaddlePointSolver : public Solver { public: @@ -28,13 +28,13 @@ namespace polysolve ////////////////////// // Set solver parameters - virtual void setParameters(const json ¶ms) override; + virtual void set_parameters(const json ¶ms) override; // Retrieve memory information from Pardiso - virtual void getInfo(json ¶ms) const override; + virtual void get_info(json ¶ms) const override; // Analyze sparsity pattern - virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } // Factorize system matrix virtual void factorize(const StiffnessMatrix &A) override; @@ -43,7 +43,7 @@ namespace polysolve virtual void solve(const Ref b, Ref x) override; // Name of the solver type (for debugging purposes) - virtual std::string name() const override { return "SaddleProblemSolver"; } + virtual std::string name() const override { return "SaddlePointSolver"; } private: StiffnessMatrix Ain_; @@ -72,4 +72,4 @@ namespace polysolve int num_iterations_; }; -} // namespace polysolve +} // namespace polysolve::linear diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/linear/Solver.cpp similarity index 70% rename from src/polysolve/LinearSolver.cpp rename to src/polysolve/linear/Solver.cpp index d9f0a47..19369c0 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/linear/Solver.cpp @@ -1,7 +1,14 @@ //////////////////////////////////////////////////////////////////////////////// -#include -#include -#include +#include + +#include "Solver.hpp" +#include "EigenSolver.hpp" +#include "SaddlePointSolver.hpp" + +#include +#include + +#include // ----------------------------------------------------------------------------- #include @@ -21,30 +28,116 @@ #include #endif #ifdef POLYSOLVE_WITH_PARDISO -#include +#include "Pardiso.hpp" #endif #ifdef POLYSOLVE_WITH_HYPRE -#include +#include "HypreSolver.hpp" #endif #ifdef POLYSOLVE_WITH_AMGCL -#include +#include "AMGCL.hpp" #endif #ifdef POLYSOLVE_WITH_CUSOLVER -#include +#include "CuSolverDN.cuh" #endif #include //////////////////////////////////////////////////////////////////////////////// -namespace polysolve +namespace polysolve::linear { + using polysolve::StiffnessMatrix; + + std::unique_ptr Solver::create(const json ¶ms_in, spdlog::logger &logger, const bool strict_validation) + { + json params = params_in; // mutable copy + + json rules; + jse::JSE jse; + + jse.strict = strict_validation; + const std::string input_spec = POLYSOLVE_LINEAR_SPEC; + std::ifstream file(input_spec); + + if (file.is_open()) + file >> rules; + else + log_and_throw_error(logger, "unable to open {} rules", input_spec); + + // set default wrt availability + for (int i = 0; i < rules.size(); i++) + { + if (rules[i]["pointer"] == "/solver") + { + rules[i]["default"] = default_solver(); + rules[i]["options"] = available_solvers(); + } + else if (rules[i]["pointer"] == "/precond") + { + rules[i]["default"] = default_precond(); + rules[i]["options"] = available_preconds(); + } + } + + // if solver is an array, pick the first available + const auto lin_solver_ptr = "/solver"_json_pointer; + if (params.contains(lin_solver_ptr) && params[lin_solver_ptr].is_array()) + { + const std::vector solvers = params[lin_solver_ptr]; + const std::vector available_solvers = Solver::available_solvers(); + std::string accepted_solver = ""; + for (const std::string &solver : solvers) + { + if (std::find(available_solvers.begin(), available_solvers.end(), solver) != available_solvers.end()) + { + accepted_solver = solver; + break; + } + } + if (!accepted_solver.empty()) + logger.info("Solver {} is the highest priority availble solver; using it.", accepted_solver); + else + logger.warn("No valid solver found in the list of specified solvers!"); + params[lin_solver_ptr] = accepted_solver; + } + + // Fallback to default linear solver if the specified solver is invalid + // NOTE: I do not know why .value() causes a segfault only on Windows + // const bool fallback_solver = params.value("/enable_overwrite_solver"_json_pointer, false); + const bool fallback_solver = + params.contains("/enable_overwrite_solver"_json_pointer) + ? params.at("/enable_overwrite_solver"_json_pointer).get() + : false; + if (fallback_solver) + { + const std::vector ss = Solver::available_solvers(); + std::string s_json = "null"; + if (!params.contains(lin_solver_ptr) || !params[lin_solver_ptr].is_string() + || std::find(ss.begin(), ss.end(), s_json = params[lin_solver_ptr].get()) == ss.end()) + { + logger.warn("Solver {} is invalid, falling back to {}", s_json, Solver::default_solver()); + params[lin_solver_ptr] = Solver::default_solver(); + } + } + + const bool valid_input = jse.verify_json(params, rules); + + if (!valid_input) + log_and_throw_error(logger, "invalid input json:\n{}", jse.log2str()); + + params = jse.inject_defaults(params, rules); + + auto res = create(params["solver"], params["precond"]); + res->set_parameters(params); + + return res; + } //////////////////////////////////////////////////////////////////////////////// #if EIGEN_VERSION_AT_LEAST(3, 3, 0) // Magic macro because C++ has no introspection -#define ENUMERATE_PRECOND(HelperFunctor, SolverType, DefaultPrecond, precond, name) \ +#define ENUMERATE_PRECOND(HelperFunctor, SolverType, default_precond, precond, name) \ do \ { \ using namespace Eigen; \ @@ -76,14 +169,14 @@ namespace polysolve else \ { \ return std::make_unique::type>(name); \ + default_precond>::type>(name); \ } \ } while (0) #else -// Magic macro because C++ has no introspection -#define ENUMERATE_PRECOND(HelperFunctor, SolverType, DefaultPrecond, precond) \ + // Magic macro because C++ has no introspection +#define ENUMERATE_PRECOND(HelperFunctor, SolverType, default_precond, precond) \ do \ { \ using namespace Eigen; \ @@ -110,7 +203,7 @@ namespace polysolve else \ { \ return std::make_unique::type>(); \ + default_precond>::type>(); \ } \ } while (0) @@ -118,18 +211,18 @@ namespace polysolve // ----------------------------------------------------------------------------- -#define RETURN_DIRECT_SOLVER_PTR(EigenSolver, Name) \ - do \ - { \ - return std::make_unique>>(Name); \ +#define RETURN_DIRECT_SOLVER_PTR(EigenSolver, Name) \ + do \ + { \ + return std::make_unique>>(Name); \ } while (0) -#define RETURN_DIRECT_DENSE_SOLVER_PTR(EigenSolver, Name) \ - do \ - { \ - return std::make_unique>>(Name); \ +#define RETURN_DIRECT_DENSE_SOLVER_PTR(EigenSolver, Name) \ + do \ + { \ + return std::make_unique>>(Name); \ } while (0) //////////////////////////////////////////////////////////////////////////////// @@ -140,14 +233,14 @@ namespace polysolve template