From bc6f295b4705f20f961e71c116eb5e43346ed51a Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 14 Jun 2022 13:10:40 -0400 Subject: [PATCH 01/29] Update Hypre Block solver --- src/polysolve/LinearSolverHypre.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/polysolve/LinearSolverHypre.cpp b/src/polysolve/LinearSolverHypre.cpp index 243bd80..1b5c9c9 100644 --- a/src/polysolve/LinearSolverHypre.cpp +++ b/src/polysolve/LinearSolverHypre.cpp @@ -32,6 +32,13 @@ namespace polysolve { if (params.contains("Hypre")) { + if (params["Hypre"].contains("block_size")) + { + if (params["Hypre"]["block_size"]==2 || params["Hypre"]["block_size"]==3) + { + dimension_ = params["Hypre"]["block_size"]; + } + } if (params["Hypre"].contains("max_iter")) { max_iter_ = params["Hypre"]["max_iter"]; From 252b87fd9f10830761a31beca526e1eaf16a3482 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 14 Jun 2022 15:37:40 -0400 Subject: [PATCH 02/29] Hypre test --- tests/test_solver.cpp | 179 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index df5478b..79f4060 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -571,3 +571,182 @@ TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); } #endif + +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("Hyprel_b2", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "gr_30_30.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(b.size()); + x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["block_size"] = 2; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); +} +#endif + +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("Hybre_crystm03", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "crystm03.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(b.size()); + x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["block_size"] = 3; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); +} +#endif + +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("hypre_smallscale", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(b.size()); + x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["block_size"] = 3; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + end = clock(); + solver->getInfo(solver_info); + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); +} +#endif From 6b301cf375600cf044230545d89ede643311e105 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Wed, 22 Jun 2022 14:10:47 -0400 Subject: [PATCH 03/29] Update the recommend parameters --- README.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 1240da1..c763030 100644 --- a/README.md +++ b/README.md @@ -45,27 +45,27 @@ Polysolve uses a json file to provide parameters to the individual solvers. The { "Eigen::LeastSquaresConjugateGradient": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Eigen::DGMRES": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Eigen::ConjugateGradient": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Eigen::BiCGSTAB": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Eigen::GMRES": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Eigen::MINRES": { "max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "Pardiso": { "mtype": -1 @@ -73,7 +73,7 @@ Polysolve uses a json file to provide parameters to the individual solvers. The "Hypre": { "max_iter": 1000, "pre_max_iter": 1000, - "tolerance": 1e-6 + "tolerance": 1e-10 }, "AMGCL": { "precond": { From 09249d764895934a6f322025ed07fef671e0d9b3 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 5 Jul 2022 19:51:10 -0400 Subject: [PATCH 04/29] Enable Hypre openmp --- cmake/recipes/hypre.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 498bddb..e46e633 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -26,7 +26,7 @@ set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) # set(HYPRE_LONG_DOUBLE ON) From d1df567339f997cdc79bf410b2985a564932380c Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 5 Jul 2022 20:24:59 -0400 Subject: [PATCH 05/29] Old json interface Old json interface and new solver option, only use for local test --- CMakeLists.txt | 2 +- cmake/recipes/json.cmake | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 243df7f..529396c 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,7 +137,7 @@ endif() # Json (MIT) include(json) -target_link_libraries(polysolve PUBLIC nlohmann_json::nlohmann_json) +target_link_libraries(polysolve PUBLIC nlohmann_json) # CHOLMOD solver if(POLYSOLVE_WITH_CHOLMOD) diff --git a/cmake/recipes/json.cmake b/cmake/recipes/json.cmake index 38c7619..73001ec 100644 --- a/cmake/recipes/json.cmake +++ b/cmake/recipes/json.cmake @@ -11,11 +11,11 @@ # # JSON MIT -if(TARGET nlohmann_json::nlohmann_json) +if(TARGET nlohmann_json) return() endif() -message(STATUS "Third-party: creating target 'nlohmann_json::nlohmann_json'") +message(STATUS "Third-party: creating target 'nlohmann_json'") # nlohmann_json is a big repo for a single header, so we just download the release archive set(NLOHMANNJSON_VERSION "v3.10.2") @@ -29,7 +29,7 @@ FetchContent_Declare( FetchContent_MakeAvailable(nlohmann_json) add_library(nlohmann_json INTERFACE) -add_library(nlohmann_json::nlohmann_json ALIAS nlohmann_json) +add_library(nlohmann_json ALIAS nlohmann_json) include(GNUInstallDirs) target_include_directories(nlohmann_json INTERFACE From 135bf3bddb56775b98b15c6081e21af079dc477b Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 5 Jul 2022 20:30:47 -0400 Subject: [PATCH 06/29] Old json interface --- CMakeLists.txt | 2 +- cmake/recipes/json.cmake | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 529396c..e522be1 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,7 +137,7 @@ endif() # Json (MIT) include(json) -target_link_libraries(polysolve PUBLIC nlohmann_json) +target_link_libraries(polysolve PUBLIC nlohmann::json) # CHOLMOD solver if(POLYSOLVE_WITH_CHOLMOD) diff --git a/cmake/recipes/json.cmake b/cmake/recipes/json.cmake index 73001ec..8aa5b09 100644 --- a/cmake/recipes/json.cmake +++ b/cmake/recipes/json.cmake @@ -11,7 +11,7 @@ # # JSON MIT -if(TARGET nlohmann_json) +if(TARGET nlohmann::json) return() endif() @@ -29,7 +29,7 @@ FetchContent_Declare( FetchContent_MakeAvailable(nlohmann_json) add_library(nlohmann_json INTERFACE) -add_library(nlohmann_json ALIAS nlohmann_json) +add_library(nlohmann::json ALIAS nlohmann_json) include(GNUInstallDirs) target_include_directories(nlohmann_json INTERFACE From 001e95b92e75667c19e29cea8888e044ede04a30 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Tue, 5 Jul 2022 21:29:12 -0400 Subject: [PATCH 07/29] Try to enable parallelism --- cmake/recipes/hypre.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index e46e633..6bf2e50 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -22,7 +22,7 @@ endif() ################################################################################ -set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) From 42ae898bcc1f27e3ec5319c615863c49e07ee9d1 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Wed, 6 Jul 2022 09:43:13 -0400 Subject: [PATCH 08/29] Turn Sequential ON(MPI ON),OPENMP off The same cmake setting as the Hypre official https://github.com/hypre-space/hypre/blob/master/src/CMakeLists.txt --- cmake/recipes/hypre.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 6bf2e50..498bddb 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -22,11 +22,11 @@ endif() ################################################################################ -set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) # set(HYPRE_LONG_DOUBLE ON) From 1afecaa79fc50ccdf060a054945ac4129f8317f4 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Wed, 6 Jul 2022 11:16:39 -0400 Subject: [PATCH 09/29] Revert "Turn Sequential OFF(MPI ON),OPENMP OFF" The same cmake setting as the Hypre official https://github.com/hypre-space/hypre/blob/master/src/CMakeLists.txt --- cmake/recipes/hypre.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 498bddb..6bf2e50 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -22,11 +22,11 @@ endif() ################################################################################ -set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) # set(HYPRE_LONG_DOUBLE ON) From 544225d9f1a5e97a33171a23a5f31a594f5a3a49 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Wed, 6 Jul 2022 11:18:08 -0400 Subject: [PATCH 10/29] Turn Sequential OFF(MPI ON),OPENMP OFF --- cmake/recipes/hypre.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 6bf2e50..53cf501 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -26,7 +26,7 @@ set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) # set(HYPRE_LONG_DOUBLE ON) From 8ab6e7dec5303fd2a3e8ccf30c7bca54c413d251 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Wed, 6 Jul 2022 15:52:37 -0400 Subject: [PATCH 11/29] Upgrade hypre to the newest, turn on Openmp trying for parallel --- cmake/recipes/hypre.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 53cf501..7a03827 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -10,7 +10,7 @@ include(FetchContent) FetchContent_Declare( hypre GIT_REPOSITORY https://github.com/hypre-space/hypre.git - GIT_TAG v2.15.1 + GIT_TAG v2.25.0 GIT_SHALLOW TRUE ) @@ -26,7 +26,7 @@ set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) # set(HYPRE_LONG_DOUBLE ON) From 68a68b5bf4cac30358c70522b8f8d2172d1c8588 Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Thu, 14 Jul 2022 18:30:32 -0400 Subject: [PATCH 12/29] Update hypre to 2.25.0, turn on OPENMP --- CMakeLists.txt | 3 ++ cmake/recipes/hypre.cmake | 66 +++++------------------------ src/polysolve/LinearSolverHypre.cpp | 45 ++++++++++++++------ src/polysolve/LinearSolverHypre.hpp | 1 - 4 files changed, 47 insertions(+), 68 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e522be1..632b64c 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -133,6 +133,9 @@ if(POLYSOLVE_WITH_HYPRE) include(hypre) target_link_libraries(polysolve PUBLIC HYPRE::HYPRE) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_HYPRE) + if(HYPRE_WITH_MPI) + target_compile_definitions(polysolve PUBLIC HYPRE_WITH_MPI) + endif() endif() # Json (MIT) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 7a03827..48bfe5a 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -6,67 +6,23 @@ endif() message(STATUS "Third-party: creating target 'HYPRE::HYPRE'") -include(FetchContent) -FetchContent_Declare( - hypre - GIT_REPOSITORY https://github.com/hypre-space/hypre.git - GIT_TAG v2.25.0 - GIT_SHALLOW TRUE -) - -FetchContent_GetProperties(hypre) -if(NOT hypre_POPULATED) - FetchContent_Populate(hypre) - file(REMOVE ${hypre_SOURCE_DIR}/src/utilities/version) -endif() - -################################################################################ - -set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) -# set(HYPRE_LONG_DOUBLE ON) -set(HYPRE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE INTERNAL "" FORCE) -add_subdirectory(${hypre_SOURCE_DIR}/src ${hypre_BINARY_DIR}) -add_library(HYPRE::HYPRE ALIAS HYPRE) - -set_property(TARGET HYPRE PROPERTY FOLDER "dependencies") - -target_include_directories(HYPRE PUBLIC ${hypre_BINARY_DIR}) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/blas) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/lapack) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/utilities) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/multivector) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/krylov) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/seq_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/parcsr_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/parcsr_block_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/distributed_matrix) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/IJ_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/matrix_matrix) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/distributed_ls) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/distributed_ls/Euclid) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/distributed_ls/ParaSails) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/parcsr_ls) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/struct_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/struct_ls) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/sstruct_mv) -target_include_directories(HYPRE PUBLIC ${hypre_SOURCE_DIR}/src/sstruct_ls) +include(FetchContent) +FetchContent_Declare( + hypre + GIT_REPOSITORY https://github.com/hypre-space/hypre.git + GIT_TAG v2.25.0 + GIT_SHALLOW TRUE +) -if(HYPRE_USING_OPENMP) - find_package(OpenMP QUIET REQUIRED) - target_link_libraries(HYPRE PUBLIC OpenMP::OpenMP_CXX) -endif() +FetchContent_MakeAvailable(hypre) -if(NOT HYPRE_SEQUENTIAL) - find_package(MPI) - if(MPI_CXX_FOUND) - target_link_libraries(HYPRE PUBLIC MPI::MPI_CXX) - endif() -endif() +add_subdirectory("${hypre_SOURCE_DIR}/src") +file(REMOVE "${hypre_SOURCE_DIR}/src/utilities/version") \ No newline at end of file diff --git a/src/polysolve/LinearSolverHypre.cpp b/src/polysolve/LinearSolverHypre.cpp index 1b5c9c9..912a7f1 100644 --- a/src/polysolve/LinearSolverHypre.cpp +++ b/src/polysolve/LinearSolverHypre.cpp @@ -4,6 +4,7 @@ #include #include +#include //////////////////////////////////////////////////////////////////////////////// namespace polysolve @@ -14,16 +15,22 @@ namespace polysolve LinearSolverHypre::LinearSolverHypre() { precond_num_ = 0; -#ifdef MPI_VERSION - /* Initialize MPI */ - int argc = 1; - char name[] = ""; - char *argv[] = {name}; - char **argvv = &argv[0]; - int myid, num_procs; - MPI_Init(&argc, &argvv); - MPI_Comm_rank(MPI_COMM_WORLD, &myid); - MPI_Comm_size(MPI_COMM_WORLD, &num_procs); +#ifdef HYPRE_WITH_MPI + int done_already; + + MPI_Initialized(&done_already); + if (!done_already) + { + /* Initialize MPI */ + int argc = 1; + char name[] = ""; + char *argv[] = {name}; + char **argvv = &argv[0]; + int myid, num_procs; + MPI_Init(&argc, &argvv); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + } #endif } @@ -75,8 +82,11 @@ namespace polysolve has_matrix_ = true; const HYPRE_Int rows = Ain.rows(); const HYPRE_Int cols = Ain.cols(); - +#ifdef HYPRE_WITH_MPI HYPRE_IJMatrixCreate(MPI_COMM_WORLD, 0, rows - 1, 0, cols - 1, &A); +#else + HYPRE_IJMatrixCreate(hypre_MPI_COMM_WORLD, 0, rows - 1, 0, cols - 1, &A); +#endif // HYPRE_IJMatrixSetPrintLevel(A, 2); HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR); HYPRE_IJMatrixInitialize(A); @@ -187,11 +197,18 @@ namespace polysolve HYPRE_IJVector x; HYPRE_ParVector par_x; +#ifdef HYPRE_WITH_MPI HYPRE_IJVectorCreate(MPI_COMM_WORLD, 0, rhs.size() - 1, &b); +#else + HYPRE_IJVectorCreate(hypre_MPI_COMM_WORLD, 0, rhs.size() - 1, &b); +#endif HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR); HYPRE_IJVectorInitialize(b); - +#ifdef HYPRE_WITH_MPI HYPRE_IJVectorCreate(MPI_COMM_WORLD, 0, rhs.size() - 1, &x); +#else + HYPRE_IJVectorCreate(hypre_MPI_COMM_WORLD, 0, rhs.size() - 1, &b); +#endif HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR); HYPRE_IJVectorInitialize(x); @@ -217,7 +234,11 @@ namespace polysolve /* Create solver */ HYPRE_Solver solver, precond; +#ifdef HYPRE_WITH_MPI HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver); +#else + HYPRE_ParCSRPCGCreate(hypre_MPI_COMM_WORLD, &solver); +#endif /* Set some parameters (See Reference Manual for more parameters) */ HYPRE_PCGSetMaxIter(solver, max_iter_); /* max iterations */ diff --git a/src/polysolve/LinearSolverHypre.hpp b/src/polysolve/LinearSolverHypre.hpp index ffc4e6f..6458296 100644 --- a/src/polysolve/LinearSolverHypre.hpp +++ b/src/polysolve/LinearSolverHypre.hpp @@ -12,7 +12,6 @@ #include #include #include -#include //////////////////////////////////////////////////////////////////////////////// // From c54adbf7c26393a098bdf88974ee0cfca000cac7 Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Tue, 19 Jul 2022 21:52:21 -0400 Subject: [PATCH 13/29] Update cuda, cannot resolve unique_ptr --- CMakeLists.txt | 40 + cmake/recipes/hypre.cmake | 2 +- src/CMakeLists.txt | 9 + src/polysolve/LinearSolver.cpp | 7 + src/polysolve/LinearSolverAMGCL_cuda.cu | 175 +++ src/polysolve/LinearSolverAMGCL_cuda.hpp | 105 ++ tests/test_solver.cpp | 1324 ++++++++++++---------- 7 files changed, 1033 insertions(+), 629 deletions(-) create mode 100644 src/polysolve/LinearSolverAMGCL_cuda.cu create mode 100644 src/polysolve/LinearSolverAMGCL_cuda.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 632b64c..76f9405 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) +option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) @@ -61,6 +62,10 @@ if(POLYSOLVE_TOPLEVEL_PROJECT) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) endif() +if(NOT DEFINED CMAKE_CUDA_STANDARD) + set(CMAKE_CUDA_STANDARD 14) + set(CMAKE_CUDA_STANDARD_REQUIRED ON) +endif() if (MSVC) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /bigobj") @@ -128,6 +133,41 @@ endif() include(eigen) target_link_libraries(polysolve PUBLIC Eigen3::Eigen) +# CUDA +if(POLYSOLVE_WITH_CUDA) + find_package(CUDAToolkit) + if (CUDAToolkit_FOUND) + set_target_properties(polysolve PROPERTIES + CUDA_SEPARABLE_COMPILATION ON CUDA_ARCHITECTURES Pascal) + include_directories(${CUDAToolkit_INCLUDE_DIRS}) + link_directories(${CUDAToolkit_LIBRARY_DIR}) + target_link_libraries(polysolve INTERFACE CUDA::cusparse CUDA::cudart_static) + endif() +endif() + +# if(POLYSOLVE_WITH_CUDA) +# find_package(CUDA) +# if (CUDA_FOUND) +# set(CUDA_TARGET_ARCH "Ampere" CACHE STRING "Target architecture(s) for CUDA") +# cuda_select_nvcc_arch_flags(CUDA_ARCH_FLAGS ${CUDA_TARGET_ARCH}) + +# if (OPENMP_FOUND) +# list(APPEND CUDA_NVCC_FLAGS -Xcompiler ${OpenMP_CXX_FLAGS}) +# endif() + +# if (CMAKE_CXX_COMPILER_ID MATCHES "GNU") +# list(APPEND CUDA_NVCC_FLAGS +# ${CUDA_ARCH_FLAGS} -Wno-deprecated-gpu-targets) + +# list(APPEND CUDA_NVCC_FLAGS -Xcompiler -fPIC -Xcompiler -Wno-vla) +# endif() + +# # add_library(cuda_target INTERFACE) +# include_directories(${CUDA_INCLUDE_DIRS}) +# link_directories(${CUDA_LIBRARY_DIR}) +# target_link_libraries(polysolve INTERFACE ${CUDA_cusparse_LIBRARY}) +# endif() +# endif() # Hypre (GNU Lesser General Public License) if(POLYSOLVE_WITH_HYPRE) include(hypre) diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 48bfe5a..904d1f7 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -6,7 +6,7 @@ endif() message(STATUS "Third-party: creating target 'HYPRE::HYPRE'") -set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e43672b..59ba6dc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,8 +14,17 @@ set(SOURCES polysolve/LinearSolverPardiso.hpp polysolve/SaddlePointSolver.cpp polysolve/SaddlePointSolver.hpp + polysolve/LinearSolverAMGCL_cuda.cu + polysolve/LinearSolverAMGCL_cuda.hpp ) +# set(CUDA_SOURCES +# polysolve/LinearSolverAMGCL_cuda.cu +# polysolve/LinearSolverAMGCL_cuda.hpp +# ) polysolve_prepend_current_path(SOURCES) polysolve_set_source_group(${SOURCES}) +# polysolve_prepend_current_path(CUDA_SOURCES) +# polysolve_set_source_group(${CUDA_SOURCES}) + target_sources(polysolve PRIVATE ${SOURCES}) diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp index 5dce34c..dd3bdaf 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/LinearSolver.cpp @@ -25,6 +25,7 @@ #endif #ifdef POLYSOLVE_WITH_AMGCL #include +#include #endif #include @@ -225,7 +226,12 @@ namespace polysolve else if (solver == "AMGCL") { return std::make_unique(); + } + else if (solver == "AMGCL_cuda") + { + return std::make_unique(); #endif + #if EIGEN_VERSION_AT_LEAST(3, 3, 0) // Available only with Eigen 3.3.0 and newer #ifndef POLYSOLVE_LARGE_INDEX @@ -294,6 +300,7 @@ namespace polysolve #endif #ifdef POLYSOLVE_WITH_AMGCL "AMGCL", + "AMGCL_cuda", #endif #if EIGEN_VERSION_AT_LEAST(3, 3, 0) #ifndef POLYSOLVE_LARGE_INDEX diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu new file mode 100644 index 0000000..e768eaf --- /dev/null +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -0,0 +1,175 @@ +#ifdef POLYSOLVE_WITH_AMGCL + +//////////////////////////////////////////////////////////////////////////////// +#include + +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace polysolve +{ + + namespace + { + /* https://stackoverflow.com/questions/15904896/range-based-for-loop-on-a-dynamic-array */ + template + struct WrappedArray + { + WrappedArray(const T *first, const T *last) + : begin_{first}, end_{last} {} + WrappedArray(const T *first, std::ptrdiff_t size) + : WrappedArray{first, first + size} {} + + const T *begin() const noexcept { return begin_; } + const T *end() const noexcept { return end_; } + const T &operator[](const size_t i) const { return begin_[i]; } + + const T *begin_; + const T *end_; + }; + + json default_params() + { + json params = R"({ + "precond": { + "relax": { + "degree": 16, + "type": "chebyshev", + "power_iters": 100, + "higher": 2, + "lower": 0.008333333333, + "scale": true + }, + "class": "amg", + "max_levels": 6, + "direct_coarse": false, + "ncycle": 2, + "coarsening": { + "type": "smoothed_aggregation", + "estimate_spectral_radius": true, + "relax": 1, + "aggr": { + "eps_strong": 0 + } + } + }, + "solver": { + "tol": 1e-10, + "maxiter": 1000, + "type": "cg" + } + })"_json; + + return params; + } + + void set_params(const json ¶ms, json &out) + { + if (params.contains("AMGCL")) + { + // Patch the stored params with input ones + if (params["AMGCL"].contains("precond")) + out["precond"].merge_patch(params["AMGCL"]["precond"]); + if (params["AMGCL"].contains("solver")) + out["solver"].merge_patch(params["AMGCL"]["solver"]); + + if (out["precond"]["class"] == "schur_pressure_correction") + { + // Initialize the u and p solvers with a tolerance that is comparable to the main solver's + if (!out["precond"].contains("usolver")) + { + out["precond"]["usolver"] = R"({"solver": {"maxiter": 100}})"_json; + out["precond"]["usolver"]["solver"]["tol"] = 10 * out["solver"]["tol"].get(); + } + if (!out["precond"].contains("usolver")) + { + out["precond"]["psolver"] = R"({"solver": {"maxiter": 100}})"_json; + out["precond"]["psolver"]["solver"]["tol"] = 10 * out["solver"]["tol"].get(); + } + } + } + } + } // namespace + + //////////////////////////////////////////////////////////////////////////////// + + LinearSolverAMGCL_cuda::LinearSolverAMGCL_cuda() + { + params_ = default_params(); + // NOTE: usolver and psolver parameters are only used if the + // preconditioner class is "schur_pressure_correction" + precond_num_ = 0; + cusparseCreate(&backend_params_.cusparse_handle); + } + + // Set solver parameters + void LinearSolverAMGCL_cuda::setParameters(const json ¶ms) + { + if (params.contains("AMGCL")) + { + set_params(params, params_); + } + } + + void LinearSolverAMGCL_cuda::getInfo(json ¶ms) const + { + params["num_iterations"] = iterations_; + params["final_res_norm"] = residual_error_; + } + + //////////////////////////////////////////////////////////////////////////////// + + void LinearSolverAMGCL_cuda::factorize(const StiffnessMatrix &Ain) + { + assert(precond_num_ > 0); + + int numRows = Ain.rows(); + + WrappedArray ia(Ain.outerIndexPtr(), numRows + 1); + WrappedArray ja(Ain.innerIndexPtr(), Ain.nonZeros()); + WrappedArray a(Ain.valuePtr(), Ain.nonZeros()); + if (params_["precond"]["class"] == "schur_pressure_correction") + { + std::vector pmask(numRows, 0); + for (size_t i = precond_num_; i < numRows; ++i) + pmask[i] = 1; + params_["precond"]["pmask"] = pmask; + } + + // AMGCL takes the parameters as a Boost property_tree (i.e., another JSON data structure) + std::stringstream ss_params; + ss_params << params_; + boost::property_tree::ptree pt_params; + boost::property_tree::read_json(ss_params, pt_params); + auto A = std::tie(numRows, ia, ja, a); + solver_ = std::make_unique(A, pt_params,backend_params_); + iterations_ = 0; + residual_error_ = 0; + } + + //////////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////////////////////////////////////////// + + void LinearSolverAMGCL_cuda::solve(const Eigen::Ref rhs, Eigen::Ref result) + { + assert(result.size() == rhs.size()); + std::vector _rhs(rhs.data(), rhs.data() + rhs.size()); + std::vector x(result.data(), result.data() + result.size()); + auto rhs_b = Backend::copy_vector(_rhs, backend_params_); + auto x_b = Backend::copy_vector(x, backend_params_); + + std::tie(iterations_, residual_error_) = (*solver_)(*rhs_b, *x_b); + thrust::copy(x_b->begin(), x_b->end(),result.data()); + } + + //////////////////////////////////////////////////////////////////////////////// + + LinearSolverAMGCL_cuda::~LinearSolverAMGCL_cuda() + { + } + +} // namespace polysolve + +#endif diff --git a/src/polysolve/LinearSolverAMGCL_cuda.hpp b/src/polysolve/LinearSolverAMGCL_cuda.hpp new file mode 100644 index 0000000..5d9593b --- /dev/null +++ b/src/polysolve/LinearSolverAMGCL_cuda.hpp @@ -0,0 +1,105 @@ +#pragma once + +#ifdef POLYSOLVE_WITH_AMGCL + +//////////////////////////////////////////////////////////////////////////////// +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// #include +#include +#include +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////// +// +// WARNING: +// The matrix is assumed to be in row-major format, since AMGCL assumes that the +// outer index is for row. If the matrix is symmetric, you are fine, because CSR +// and CSC are the same. If the matrix is not symmetric and you pass in a +// column-major matrix, the solver will actually solve A^T x = b. +// + +namespace polysolve +{ + class LinearSolverAMGCL_cuda : public LinearSolver + { + + public: + LinearSolverAMGCL_cuda(); + ~LinearSolverAMGCL_cuda(); + + private: + POLYSOLVE_DELETE_MOVE_COPY(LinearSolverAMGCL_cuda) + + public: + ////////////////////// + // Public interface // + ////////////////////// + + // Set solver parameters + virtual void setParameters(const json ¶ms) override; + + // Retrieve information + virtual void getInfo(json ¶ms) const override; + + // Analyze sparsity pattern + virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override + { + precond_num_ = precond_num; + } + + // Factorize system matrix + virtual void factorize(const StiffnessMatrix &A) override; + + // Solve the linear system Ax = b + virtual void solve(const Ref b, Ref x) override; + + // Name of the solver type (for debugging purposes) + virtual std::string name() const override { return "AMGCL_cuda"; } + + private: + using Backend = amgcl::backend::cuda; + using Solver = amgcl::make_solver< + amgcl::runtime::preconditioner, + amgcl::runtime::solver::wrapper>; + std::unique_ptr solver_; + json params_; + typename Backend::params backend_params_; + + int precond_num_; + int block_size_ = 1; + + // Output info + size_t iterations_; + double residual_error_; + }; + +} // namespace polysolve + +#endif diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index 79f4060..fcbb93c 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include ////////////////////////////////////////////////////////////////////////// using namespace polysolve; @@ -37,549 +37,699 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) fin.close(); A.setFromTriplets(triple.begin(), triple.end()); }; -TEST_CASE("all", "[solver]") +void loadSymmetric_test(Eigen::SparseMatrix &A, std::string PATH) { - const std::string path = POLYSOLVE_DATA_DIR; - Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - auto solvers = LinearSolver::availableSolvers(); - - for (const auto &s : solvers) + std::ifstream fin(PATH); + long int M, N, L; + while (fin.peek() == '%') { - if (s == "Eigen::DGMRES") - continue; -#ifdef WIN32 - if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") - continue; -#endif - auto solver = LinearSolver::create(s, ""); - json params; - params[s]["tolerance"] = 1e-10; - solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(b.size()); - x.setZero(); - - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - - // solver->getInfo(solver_info); - - // std::cout<<"Solver error: "< A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - auto solvers = LinearSolver::availableSolvers(); - - for (const auto &s : solvers) + fin >> M >> N >> L; + A.resize(M, N); + A.reserve(L * 2 - M); + std::vector> triple; + for (size_t i = 0; i < L; i++) { - if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES" || s == "Eigen::LeastSquaresConjugateGradient" || s == "Eigen::DGMRES") + int m, n; + double data; + fin >> m >> n >> data; + if (m != n) { - auto solver = LinearSolver::create(s, ""); - json params; - params[s]["max_iter"] = 1000; - params[s]["tolerance"] = 1e-10; - solver->setParameters(params); - - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(b.size()); - x.setZero(); - - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - - // solver->getInfo(solver_info); - - // std::cout<<"Solver error: "<(m - 1, n - 1, data/2)); + triple.push_back(Eigen::Triplet(n - 1, m - 1, data/2)); } - } -} - -TEST_CASE("pre_factor", "[solver]") -{ - const std::string path = POLYSOLVE_DATA_DIR; - Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - auto solvers = LinearSolver::availableSolvers(); - - for (const auto &s : solvers) - { - if (s == "Eigen::DGMRES") - continue; -#ifdef WIN32 - if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") - continue; -#endif - auto solver = LinearSolver::create(s, ""); - solver->analyzePattern(A, A.rows()); - - std::default_random_engine eng{42}; - std::uniform_real_distribution urd(0.1, 5); - - for (int i = 0; i < 10; ++i) - { - std::vector> tripletList; - - for (int k = 0; k < A.outerSize(); ++k) - { - for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) - { - if (it.row() == it.col()) - { - tripletList.emplace_back(it.row(), it.col(), urd(eng) * 100); - } - else if (it.row() < it.col()) - { - const double val = -urd(eng); - tripletList.emplace_back(it.row(), it.col(), val); - tripletList.emplace_back(it.col(), it.row(), val); - } - } - } - - Eigen::SparseMatrix Atmp(A.rows(), A.cols()); - Atmp.setFromTriplets(tripletList.begin(), tripletList.end()); - - Eigen::VectorXd b(Atmp.rows()); - b.setRandom(); - Eigen::VectorXd x(b.size()); - x.setZero(); - - solver->factorize(Atmp); - solver->solve(b, x); - - // solver->getInfo(solver_info); - - // std::cout<<"Solver error: "<(m - 1, n - 1, data)); } } -} - -#ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("hypre", "[solver]") -{ - const std::string path = POLYSOLVE_DATA_DIR; - Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - auto solver = LinearSolver::create("Hypre", ""); - // solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(b.size()); - x.setZero(); - - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - - // solver->getInfo(solver_info); - - // std::cout<<"Solver error: "< A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - // solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(A.rows()); - x.setZero(); - { - json solver_info; - - auto solver = LinearSolver::create("Hypre", ""); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] > 1); - } - - { - json solver_info; - auto solver = LinearSolver::create("Hypre", ""); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] == 1); - } - - // std::cout<<"Solver error: "< A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - // solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(A.rows()); - x.setZero(); - { - json solver_info; - - auto solver = LinearSolver::create("AMGCL", ""); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] > 0); - } - - { - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] == 0); - } - - // std::cout<<"Solver error: "< A; - bool ok = loadMarket(A, path + "/A0.mat"); - REQUIRE(ok); - - Eigen::VectorXd b; - ok = loadMarketVector(b, path + "/b0.mat"); - REQUIRE(ok); - - auto solver = LinearSolver::create("SaddlePointSolver", ""); - solver->analyzePattern(A, 9934); - solver->factorize(A); - Eigen::VectorXd x(A.rows()); - solver->solve(b, x); - const double err = (A * x - b).norm(); - REQUIRE(err < 1e-8); -} - -#ifdef POLYSOLVE_WITH_AMGCL -TEST_CASE("amgcl_blocksolver_small_scale", "[solver]") -{ -#ifndef NDEBUG - return; -#endif - const std::string path = POLYSOLVE_DATA_DIR; - Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - // solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(A.rows()); - Eigen::VectorXd x_b(A.rows()); - x.setZero(); - x_b.setZero(); - { - json solver_info; - - auto solver = LinearSolver::create("AMGCL", ""); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] > 0); - const double err = (A * x - b).norm(); - REQUIRE(err < 1e-5); - } - - { - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["block_size"] = 3; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x_b); - solver->getInfo(solver_info); - - REQUIRE(solver_info["num_iterations"] > 0); - const double err = (A * x_b - b).norm(); - REQUIRE(err < 1e-5); - } -} -#endif - -#ifdef POLYSOLVE_WITH_AMGCL -TEST_CASE("amgcl_blocksolver_b2", "[solver]") -{ -#ifndef NDEBUG - return; -#endif - const std::string path = POLYSOLVE_DATA_DIR; - std::string MatrixName = "gr_30_30.mtx"; - Eigen::SparseMatrix A; - loadSymmetric(A, path + "/" + MatrixName); - - std::cout << "Matrix Load OK" << std::endl; - - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(A.rows()); - Eigen::VectorXd x_b(A.rows()); - x.setOnes(); - x_b.setOnes(); - { - amgcl::profiler<> prof("gr_30_30_Scalar"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 1000; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - { - amgcl::profiler<> prof("gr_30_30_Block"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 1000; - params["AMGCL"]["block_size"] = 2; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x_b); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); - REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -} -#endif - -#ifdef POLYSOLVE_WITH_AMGCL -TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]") -{ -#ifndef NDEBUG - return; -#endif - std::cout << "Polysolve AMGCL Solver" << std::endl; - const std::string path = POLYSOLVE_DATA_DIR; - std::string MatrixName = "crystm03.mtx"; - Eigen::SparseMatrix A; - loadSymmetric(A, path + "/" + MatrixName); - std::cout << "Matrix Load OK" << std::endl; - Eigen::VectorXd b(A.rows()); - b.setOnes(); - Eigen::VectorXd x_b(A.rows()); - x_b.setZero(); - Eigen::VectorXd x(A.rows()); - x.setZero(); - { - amgcl::profiler<> prof("crystm03_Block"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 1000; - params["AMGCL"]["block_size"] = 3; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x_b); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - { - amgcl::profiler<> prof("crystm03_Scalar"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 10000; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); - REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -} -#endif - -#ifdef POLYSOLVE_WITH_AMGCL -TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") -{ -#ifndef NDEBUG - return; -#endif - std::cout << "Polysolve AMGCL Solver" << std::endl; - const std::string path = POLYSOLVE_DATA_DIR; - std::string MatrixName = "crystm03.mtx"; - Eigen::SparseMatrix A; - loadSymmetric(A, path + "/" + MatrixName); - - std::cout << "Matrix Load OK" << std::endl; - - Eigen::VectorXd b(A.rows()); - b.setOnes(); - Eigen::VectorXd x_b(A.rows()); - x_b.setZero(); - Eigen::VectorXd x(A.rows()); - x.setZero(); - { - amgcl::profiler<> prof("crystm03_Block"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 10000; - params["AMGCL"]["block_size"] = 3; - params["AMGCL"]["solver_type"] = "bicgstab"; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x_b); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - { - amgcl::profiler<> prof("crystm03_Scalar"); - json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); - json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 10000; - params["AMGCL"]["solver_type"] = "bicgstab"; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x); - prof.toc("solve"); - solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl - << prof << std::endl; - } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); - REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -} -#endif + fin.close(); + A.setFromTriplets(triple.begin(), triple.end()); +}; +// TEST_CASE("all", "[solver]") +// { +// const std::string path = POLYSOLVE_DATA_DIR; +// Eigen::SparseMatrix A; +// const bool ok = loadMarket(A, path + "/A_2.mat"); +// REQUIRE(ok); + +// auto solvers = LinearSolver::availableSolvers(); + +// for (const auto &s : solvers) +// { +// if (s == "Eigen::DGMRES") +// continue; +// #ifdef WIN32 +// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") +// continue; +// #endif +// auto solver = LinearSolver::create(s, ""); +// json params; +// params[s]["tolerance"] = 1e-10; +// solver->setParameters(params); +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); + +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); + +// // solver->getInfo(solver_info); + +// // std::cout<<"Solver error: "< A; +// const bool ok = loadMarket(A, path + "/A_2.mat"); +// REQUIRE(ok); + +// auto solvers = LinearSolver::availableSolvers(); + +// for (const auto &s : solvers) +// { +// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES" || s == "Eigen::LeastSquaresConjugateGradient" || s == "Eigen::DGMRES") +// { +// auto solver = LinearSolver::create(s, ""); +// json params; +// params[s]["max_iter"] = 1000; +// params[s]["tolerance"] = 1e-10; +// solver->setParameters(params); + +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); + +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); + +// // solver->getInfo(solver_info); + +// // std::cout<<"Solver error: "< A; +// const bool ok = loadMarket(A, path + "/A_2.mat"); +// REQUIRE(ok); + +// auto solvers = LinearSolver::availableSolvers(); + +// for (const auto &s : solvers) +// { +// if (s == "Eigen::DGMRES") +// continue; +// #ifdef WIN32 +// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") +// continue; +// #endif +// auto solver = LinearSolver::create(s, ""); +// solver->analyzePattern(A, A.rows()); + +// std::default_random_engine eng{42}; +// std::uniform_real_distribution urd(0.1, 5); + +// for (int i = 0; i < 10; ++i) +// { +// std::vector> tripletList; + +// for (int k = 0; k < A.outerSize(); ++k) +// { +// for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) +// { +// if (it.row() == it.col()) +// { +// tripletList.emplace_back(it.row(), it.col(), urd(eng) * 100); +// } +// else if (it.row() < it.col()) +// { +// const double val = -urd(eng); +// tripletList.emplace_back(it.row(), it.col(), val); +// tripletList.emplace_back(it.col(), it.row(), val); +// } +// } +// } + +// Eigen::SparseMatrix Atmp(A.rows(), A.cols()); +// Atmp.setFromTriplets(tripletList.begin(), tripletList.end()); + +// Eigen::VectorXd b(Atmp.rows()); +// b.setRandom(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); + +// solver->factorize(Atmp); +// solver->solve(b, x); + +// // solver->getInfo(solver_info); + +// // std::cout<<"Solver error: "< A; +// // const bool ok = loadMarket(A, path + "/A_2.mat"); +// // REQUIRE(ok); +// loadSymmetric_test(A,path + "/A_2.mat"); +// auto solver = LinearSolver::create("Hypre", ""); +// // solver->setParameters(params); +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); + +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); + +// // solver->getInfo(solver_info); + +// // std::cout<<"Solver error: "< A; +// // const bool ok = loadMarket(A, path + "/A_2.mat"); +// // REQUIRE(ok); + +// loadSymmetric_test(A,path + "/A_2.mat"); +// // solver->setParameters(params); +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(A.rows()); +// x.setZero(); +// { +// json solver_info; + +// auto solver = LinearSolver::create("Hypre", ""); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] > 1); +// } + +// { +// json solver_info; +// auto solver = LinearSolver::create("Hypre", ""); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); + +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] == 1); +// } + +// // std::cout<<"Solver error: "< A; +// const bool ok = loadMarket(A, path + "/A_2.mat"); +// REQUIRE(ok); + +// // solver->setParameters(params); +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(A.rows()); +// x.setZero(); +// { +// json solver_info; + +// auto solver = LinearSolver::create("AMGCL", ""); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] > 0); +// } + +// { +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); + +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] == 0); +// } + +// // std::cout<<"Solver error: "< A; +// bool ok = loadMarket(A, path + "/A0.mat"); +// REQUIRE(ok); + +// Eigen::VectorXd b; +// ok = loadMarketVector(b, path + "/b0.mat"); +// REQUIRE(ok); + +// auto solver = LinearSolver::create("SaddlePointSolver", ""); +// solver->analyzePattern(A, 9934); +// solver->factorize(A); +// Eigen::VectorXd x(A.rows()); +// solver->solve(b, x); +// const double err = (A * x - b).norm(); +// REQUIRE(err < 1e-8); +// } + +// #ifdef POLYSOLVE_WITH_AMGCL +// TEST_CASE("amgcl_blocksolver_small_scale", "[solver]") +// { +// #ifndef NDEBUG +// return; +// #endif +// const std::string path = POLYSOLVE_DATA_DIR; +// Eigen::SparseMatrix A; +// const bool ok = loadMarket(A, path + "/A_2.mat"); +// REQUIRE(ok); + +// // solver->setParameters(params); +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(A.rows()); +// Eigen::VectorXd x_b(A.rows()); +// x.setZero(); +// x_b.setZero(); +// { +// json solver_info; + +// auto solver = LinearSolver::create("AMGCL", ""); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] > 0); +// const double err = (A * x - b).norm(); +// REQUIRE(err < 1e-5); +// } + +// { +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["block_size"] = 3; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x_b); +// solver->getInfo(solver_info); + +// REQUIRE(solver_info["num_iterations"] > 0); +// const double err = (A * x_b - b).norm(); +// REQUIRE(err < 1e-5); +// } +// } +// #endif + +// #ifdef POLYSOLVE_WITH_AMGCL +// TEST_CASE("amgcl_blocksolver_b2", "[solver]") +// { +// #ifndef NDEBUG +// return; +// #endif +// const std::string path = POLYSOLVE_DATA_DIR; +// std::string MatrixName = "gr_30_30.mtx"; +// Eigen::SparseMatrix A; +// loadSymmetric(A, path + "/" + MatrixName); + +// std::cout << "Matrix Load OK" << std::endl; + +// Eigen::VectorXd b(A.rows()); +// b.setRandom(); +// Eigen::VectorXd x(A.rows()); +// Eigen::VectorXd x_b(A.rows()); +// x.setOnes(); +// x_b.setOnes(); +// { +// amgcl::profiler<> prof("gr_30_30_Scalar"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 1000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// { +// amgcl::profiler<> prof("gr_30_30_Block"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 1000; +// params["AMGCL"]["block_size"] = 2; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x_b); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); +// } +// #endif + +// #ifdef POLYSOLVE_WITH_AMGCL +// TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]") +// { +// #ifndef NDEBUG +// return; +// #endif +// std::cout << "Polysolve AMGCL Solver" << std::endl; +// const std::string path = POLYSOLVE_DATA_DIR; +// std::string MatrixName = "crystm03.mtx"; +// Eigen::SparseMatrix A; +// loadSymmetric(A, path + "/" + MatrixName); +// std::cout << "Matrix Load OK" << std::endl; +// Eigen::VectorXd b(A.rows()); +// b.setOnes(); +// Eigen::VectorXd x_b(A.rows()); +// x_b.setZero(); +// Eigen::VectorXd x(A.rows()); +// x.setZero(); +// { +// amgcl::profiler<> prof("crystm03_Block"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 1000; +// params["AMGCL"]["block_size"] = 3; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x_b); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// { +// amgcl::profiler<> prof("crystm03_Scalar"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 10000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); +// } +// #endif + +// #ifdef POLYSOLVE_WITH_AMGCL +// TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") +// { +// #ifndef NDEBUG +// return; +// #endif +// std::cout << "Polysolve AMGCL Solver" << std::endl; +// const std::string path = POLYSOLVE_DATA_DIR; +// std::string MatrixName = "crystm03.mtx"; +// Eigen::SparseMatrix A; +// loadSymmetric(A, path + "/" + MatrixName); + +// std::cout << "Matrix Load OK" << std::endl; + +// Eigen::VectorXd b(A.rows()); +// b.setOnes(); +// Eigen::VectorXd x_b(A.rows()); +// x_b.setZero(); +// Eigen::VectorXd x(A.rows()); +// x.setZero(); +// { +// amgcl::profiler<> prof("crystm03_Block"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 10000; +// params["AMGCL"]["block_size"] = 3; +// params["AMGCL"]["solver_type"] = "bicgstab"; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x_b); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// { +// amgcl::profiler<> prof("crystm03_Scalar"); +// json solver_info; +// auto solver = LinearSolver::create("AMGCL", ""); +// prof.tic("setup"); +// json params; +// params["AMGCL"]["tolerance"] = 1e-8; +// params["AMGCL"]["max_iter"] = 10000; +// params["AMGCL"]["solver_type"] = "bicgstab"; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// prof.toc("setup"); +// prof.tic("solve"); +// solver->solve(b, x); +// prof.toc("solve"); +// solver->getInfo(solver_info); +// REQUIRE(solver_info["num_iterations"] > 0); +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl +// << prof << std::endl; +// } +// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); +// } +// #endif + +// #ifdef POLYSOLVE_WITH_HYPRE +// TEST_CASE("Hypre_b2", "[solver]") +// { +// const std::string path = POLYSOLVE_DATA_DIR; +// std::string MatrixName = "gr_30_30.mtx"; +// Eigen::SparseMatrix A; +// loadSymmetric(A, path + "/" + MatrixName); +// std::cout << "Matrix Load OK" << std::endl; +// Eigen::VectorXd b(A.rows()); +// b.setOnes(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); +// Eigen::VectorXd x_b(b.size()); +// x_b.setZero(); +// { +// clock_t start, end; +// json solver_info; +// start = clock(); +// auto solver = LinearSolver::create("Hypre", ""); +// json params; +// params["Hypre"]["tolerance"] = 1e-8; +// params["Hypre"]["max_iter"] = 1000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); +// end = clock(); +// solver->getInfo(solver_info); +// std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl; +// } +// { +// clock_t start, end; +// json solver_info; +// start = clock(); +// auto solver = LinearSolver::create("Hypre", ""); +// json params; +// params["Hypre"]["block_size"] = 2; +// params["Hypre"]["tolerance"] = 1e-8; +// params["Hypre"]["max_iter"] = 1000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x_b); +// end = clock(); +// solver->getInfo(solver_info); +// std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl; +// } +// const double err = (A * x - b).norm() / b.norm(); +// const double err_b = (A * x_b - b).norm() / b.norm(); +// std::cout << "Scalar relative error " << err << std::endl; +// std::cout << "Block relative error " << err_b << std::endl; +// REQUIRE(err < 1e-8); +// REQUIRE(err_b < 1e-8); +// } +// #endif + +// #ifdef POLYSOLVE_WITH_HYPRE +// TEST_CASE("Hypre_crystm03", "[solver]") +// { +// const std::string path = POLYSOLVE_DATA_DIR; +// std::string MatrixName = "crystm03.mtx"; +// Eigen::SparseMatrix A; +// loadSymmetric(A, path + "/" + MatrixName); +// std::cout << "Matrix Load OK" << std::endl; +// Eigen::VectorXd b(A.rows()); +// b.setOnes(); +// Eigen::VectorXd x(b.size()); +// x.setZero(); +// Eigen::VectorXd x_b(b.size()); +// x_b.setZero(); +// { +// clock_t start, end; +// json solver_info; +// start = clock(); +// auto solver = LinearSolver::create("Hypre", ""); +// json params; +// params["Hypre"]["tolerance"] = 1e-8; +// params["Hypre"]["max_iter"] = 1000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x); +// end = clock(); +// solver->getInfo(solver_info); +// std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl; +// } +// { +// clock_t start, end; +// json solver_info; +// start = clock(); +// auto solver = LinearSolver::create("Hypre", ""); +// json params; +// params["Hypre"]["block_size"] = 3; +// params["Hypre"]["tolerance"] = 1e-8; +// params["Hypre"]["max_iter"] = 1000; +// solver->setParameters(params); +// solver->analyzePattern(A, A.rows()); +// solver->factorize(A); +// solver->solve(b, x_b); +// end = clock(); +// solver->getInfo(solver_info); +// std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; +// std::cout << solver_info["num_iterations"] << std::endl; +// std::cout << solver_info["final_res_norm"] << std::endl; +// } +// const double err = (A * x - b).norm() / b.norm(); +// const double err_b = (A * x_b - b).norm() / b.norm(); +// std::cout << "Scalar relative error " << err << std::endl; +// std::cout << "Block relative error " << err_b << std::endl; +// REQUIRE(err < 1e-8); +// REQUIRE(err_b < 1e-8); +// } +// #endif #ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("Hyprel_b2", "[solver]") +TEST_CASE("hypre_smallscale", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; - std::string MatrixName = "gr_30_30.mtx"; Eigen::SparseMatrix A; - loadSymmetric(A, path + "/" + MatrixName); - std::cout << "Matrix Load OK" << std::endl; + // const bool ok = loadMarket(A, path + "/A_2.mat"); + // REQUIRE(ok); + loadSymmetric_test(A,path + "/A_2.mat"); Eigen::VectorXd b(A.rows()); b.setOnes(); Eigen::VectorXd x(b.size()); @@ -610,7 +760,7 @@ TEST_CASE("Hyprel_b2", "[solver]") start = clock(); auto solver = LinearSolver::create("Hypre", ""); json params; - params["Hypre"]["block_size"] = 2; + params["Hypre"]["block_size"] = 3; params["Hypre"]["tolerance"] = 1e-8; params["Hypre"]["max_iter"] = 1000; solver->setParameters(params); @@ -632,121 +782,39 @@ TEST_CASE("Hyprel_b2", "[solver]") } #endif -#ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("Hybre_crystm03", "[solver]") +TEST_CASE("amgcl_blocksolver_Bicgstab", "[solver]") { +#ifndef NDEBUG + return; +#endif + std::cout << "Polysolve AMGCL Solver" << std::endl; const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "crystm03.mtx"; Eigen::SparseMatrix A; loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; - Eigen::VectorXd b(A.rows()); - b.setOnes(); - Eigen::VectorXd x(b.size()); - x.setZero(); - Eigen::VectorXd x_b(b.size()); - x_b.setZero(); - { - clock_t start, end; - json solver_info; - start = clock(); - auto solver = LinearSolver::create("Hypre", ""); - json params; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x); - end = clock(); - solver->getInfo(solver_info); - std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; - } - { - clock_t start, end; - json solver_info; - start = clock(); - auto solver = LinearSolver::create("Hypre", ""); - json params; - params["Hypre"]["block_size"] = 3; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x_b); - end = clock(); - solver->getInfo(solver_info); - std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; - } - const double err = (A * x - b).norm() / b.norm(); - const double err_b = (A * x_b - b).norm() / b.norm(); - std::cout << "Scalar relative error " << err << std::endl; - std::cout << "Block relative error " << err_b << std::endl; - REQUIRE(err < 1e-8); - REQUIRE(err_b < 1e-8); -} -#endif -#ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("hypre_smallscale", "[solver]") -{ - const std::string path = POLYSOLVE_DATA_DIR; - Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); Eigen::VectorXd b(A.rows()); b.setOnes(); - Eigen::VectorXd x(b.size()); + Eigen::VectorXd x(A.rows()); x.setZero(); - Eigen::VectorXd x_b(b.size()); - x_b.setZero(); { - clock_t start, end; json solver_info; - start = clock(); - auto solver = LinearSolver::create("Hypre", ""); + auto solver = LinearSolver::create("AMGCL_cuda", ""); json params; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + params["AMGCL"]["solver_type"] = "cg"; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); solver->solve(b, x); - end = clock(); - solver->getInfo(solver_info); - std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; - } - { - clock_t start, end; - json solver_info; - start = clock(); - auto solver = LinearSolver::create("Hypre", ""); - json params; - params["Hypre"]["block_size"] = 3; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; - solver->setParameters(params); - solver->analyzePattern(A, A.rows()); - solver->factorize(A); - solver->solve(b, x_b); - end = clock(); solver->getInfo(solver_info); - std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + REQUIRE(solver_info["num_iterations"] > 0); std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } - const double err = (A * x - b).norm() / b.norm(); - const double err_b = (A * x_b - b).norm() / b.norm(); - std::cout << "Scalar relative error " << err << std::endl; - std::cout << "Block relative error " << err_b << std::endl; - REQUIRE(err < 1e-8); - REQUIRE(err_b < 1e-8); + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); } -#endif + From c215d8ca13f14954480a6dd268f042fc16e6deca Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Wed, 20 Jul 2022 11:24:31 -0400 Subject: [PATCH 14/29] Update AMGCL CUDA solver --- CMakeLists.txt | 49 ++- ...lysolve_target_link_system_libraries.cmake | 28 ++ cmake/polysolve/polysolve_warnings.cmake | 1 - cudatest/CMakeLists.txt | 10 + cudatest/cudatest.cu | 73 +++++ cudatest/main.cpp | 6 + src/polysolve/LinearSolver.cpp | 2 + src/polysolve/LinearSolverAMGCL_cuda.cu | 4 +- src/polysolve/LinearSolverAMGCL_cuda.hpp | 2 + tests/CMakeLists.txt | 1 - tests/test_solver.cpp | 282 +++++++++--------- 11 files changed, 313 insertions(+), 145 deletions(-) create mode 100644 cmake/polysolve/polysolve_target_link_system_libraries.cmake create mode 100644 cudatest/CMakeLists.txt create mode 100644 cudatest/cudatest.cu create mode 100644 cudatest/main.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 76f9405..e9f0f65 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,6 +24,7 @@ if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/PolySolveOptions.cmake) include(${CMAKE_CURRENT_SOURCE_DIR}/PolySolveOptions.cmake) endif() + ################################################################################ project(PolySolve @@ -82,6 +83,7 @@ include(polysolve_use_colors) # IPC Toolkit utils include(polysolve_prepend_current_path) include(polysolve_set_source_group) +include(polysolve_target_link_system_libraries) # Sort projects inside the solution set_property(GLOBAL PROPERTY USE_FOLDERS ON) @@ -134,17 +136,45 @@ include(eigen) target_link_libraries(polysolve PUBLIC Eigen3::Eigen) # CUDA + +set(CMAKE_CUDA_COMPILER "/usr/local/cuda/bin/nvcc") if(POLYSOLVE_WITH_CUDA) - find_package(CUDAToolkit) - if (CUDAToolkit_FOUND) - set_target_properties(polysolve PROPERTIES - CUDA_SEPARABLE_COMPILATION ON CUDA_ARCHITECTURES Pascal) - include_directories(${CUDAToolkit_INCLUDE_DIRS}) - link_directories(${CUDAToolkit_LIBRARY_DIR}) - target_link_libraries(polysolve INTERFACE CUDA::cusparse CUDA::cudart_static) + include(CheckLanguage) + check_language(CUDA) + set(CMAKE_CUDA_ARCHITECTURES "70;75;86") + if(CMAKE_CUDA_COMPILER) + target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) + enable_language(CUDA) + else() + message(FATAL_ERROR "No CUDA support found!") + endif() + # We need to explicitly state that we need all CUDA files in the particle + # library to be built with -dc as the member functions could be called by + # other libraries and executables. + set_target_properties(polysolve PROPERTIES CUDA_SEPARABLE_COMPILATION ON) + + # Nvidia RTX8000 -> compute_75 + # Nvidia V100 -> compute_70 + # Nvidia 1080/1080Ti -> compute_61 + # Nvidia 3080Ti -> compute_86 + if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) + set(CMAKE_CUDA_ARCHITECTURES 70 75 86) + endif() + set_target_properties(polysolve PROPERTIES CUDA_ARCHITECTURES "70;75;86") + + if(APPLE) + # We need to add the path to the driver (libcuda.dylib) as an rpath, + # so that the static cuda runtime can find it at runtime. + set_property(TARGET polysolve + PROPERTY + BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) endif() + find_package(CUDAToolkit) + polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cudart) + polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cusparse) endif() + # if(POLYSOLVE_WITH_CUDA) # find_package(CUDA) # if (CUDA_FOUND) @@ -168,6 +198,7 @@ endif() # target_link_libraries(polysolve INTERFACE ${CUDA_cusparse_LIBRARY}) # endif() # endif() + # Hypre (GNU Lesser General Public License) if(POLYSOLVE_WITH_HYPRE) include(hypre) @@ -267,3 +298,7 @@ if(POLYSOLVE_WITH_TESTS) add_subdirectory(tests) endif() + +if(POLYSOLVE_WITH_CUDA) + add_subdirectory(cudatest) +endif() diff --git a/cmake/polysolve/polysolve_target_link_system_libraries.cmake b/cmake/polysolve/polysolve_target_link_system_libraries.cmake new file mode 100644 index 0000000..43f8a14 --- /dev/null +++ b/cmake/polysolve/polysolve_target_link_system_libraries.cmake @@ -0,0 +1,28 @@ +function(polysolve_target_link_system_libraries target) + set(options PRIVATE PUBLIC INTERFACE) + cmake_parse_arguments(TLLS "${options}" "" "" ${ARGN}) + foreach(op ${options}) + if(TLLS_${op}) + set(scope ${op}) + endif() + endforeach(op) + set(libs ${TLLS_UNPARSED_ARGUMENTS}) + + foreach(lib ${libs}) + get_target_property(lib_include_dirs ${lib} INTERFACE_INCLUDE_DIRECTORIES) + if(lib_include_dirs) + if(scope) + target_include_directories(${target} SYSTEM ${scope} ${lib_include_dirs}) + else() + target_include_directories(${target} SYSTEM PRIVATE ${lib_include_dirs}) + endif() + else() + message(WARNING "${lib} doesn't set INTERFACE_INCLUDE_DIRECTORIES. No include_directories set.") + endif() + if(scope) + target_link_libraries(${target} ${scope} ${lib}) + else() + target_link_libraries(${target} ${lib}) + endif() + endforeach() +endfunction(polysolve_target_link_system_libraries) \ No newline at end of file diff --git a/cmake/polysolve/polysolve_warnings.cmake b/cmake/polysolve/polysolve_warnings.cmake index aafade8..16bb62d 100644 --- a/cmake/polysolve/polysolve_warnings.cmake +++ b/cmake/polysolve/polysolve_warnings.cmake @@ -10,7 +10,6 @@ endif() set(POLYSOLVE_FLAGS -Wall -Wextra - -pedantic # -Wconversion #-Wunsafe-loop-optimizations # broken with C++11 loops diff --git a/cudatest/CMakeLists.txt b/cudatest/CMakeLists.txt new file mode 100644 index 0000000..bd85295 --- /dev/null +++ b/cudatest/CMakeLists.txt @@ -0,0 +1,10 @@ +################################################################################ +# Tests +################################################################################ +add_executable(cudatest main.cpp cudatest.cu ) +target_link_libraries(cudatest PUBLIC polysolve::polysolve) + +include(catch2) +target_link_libraries(cudatest PUBLIC Catch2::Catch2) +set(PARSE_CATCH_TESTS_ADD_TO_CONFIGURE_DEPENDS ON) +catch_discover_tests(cudatest) \ No newline at end of file diff --git a/cudatest/cudatest.cu b/cudatest/cudatest.cu new file mode 100644 index 0000000..da145c5 --- /dev/null +++ b/cudatest/cudatest.cu @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#include +#include +#include +////////////////////////////////////////////////////////////////////////// +#include + +using namespace polysolve; + +void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) +{ + std::ifstream fin(PATH); + long int M, N, L; + while (fin.peek() == '%') + { + fin.ignore(2048, '\n'); + } + fin >> M >> N >> L; + A.resize(M, N); + A.reserve(L * 2 - M); + std::vector> triple; + for (size_t i = 0; i < L; i++) + { + int m, n; + double data; + fin >> m >> n >> data; + triple.push_back(Eigen::Triplet(m - 1, n - 1, data)); + if (m != n) + { + triple.push_back(Eigen::Triplet(n - 1, m - 1, data)); + } + } + fin.close(); + A.setFromTriplets(triple.begin(), triple.end()); +}; + +TEST_CASE("AMGCL_GPU", "[solver]"){ + const std::string path = "/home/yiwei/polysolve/tests/data"; + std::string MatrixName = "Serena.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + + std::cout << "Matrix Load OK" << std::endl; + + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(A.rows()); + x.setZero(); + { + amgcl::profiler<> prof("Serena_GPU"); + json solver_info; + auto solver = LinearSolver::create("AMGCL_cuda", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + params["AMGCL"]["solver_type"] = "cg"; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout<< prof << std::endl; + } + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +} diff --git a/cudatest/main.cpp b/cudatest/main.cpp new file mode 100644 index 0000000..f30959d --- /dev/null +++ b/cudatest/main.cpp @@ -0,0 +1,6 @@ +//////////////////////////////////////////////////////////////////////////////// +// Keep this file empty, and implement unit tests in separate compilation units! +//////////////////////////////////////////////////////////////////////////////// + +#define CATCH_CONFIG_MAIN +#include \ No newline at end of file diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp index dd3bdaf..b4bd18b 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/LinearSolver.cpp @@ -25,8 +25,10 @@ #endif #ifdef POLYSOLVE_WITH_AMGCL #include +#ifdef POLYSOLVE_WITH_CUDA #include #endif +#endif #include //////////////////////////////////////////////////////////////////////////////// diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu index e768eaf..5b09387 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.cu +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -1,7 +1,8 @@ #ifdef POLYSOLVE_WITH_AMGCL +#ifdef POLYSOLVE_WITH_CUDA //////////////////////////////////////////////////////////////////////////////// -#include +#include #include #include @@ -173,3 +174,4 @@ namespace polysolve } // namespace polysolve #endif +#endif diff --git a/src/polysolve/LinearSolverAMGCL_cuda.hpp b/src/polysolve/LinearSolverAMGCL_cuda.hpp index 5d9593b..b3d0dce 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.hpp +++ b/src/polysolve/LinearSolverAMGCL_cuda.hpp @@ -1,6 +1,7 @@ #pragma once #ifdef POLYSOLVE_WITH_AMGCL +#ifdef POLYSOLVE_WITH_CUDA //////////////////////////////////////////////////////////////////////////////// #include @@ -103,3 +104,4 @@ namespace polysolve } // namespace polysolve #endif +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f9e005f..226578e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -7,7 +7,6 @@ set(test_sources test_solver.cpp ) add_executable(unit_tests ${test_sources}) - ################################################################################ # Required Libraries ################################################################################ diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index fcbb93c..50b41c4 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include ////////////////////////////////////////////////////////////////////////// using namespace polysolve; @@ -602,126 +602,6 @@ void loadSymmetric_test(Eigen::SparseMatrix &A, std::string PATH) // } // #endif -// #ifdef POLYSOLVE_WITH_HYPRE -// TEST_CASE("Hypre_b2", "[solver]") -// { -// const std::string path = POLYSOLVE_DATA_DIR; -// std::string MatrixName = "gr_30_30.mtx"; -// Eigen::SparseMatrix A; -// loadSymmetric(A, path + "/" + MatrixName); -// std::cout << "Matrix Load OK" << std::endl; -// Eigen::VectorXd b(A.rows()); -// b.setOnes(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); -// Eigen::VectorXd x_b(b.size()); -// x_b.setZero(); -// { -// clock_t start, end; -// json solver_info; -// start = clock(); -// auto solver = LinearSolver::create("Hypre", ""); -// json params; -// params["Hypre"]["tolerance"] = 1e-8; -// params["Hypre"]["max_iter"] = 1000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); -// end = clock(); -// solver->getInfo(solver_info); -// std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl; -// } -// { -// clock_t start, end; -// json solver_info; -// start = clock(); -// auto solver = LinearSolver::create("Hypre", ""); -// json params; -// params["Hypre"]["block_size"] = 2; -// params["Hypre"]["tolerance"] = 1e-8; -// params["Hypre"]["max_iter"] = 1000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x_b); -// end = clock(); -// solver->getInfo(solver_info); -// std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl; -// } -// const double err = (A * x - b).norm() / b.norm(); -// const double err_b = (A * x_b - b).norm() / b.norm(); -// std::cout << "Scalar relative error " << err << std::endl; -// std::cout << "Block relative error " << err_b << std::endl; -// REQUIRE(err < 1e-8); -// REQUIRE(err_b < 1e-8); -// } -// #endif - -// #ifdef POLYSOLVE_WITH_HYPRE -// TEST_CASE("Hypre_crystm03", "[solver]") -// { -// const std::string path = POLYSOLVE_DATA_DIR; -// std::string MatrixName = "crystm03.mtx"; -// Eigen::SparseMatrix A; -// loadSymmetric(A, path + "/" + MatrixName); -// std::cout << "Matrix Load OK" << std::endl; -// Eigen::VectorXd b(A.rows()); -// b.setOnes(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); -// Eigen::VectorXd x_b(b.size()); -// x_b.setZero(); -// { -// clock_t start, end; -// json solver_info; -// start = clock(); -// auto solver = LinearSolver::create("Hypre", ""); -// json params; -// params["Hypre"]["tolerance"] = 1e-8; -// params["Hypre"]["max_iter"] = 1000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); -// end = clock(); -// solver->getInfo(solver_info); -// std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl; -// } -// { -// clock_t start, end; -// json solver_info; -// start = clock(); -// auto solver = LinearSolver::create("Hypre", ""); -// json params; -// params["Hypre"]["block_size"] = 3; -// params["Hypre"]["tolerance"] = 1e-8; -// params["Hypre"]["max_iter"] = 1000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x_b); -// end = clock(); -// solver->getInfo(solver_info); -// std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl; -// } -// const double err = (A * x - b).norm() / b.norm(); -// const double err_b = (A * x_b - b).norm() / b.norm(); -// std::cout << "Scalar relative error " << err << std::endl; -// std::cout << "Block relative error " << err_b << std::endl; -// REQUIRE(err < 1e-8); -// REQUIRE(err_b < 1e-8); -// } -// #endif - #ifdef POLYSOLVE_WITH_HYPRE TEST_CASE("hypre_smallscale", "[solver]") { @@ -750,7 +630,6 @@ TEST_CASE("hypre_smallscale", "[solver]") solver->solve(b, x); end = clock(); solver->getInfo(solver_info); - std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } @@ -769,52 +648,185 @@ TEST_CASE("hypre_smallscale", "[solver]") solver->solve(b, x_b); end = clock(); solver->getInfo(solver_info); - std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } const double err = (A * x - b).norm() / b.norm(); const double err_b = (A * x_b - b).norm() / b.norm(); - std::cout << "Scalar relative error " << err << std::endl; - std::cout << "Block relative error " << err_b << std::endl; REQUIRE(err < 1e-8); REQUIRE(err_b < 1e-8); } #endif -TEST_CASE("amgcl_blocksolver_Bicgstab", "[solver]") +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("Hypre_b2", "[solver]") { -#ifndef NDEBUG - return; + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "gr_30_30.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(b.size()); + x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + end = clock(); + solver->getInfo(solver_info); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["block_size"] = 2; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + end = clock(); + solver->getInfo(solver_info); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); +} #endif - std::cout << "Polysolve AMGCL Solver" << std::endl; + +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("Hypre_crystm03", "[solver]") +{ const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "crystm03.mtx"; Eigen::SparseMatrix A; loadSymmetric(A, path + "/" + MatrixName); - std::cout << "Matrix Load OK" << std::endl; + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(b.size()); + x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + end = clock(); + solver->getInfo(solver_info); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + { + clock_t start, end; + json solver_info; + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + json params; + params["Hypre"]["block_size"] = 3; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + end = clock(); + solver->getInfo(solver_info); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; + } + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); +} +#endif +#ifdef POLYSOLVE_WITH_AMGCL +TEST_CASE("AMGCL_CPU", "[solver]") +{ +#ifndef NDEBUG + return; +#endif + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "Serena.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; Eigen::VectorXd b(A.rows()); b.setOnes(); + Eigen::VectorXd x_b(A.rows()); + x_b.setZero(); Eigen::VectorXd x(A.rows()); x.setZero(); { + amgcl::profiler<> prof("Serena_Block"); json solver_info; - auto solver = LinearSolver::create("AMGCL_cuda", ""); + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 1000; + params["AMGCL"]["block_size"] = 3; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x_b); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout<< prof << std::endl; + } + { + amgcl::profiler<> prof("Serena_Scalar"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); json params; params["AMGCL"]["tolerance"] = 1e-8; params["AMGCL"]["max_iter"] = 10000; - params["AMGCL"]["solver_type"] = "cg"; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); solver->solve(b, x); + prof.toc("solve"); solver->getInfo(solver_info); REQUIRE(solver_info["num_iterations"] > 0); - std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; + std::cout << prof << std::endl; } REQUIRE((A * x - b).norm() / b.norm() < 1e-7); + REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); } - +#endif From c05c48be4c60663d00f223ca4e71053f775fabd9 Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Wed, 20 Jul 2022 23:58:09 -0400 Subject: [PATCH 15/29] Update AMGCL CUDA solver and Hypre block solver --- CMakeLists.txt | 38 +- cmake/polysolve/polysolve_warnings.cmake | 1 + cmake/recipes/hypre.cmake | 4 +- cmake/recipes/json.cmake | 6 +- cudatest/CMakeLists.txt | 33 +- cudatest/cudatest.cu | 43 +- src/CMakeLists.txt | 7 +- tests/test_solver.cpp | 1160 ++++++++++------------ 8 files changed, 623 insertions(+), 669 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e9f0f65..47b27f2 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -136,8 +136,8 @@ include(eigen) target_link_libraries(polysolve PUBLIC Eigen3::Eigen) # CUDA - set(CMAKE_CUDA_COMPILER "/usr/local/cuda/bin/nvcc") + if(POLYSOLVE_WITH_CUDA) include(CheckLanguage) check_language(CUDA) @@ -174,31 +174,6 @@ if(POLYSOLVE_WITH_CUDA) polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cusparse) endif() - -# if(POLYSOLVE_WITH_CUDA) -# find_package(CUDA) -# if (CUDA_FOUND) -# set(CUDA_TARGET_ARCH "Ampere" CACHE STRING "Target architecture(s) for CUDA") -# cuda_select_nvcc_arch_flags(CUDA_ARCH_FLAGS ${CUDA_TARGET_ARCH}) - -# if (OPENMP_FOUND) -# list(APPEND CUDA_NVCC_FLAGS -Xcompiler ${OpenMP_CXX_FLAGS}) -# endif() - -# if (CMAKE_CXX_COMPILER_ID MATCHES "GNU") -# list(APPEND CUDA_NVCC_FLAGS -# ${CUDA_ARCH_FLAGS} -Wno-deprecated-gpu-targets) - -# list(APPEND CUDA_NVCC_FLAGS -Xcompiler -fPIC -Xcompiler -Wno-vla) -# endif() - -# # add_library(cuda_target INTERFACE) -# include_directories(${CUDA_INCLUDE_DIRS}) -# link_directories(${CUDA_LIBRARY_DIR}) -# target_link_libraries(polysolve INTERFACE ${CUDA_cusparse_LIBRARY}) -# endif() -# endif() - # Hypre (GNU Lesser General Public License) if(POLYSOLVE_WITH_HYPRE) include(hypre) @@ -211,7 +186,7 @@ endif() # Json (MIT) include(json) -target_link_libraries(polysolve PUBLIC nlohmann::json) +target_link_libraries(polysolve PUBLIC nlohmann_json::nlohmann_json) # CHOLMOD solver if(POLYSOLVE_WITH_CHOLMOD) @@ -297,8 +272,11 @@ if(POLYSOLVE_WITH_TESTS) include("${catch2_SOURCE_DIR}/contrib/Catch.cmake") add_subdirectory(tests) + + # Cuda test + if(POLYSOLVE_WITH_CUDA) + add_subdirectory(cudatest) + endif() endif() -if(POLYSOLVE_WITH_CUDA) - add_subdirectory(cudatest) -endif() + diff --git a/cmake/polysolve/polysolve_warnings.cmake b/cmake/polysolve/polysolve_warnings.cmake index 16bb62d..aafade8 100644 --- a/cmake/polysolve/polysolve_warnings.cmake +++ b/cmake/polysolve/polysolve_warnings.cmake @@ -10,6 +10,7 @@ endif() set(POLYSOLVE_FLAGS -Wall -Wextra + -pedantic # -Wconversion #-Wunsafe-loop-optimizations # broken with C++11 loops diff --git a/cmake/recipes/hypre.cmake b/cmake/recipes/hypre.cmake index 904d1f7..f0fd537 100644 --- a/cmake/recipes/hypre.cmake +++ b/cmake/recipes/hypre.cmake @@ -6,11 +6,11 @@ endif() message(STATUS "Third-party: creating target 'HYPRE::HYPRE'") -set(HYPRE_SEQUENTIAL OFF CACHE INTERNAL "" FORCE) +set(HYPRE_SEQUENTIAL ON CACHE INTERNAL "" FORCE) set(HYPRE_PRINT_ERRORS ON CACHE INTERNAL "" FORCE) set(HYPRE_BIGINT ON CACHE INTERNAL "" FORCE) set(HYPRE_USING_FEI OFF CACHE INTERNAL "" FORCE) -set(HYPRE_USING_OPENMP ON CACHE INTERNAL "" FORCE) +set(HYPRE_USING_OPENMP OFF CACHE INTERNAL "" FORCE) set(HYPRE_SHARED OFF CACHE INTERNAL "" FORCE) diff --git a/cmake/recipes/json.cmake b/cmake/recipes/json.cmake index 8aa5b09..38c7619 100644 --- a/cmake/recipes/json.cmake +++ b/cmake/recipes/json.cmake @@ -11,11 +11,11 @@ # # JSON MIT -if(TARGET nlohmann::json) +if(TARGET nlohmann_json::nlohmann_json) return() endif() -message(STATUS "Third-party: creating target 'nlohmann_json'") +message(STATUS "Third-party: creating target 'nlohmann_json::nlohmann_json'") # nlohmann_json is a big repo for a single header, so we just download the release archive set(NLOHMANNJSON_VERSION "v3.10.2") @@ -29,7 +29,7 @@ FetchContent_Declare( FetchContent_MakeAvailable(nlohmann_json) add_library(nlohmann_json INTERFACE) -add_library(nlohmann::json ALIAS nlohmann_json) +add_library(nlohmann_json::nlohmann_json ALIAS nlohmann_json) include(GNUInstallDirs) target_include_directories(nlohmann_json INTERFACE diff --git a/cudatest/CMakeLists.txt b/cudatest/CMakeLists.txt index bd85295..0ad97d0 100644 --- a/cudatest/CMakeLists.txt +++ b/cudatest/CMakeLists.txt @@ -1,10 +1,35 @@ ################################################################################ # Tests ################################################################################ -add_executable(cudatest main.cpp cudatest.cu ) -target_link_libraries(cudatest PUBLIC polysolve::polysolve) +set(test_sources + main.cpp + cudatest.cu +) +add_executable(cuda_test ${test_sources}) +################################################################################ +# Required Libraries +################################################################################ include(catch2) -target_link_libraries(cudatest PUBLIC Catch2::Catch2) +target_link_libraries(cuda_test PUBLIC Catch2::Catch2) + +target_link_libraries(cuda_test PUBLIC polysolve::polysolve) + +include(polysolve_warnings) +target_link_libraries(cuda_test PRIVATE polysolve::warnings) +################################################################################ +# Register tests +################################################################################ +foreach(source IN ITEMS ${test_sources}) + source_group("tests" FILES "${source}") +endforeach() + +# Register tests set(PARSE_CATCH_TESTS_ADD_TO_CONFIGURE_DEPENDS ON) -catch_discover_tests(cudatest) \ No newline at end of file +catch_discover_tests(cuda_test) + +################################################################################ +# Data +################################################################################ +set(DATA_DIR "${CMAKE_SOURCE_DIR}/tests/data") +target_compile_definitions(cuda_test PUBLIC -DPOLYSOLVE_DATA_DIR=\"${DATA_DIR}\") \ No newline at end of file diff --git a/cudatest/cudatest.cu b/cudatest/cudatest.cu index da145c5..e354725 100644 --- a/cudatest/cudatest.cu +++ b/cudatest/cudatest.cu @@ -37,9 +37,9 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) A.setFromTriplets(triple.begin(), triple.end()); }; -TEST_CASE("AMGCL_GPU", "[solver]"){ - const std::string path = "/home/yiwei/polysolve/tests/data"; - std::string MatrixName = "Serena.mtx"; +TEST_CASE("amgcl_crystm03_cg", "[solver]"){ + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "crystm03.mtx"; Eigen::SparseMatrix A; loadSymmetric(A, path + "/" + MatrixName); @@ -50,7 +50,7 @@ TEST_CASE("AMGCL_GPU", "[solver]"){ Eigen::VectorXd x(A.rows()); x.setZero(); { - amgcl::profiler<> prof("Serena_GPU"); + amgcl::profiler<> prof("crystm03_GPU"); json solver_info; auto solver = LinearSolver::create("AMGCL_cuda", ""); prof.tic("setup"); @@ -71,3 +71,38 @@ TEST_CASE("AMGCL_GPU", "[solver]"){ } REQUIRE((A * x - b).norm() / b.norm() < 1e-7); } + +TEST_CASE("amgcl_crystm03_bicgstab", "[solver]"){ + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "crystm03.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + + std::cout << "Matrix Load OK" << std::endl; + + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x(A.rows()); + x.setZero(); + { + amgcl::profiler<> prof("crystm03_GPU"); + json solver_info; + auto solver = LinearSolver::create("AMGCL_cuda", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + params["AMGCL"]["solver_type"] = "bicgstab"; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout<< prof << std::endl; + } + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 59ba6dc..20283cf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,14 +17,9 @@ set(SOURCES polysolve/LinearSolverAMGCL_cuda.cu polysolve/LinearSolverAMGCL_cuda.hpp ) -# set(CUDA_SOURCES -# polysolve/LinearSolverAMGCL_cuda.cu -# polysolve/LinearSolverAMGCL_cuda.hpp -# ) polysolve_prepend_current_path(SOURCES) polysolve_set_source_group(${SOURCES}) -# polysolve_prepend_current_path(CUDA_SOURCES) -# polysolve_set_source_group(${CUDA_SOURCES}) + target_sources(polysolve PRIVATE ${SOURCES}) diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index 50b41c4..252de03 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -37,629 +37,543 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) fin.close(); A.setFromTriplets(triple.begin(), triple.end()); }; -void loadSymmetric_test(Eigen::SparseMatrix &A, std::string PATH) +TEST_CASE("all", "[solver]") { - std::ifstream fin(PATH); - long int M, N, L; - while (fin.peek() == '%') + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solvers = LinearSolver::availableSolvers(); + + for (const auto &s : solvers) { - fin.ignore(2048, '\n'); + if (s == "Eigen::DGMRES") + continue; +#ifdef WIN32 + if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") + continue; +#endif + auto solver = LinearSolver::create(s, ""); + json params; + params[s]["tolerance"] = 1e-10; + solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // solver->getInfo(solver_info); + + // std::cout<<"Solver error: "<> M >> N >> L; - A.resize(M, N); - A.reserve(L * 2 - M); - std::vector> triple; - for (size_t i = 0; i < L; i++) +} + +TEST_CASE("eigen_params", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solvers = LinearSolver::availableSolvers(); + + for (const auto &s : solvers) { - int m, n; - double data; - fin >> m >> n >> data; - if (m != n) + if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES" || s == "Eigen::LeastSquaresConjugateGradient" || s == "Eigen::DGMRES") { - triple.push_back(Eigen::Triplet(m - 1, n - 1, data/2)); - triple.push_back(Eigen::Triplet(n - 1, m - 1, data/2)); + auto solver = LinearSolver::create(s, ""); + json params; + params[s]["max_iter"] = 1000; + params[s]["tolerance"] = 1e-10; + solver->setParameters(params); + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // solver->getInfo(solver_info); + + // std::cout<<"Solver error: "<(m - 1, n - 1, data)); + } +} + +TEST_CASE("pre_factor", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solvers = LinearSolver::availableSolvers(); + + for (const auto &s : solvers) + { + if (s == "Eigen::DGMRES") + continue; +#ifdef WIN32 + if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") + continue; +#endif + auto solver = LinearSolver::create(s, ""); + solver->analyzePattern(A, A.rows()); + + std::default_random_engine eng{42}; + std::uniform_real_distribution urd(0.1, 5); + + for (int i = 0; i < 10; ++i) + { + std::vector> tripletList; + + for (int k = 0; k < A.outerSize(); ++k) + { + for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) + { + if (it.row() == it.col()) + { + tripletList.emplace_back(it.row(), it.col(), urd(eng) * 100); + } + else if (it.row() < it.col()) + { + const double val = -urd(eng); + tripletList.emplace_back(it.row(), it.col(), val); + tripletList.emplace_back(it.col(), it.row(), val); + } + } + } + + Eigen::SparseMatrix Atmp(A.rows(), A.cols()); + Atmp.setFromTriplets(tripletList.begin(), tripletList.end()); + + Eigen::VectorXd b(Atmp.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->factorize(Atmp); + solver->solve(b, x); + + // solver->getInfo(solver_info); + + // std::cout<<"Solver error: "< A; -// const bool ok = loadMarket(A, path + "/A_2.mat"); -// REQUIRE(ok); - -// auto solvers = LinearSolver::availableSolvers(); - -// for (const auto &s : solvers) -// { -// if (s == "Eigen::DGMRES") -// continue; -// #ifdef WIN32 -// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") -// continue; -// #endif -// auto solver = LinearSolver::create(s, ""); -// json params; -// params[s]["tolerance"] = 1e-10; -// solver->setParameters(params); -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); - -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); - -// // solver->getInfo(solver_info); - -// // std::cout<<"Solver error: "< A; -// const bool ok = loadMarket(A, path + "/A_2.mat"); -// REQUIRE(ok); - -// auto solvers = LinearSolver::availableSolvers(); - -// for (const auto &s : solvers) -// { -// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES" || s == "Eigen::LeastSquaresConjugateGradient" || s == "Eigen::DGMRES") -// { -// auto solver = LinearSolver::create(s, ""); -// json params; -// params[s]["max_iter"] = 1000; -// params[s]["tolerance"] = 1e-10; -// solver->setParameters(params); - -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); - -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); - -// // solver->getInfo(solver_info); - -// // std::cout<<"Solver error: "< A; -// const bool ok = loadMarket(A, path + "/A_2.mat"); -// REQUIRE(ok); - -// auto solvers = LinearSolver::availableSolvers(); - -// for (const auto &s : solvers) -// { -// if (s == "Eigen::DGMRES") -// continue; -// #ifdef WIN32 -// if (s == "Eigen::ConjugateGradient" || s == "Eigen::BiCGSTAB" || s == "Eigen::GMRES" || s == "Eigen::MINRES") -// continue; -// #endif -// auto solver = LinearSolver::create(s, ""); -// solver->analyzePattern(A, A.rows()); - -// std::default_random_engine eng{42}; -// std::uniform_real_distribution urd(0.1, 5); - -// for (int i = 0; i < 10; ++i) -// { -// std::vector> tripletList; - -// for (int k = 0; k < A.outerSize(); ++k) -// { -// for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) -// { -// if (it.row() == it.col()) -// { -// tripletList.emplace_back(it.row(), it.col(), urd(eng) * 100); -// } -// else if (it.row() < it.col()) -// { -// const double val = -urd(eng); -// tripletList.emplace_back(it.row(), it.col(), val); -// tripletList.emplace_back(it.col(), it.row(), val); -// } -// } -// } - -// Eigen::SparseMatrix Atmp(A.rows(), A.cols()); -// Atmp.setFromTriplets(tripletList.begin(), tripletList.end()); - -// Eigen::VectorXd b(Atmp.rows()); -// b.setRandom(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); - -// solver->factorize(Atmp); -// solver->solve(b, x); - -// // solver->getInfo(solver_info); - -// // std::cout<<"Solver error: "< A; -// // const bool ok = loadMarket(A, path + "/A_2.mat"); -// // REQUIRE(ok); -// loadSymmetric_test(A,path + "/A_2.mat"); -// auto solver = LinearSolver::create("Hypre", ""); -// // solver->setParameters(params); -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(b.size()); -// x.setZero(); - -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); - -// // solver->getInfo(solver_info); - -// // std::cout<<"Solver error: "< A; -// // const bool ok = loadMarket(A, path + "/A_2.mat"); -// // REQUIRE(ok); - -// loadSymmetric_test(A,path + "/A_2.mat"); -// // solver->setParameters(params); -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(A.rows()); -// x.setZero(); -// { -// json solver_info; - -// auto solver = LinearSolver::create("Hypre", ""); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] > 1); -// } - -// { -// json solver_info; -// auto solver = LinearSolver::create("Hypre", ""); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); - -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] == 1); -// } - -// // std::cout<<"Solver error: "< A; -// const bool ok = loadMarket(A, path + "/A_2.mat"); -// REQUIRE(ok); - -// // solver->setParameters(params); -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(A.rows()); -// x.setZero(); -// { -// json solver_info; - -// auto solver = LinearSolver::create("AMGCL", ""); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] > 0); -// } - -// { -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); - -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] == 0); -// } - -// // std::cout<<"Solver error: "< A; -// bool ok = loadMarket(A, path + "/A0.mat"); -// REQUIRE(ok); - -// Eigen::VectorXd b; -// ok = loadMarketVector(b, path + "/b0.mat"); -// REQUIRE(ok); - -// auto solver = LinearSolver::create("SaddlePointSolver", ""); -// solver->analyzePattern(A, 9934); -// solver->factorize(A); -// Eigen::VectorXd x(A.rows()); -// solver->solve(b, x); -// const double err = (A * x - b).norm(); -// REQUIRE(err < 1e-8); -// } - -// #ifdef POLYSOLVE_WITH_AMGCL -// TEST_CASE("amgcl_blocksolver_small_scale", "[solver]") -// { -// #ifndef NDEBUG -// return; -// #endif -// const std::string path = POLYSOLVE_DATA_DIR; -// Eigen::SparseMatrix A; -// const bool ok = loadMarket(A, path + "/A_2.mat"); -// REQUIRE(ok); - -// // solver->setParameters(params); -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(A.rows()); -// Eigen::VectorXd x_b(A.rows()); -// x.setZero(); -// x_b.setZero(); -// { -// json solver_info; - -// auto solver = LinearSolver::create("AMGCL", ""); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x); -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] > 0); -// const double err = (A * x - b).norm(); -// REQUIRE(err < 1e-5); -// } - -// { -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["block_size"] = 3; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// solver->solve(b, x_b); -// solver->getInfo(solver_info); - -// REQUIRE(solver_info["num_iterations"] > 0); -// const double err = (A * x_b - b).norm(); -// REQUIRE(err < 1e-5); -// } -// } -// #endif - -// #ifdef POLYSOLVE_WITH_AMGCL -// TEST_CASE("amgcl_blocksolver_b2", "[solver]") -// { -// #ifndef NDEBUG -// return; -// #endif -// const std::string path = POLYSOLVE_DATA_DIR; -// std::string MatrixName = "gr_30_30.mtx"; -// Eigen::SparseMatrix A; -// loadSymmetric(A, path + "/" + MatrixName); - -// std::cout << "Matrix Load OK" << std::endl; - -// Eigen::VectorXd b(A.rows()); -// b.setRandom(); -// Eigen::VectorXd x(A.rows()); -// Eigen::VectorXd x_b(A.rows()); -// x.setOnes(); -// x_b.setOnes(); -// { -// amgcl::profiler<> prof("gr_30_30_Scalar"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 1000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// { -// amgcl::profiler<> prof("gr_30_30_Block"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 1000; -// params["AMGCL"]["block_size"] = 2; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x_b); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); -// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -// } -// #endif - -// #ifdef POLYSOLVE_WITH_AMGCL -// TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]") -// { -// #ifndef NDEBUG -// return; -// #endif -// std::cout << "Polysolve AMGCL Solver" << std::endl; -// const std::string path = POLYSOLVE_DATA_DIR; -// std::string MatrixName = "crystm03.mtx"; -// Eigen::SparseMatrix A; -// loadSymmetric(A, path + "/" + MatrixName); -// std::cout << "Matrix Load OK" << std::endl; -// Eigen::VectorXd b(A.rows()); -// b.setOnes(); -// Eigen::VectorXd x_b(A.rows()); -// x_b.setZero(); -// Eigen::VectorXd x(A.rows()); -// x.setZero(); -// { -// amgcl::profiler<> prof("crystm03_Block"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 1000; -// params["AMGCL"]["block_size"] = 3; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x_b); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// { -// amgcl::profiler<> prof("crystm03_Scalar"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 10000; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); -// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -// } -// #endif - -// #ifdef POLYSOLVE_WITH_AMGCL -// TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") -// { -// #ifndef NDEBUG -// return; -// #endif -// std::cout << "Polysolve AMGCL Solver" << std::endl; -// const std::string path = POLYSOLVE_DATA_DIR; -// std::string MatrixName = "crystm03.mtx"; -// Eigen::SparseMatrix A; -// loadSymmetric(A, path + "/" + MatrixName); - -// std::cout << "Matrix Load OK" << std::endl; - -// Eigen::VectorXd b(A.rows()); -// b.setOnes(); -// Eigen::VectorXd x_b(A.rows()); -// x_b.setZero(); -// Eigen::VectorXd x(A.rows()); -// x.setZero(); -// { -// amgcl::profiler<> prof("crystm03_Block"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 10000; -// params["AMGCL"]["block_size"] = 3; -// params["AMGCL"]["solver_type"] = "bicgstab"; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x_b); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// { -// amgcl::profiler<> prof("crystm03_Scalar"); -// json solver_info; -// auto solver = LinearSolver::create("AMGCL", ""); -// prof.tic("setup"); -// json params; -// params["AMGCL"]["tolerance"] = 1e-8; -// params["AMGCL"]["max_iter"] = 10000; -// params["AMGCL"]["solver_type"] = "bicgstab"; -// solver->setParameters(params); -// solver->analyzePattern(A, A.rows()); -// solver->factorize(A); -// prof.toc("setup"); -// prof.tic("solve"); -// solver->solve(b, x); -// prof.toc("solve"); -// solver->getInfo(solver_info); -// REQUIRE(solver_info["num_iterations"] > 0); -// std::cout << solver_info["num_iterations"] << std::endl; -// std::cout << solver_info["final_res_norm"] << std::endl -// << prof << std::endl; -// } -// REQUIRE((A * x - b).norm() / b.norm() < 1e-7); -// REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); -// } -// #endif +} #ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("hypre_smallscale", "[solver]") +TEST_CASE("hypre", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; Eigen::SparseMatrix A; - // const bool ok = loadMarket(A, path + "/A_2.mat"); - // REQUIRE(ok); - loadSymmetric_test(A,path + "/A_2.mat"); + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solver = LinearSolver::create("Hypre", ""); + // solver->setParameters(params); Eigen::VectorXd b(A.rows()); - b.setOnes(); + b.setRandom(); Eigen::VectorXd x(b.size()); x.setZero(); - Eigen::VectorXd x_b(b.size()); - x_b.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // solver->getInfo(solver_info); + + // std::cout<<"Solver error: "< A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(A.rows()); + x.setZero(); { - clock_t start, end; json solver_info; - start = clock(); + auto solver = LinearSolver::create("Hypre", ""); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] > 1); + } + + { + json solver_info; + auto solver = LinearSolver::create("Hypre", ""); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] == 1); + } + + // std::cout<<"Solver error: "< A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(A.rows()); + x.setZero(); + { + json solver_info; + + auto solver = LinearSolver::create("AMGCL", ""); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] > 0); + } + + { + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] == 0); + } + + // std::cout<<"Solver error: "< A; + bool ok = loadMarket(A, path + "/A0.mat"); + REQUIRE(ok); + + Eigen::VectorXd b; + ok = loadMarketVector(b, path + "/b0.mat"); + REQUIRE(ok); + + auto solver = LinearSolver::create("SaddlePointSolver", ""); + solver->analyzePattern(A, 9934); + solver->factorize(A); + Eigen::VectorXd x(A.rows()); + solver->solve(b, x); + const double err = (A * x - b).norm(); + REQUIRE(err < 1e-8); +} + +#ifdef POLYSOLVE_WITH_AMGCL +TEST_CASE("amgcl_blocksolver_small_scale", "[solver]") +{ +#ifndef NDEBUG + return; +#endif + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(A.rows()); + Eigen::VectorXd x_b(A.rows()); + x.setZero(); + x_b.setZero(); + { + json solver_info; + + auto solver = LinearSolver::create("AMGCL", ""); json params; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; + params["AMGCL"]["tolerance"] = 1e-8; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); solver->solve(b, x); - end = clock(); solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] > 0); + const double err = (A * x - b).norm(); + REQUIRE(err < 1e-5); + } + + { + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["block_size"] = 3; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x_b); + solver->getInfo(solver_info); + + REQUIRE(solver_info["num_iterations"] > 0); + const double err = (A * x_b - b).norm(); + REQUIRE(err < 1e-5); + } +} +#endif + +#ifdef POLYSOLVE_WITH_AMGCL +TEST_CASE("amgcl_blocksolver_b2", "[solver]") +{ +#ifndef NDEBUG + return; +#endif + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "gr_30_30.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + + std::cout << "Matrix Load OK" << std::endl; + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(A.rows()); + Eigen::VectorXd x_b(A.rows()); + x.setOnes(); + x_b.setOnes(); + { + amgcl::profiler<> prof("gr_30_30_Scalar"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 1000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; } { - clock_t start, end; + amgcl::profiler<> prof("gr_30_30_Block"); json solver_info; - start = clock(); - auto solver = LinearSolver::create("Hypre", ""); + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); json params; - params["Hypre"]["block_size"] = 3; - params["Hypre"]["tolerance"] = 1e-8; - params["Hypre"]["max_iter"] = 1000; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 1000; + params["AMGCL"]["block_size"] = 2; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); solver->solve(b, x_b); - end = clock(); + prof.toc("solve"); solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); std::cout << solver_info["num_iterations"] << std::endl; - std::cout << solver_info["final_res_norm"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; } - const double err = (A * x - b).norm() / b.norm(); - const double err_b = (A * x_b - b).norm() / b.norm(); - REQUIRE(err < 1e-8); - REQUIRE(err_b < 1e-8); + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); + REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); +} +#endif + +#ifdef POLYSOLVE_WITH_AMGCL +TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]") +{ +#ifndef NDEBUG + return; +#endif + std::cout << "Polysolve AMGCL Solver" << std::endl; + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "crystm03.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + std::cout << "Matrix Load OK" << std::endl; + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x_b(A.rows()); + x_b.setZero(); + Eigen::VectorXd x(A.rows()); + x.setZero(); + { + amgcl::profiler<> prof("crystm03_Block"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 1000; + params["AMGCL"]["block_size"] = 3; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x_b); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; + } + { + amgcl::profiler<> prof("crystm03_Scalar"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; + } + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); + REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); +} +#endif + +#ifdef POLYSOLVE_WITH_AMGCL +TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") +{ +#ifndef NDEBUG + return; +#endif + std::cout << "Polysolve AMGCL Solver" << std::endl; + const std::string path = POLYSOLVE_DATA_DIR; + std::string MatrixName = "crystm03.mtx"; + Eigen::SparseMatrix A; + loadSymmetric(A, path + "/" + MatrixName); + + std::cout << "Matrix Load OK" << std::endl; + + Eigen::VectorXd b(A.rows()); + b.setOnes(); + Eigen::VectorXd x_b(A.rows()); + x_b.setZero(); + Eigen::VectorXd x(A.rows()); + x.setZero(); + { + amgcl::profiler<> prof("crystm03_Block"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + params["AMGCL"]["block_size"] = 3; + params["AMGCL"]["solver_type"] = "bicgstab"; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x_b); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; + } + { + amgcl::profiler<> prof("crystm03_Scalar"); + json solver_info; + auto solver = LinearSolver::create("AMGCL", ""); + prof.tic("setup"); + json params; + params["AMGCL"]["tolerance"] = 1e-8; + params["AMGCL"]["max_iter"] = 10000; + params["AMGCL"]["solver_type"] = "bicgstab"; + solver->setParameters(params); + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + prof.toc("setup"); + prof.tic("solve"); + solver->solve(b, x); + prof.toc("solve"); + solver->getInfo(solver_info); + REQUIRE(solver_info["num_iterations"] > 0); + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl + << prof << std::endl; + } + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); + REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); } #endif #ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("Hypre_b2", "[solver]") +TEST_CASE("Hyprel_b2", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "gr_30_30.mtx"; @@ -686,6 +600,7 @@ TEST_CASE("Hypre_b2", "[solver]") solver->solve(b, x); end = clock(); solver->getInfo(solver_info); + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } @@ -704,18 +619,21 @@ TEST_CASE("Hypre_b2", "[solver]") solver->solve(b, x_b); end = clock(); solver->getInfo(solver_info); + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } const double err = (A * x - b).norm() / b.norm(); const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; REQUIRE(err < 1e-8); REQUIRE(err_b < 1e-8); } #endif #ifdef POLYSOLVE_WITH_HYPRE -TEST_CASE("Hypre_crystm03", "[solver]") +TEST_CASE("Hybre_crystm03", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "crystm03.mtx"; @@ -742,6 +660,7 @@ TEST_CASE("Hypre_crystm03", "[solver]") solver->solve(b, x); end = clock(); solver->getInfo(solver_info); + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } @@ -760,73 +679,74 @@ TEST_CASE("Hypre_crystm03", "[solver]") solver->solve(b, x_b); end = clock(); solver->getInfo(solver_info); + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; std::cout << solver_info["num_iterations"] << std::endl; std::cout << solver_info["final_res_norm"] << std::endl; } const double err = (A * x - b).norm() / b.norm(); const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; REQUIRE(err < 1e-8); REQUIRE(err_b < 1e-8); } #endif -#ifdef POLYSOLVE_WITH_AMGCL -TEST_CASE("AMGCL_CPU", "[solver]") +#ifdef POLYSOLVE_WITH_HYPRE +TEST_CASE("hypre_smallscale", "[solver]") { -#ifndef NDEBUG - return; -#endif const std::string path = POLYSOLVE_DATA_DIR; - std::string MatrixName = "Serena.mtx"; Eigen::SparseMatrix A; - loadSymmetric(A, path + "/" + MatrixName); - std::cout << "Matrix Load OK" << std::endl; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); Eigen::VectorXd b(A.rows()); b.setOnes(); - Eigen::VectorXd x_b(A.rows()); - x_b.setZero(); - Eigen::VectorXd x(A.rows()); + Eigen::VectorXd x(b.size()); x.setZero(); + Eigen::VectorXd x_b(b.size()); + x_b.setZero(); { - amgcl::profiler<> prof("Serena_Block"); + clock_t start, end; json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 1000; - params["AMGCL"]["block_size"] = 3; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x_b); - prof.toc("solve"); + solver->solve(b, x); + end = clock(); solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout<< prof << std::endl; + std::cout << "Scalar Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; } { - amgcl::profiler<> prof("Serena_Scalar"); + clock_t start, end; json solver_info; - auto solver = LinearSolver::create("AMGCL", ""); - prof.tic("setup"); + start = clock(); + auto solver = LinearSolver::create("Hypre", ""); json params; - params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["max_iter"] = 10000; + params["Hypre"]["block_size"] = 3; + params["Hypre"]["tolerance"] = 1e-8; + params["Hypre"]["max_iter"] = 1000; solver->setParameters(params); solver->analyzePattern(A, A.rows()); solver->factorize(A); - prof.toc("setup"); - prof.tic("solve"); - solver->solve(b, x); - prof.toc("solve"); + solver->solve(b, x_b); + end = clock(); solver->getInfo(solver_info); - REQUIRE(solver_info["num_iterations"] > 0); - std::cout << prof << std::endl; + std::cout << "Block Running time is " << double(end - start) / CLOCKS_PER_SEC << std::endl; + std::cout << solver_info["num_iterations"] << std::endl; + std::cout << solver_info["final_res_norm"] << std::endl; } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); - REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7); + const double err = (A * x - b).norm() / b.norm(); + const double err_b = (A * x_b - b).norm() / b.norm(); + std::cout << "Scalar relative error " << err << std::endl; + std::cout << "Block relative error " << err_b << std::endl; + REQUIRE(err < 1e-8); + REQUIRE(err_b < 1e-8); } -#endif +#endif \ No newline at end of file From 2fdfdddbf315daf98281d7db99429b962246488e Mon Sep 17 00:00:00 2001 From: Yiwei Shao <44545837+njsyw1997@users.noreply.github.com> Date: Thu, 21 Jul 2022 10:41:50 -0400 Subject: [PATCH 16/29] Delete the local nvcc path --- CMakeLists.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 47b27f2..68e864e 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -136,8 +136,6 @@ include(eigen) target_link_libraries(polysolve PUBLIC Eigen3::Eigen) # CUDA -set(CMAKE_CUDA_COMPILER "/usr/local/cuda/bin/nvcc") - if(POLYSOLVE_WITH_CUDA) include(CheckLanguage) check_language(CUDA) From 41adee3a8013e6d7a3c5ea5556c9e5e1d2e11fa0 Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Thu, 21 Jul 2022 11:03:31 -0400 Subject: [PATCH 17/29] Update nvcc path --- CMakeLists.txt | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 68e864e..9a0a6ec 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,18 +137,19 @@ target_link_libraries(polysolve PUBLIC Eigen3::Eigen) # CUDA if(POLYSOLVE_WITH_CUDA) - include(CheckLanguage) - check_language(CUDA) - set(CMAKE_CUDA_ARCHITECTURES "70;75;86") - if(CMAKE_CUDA_COMPILER) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) - enable_language(CUDA) + find_package(CUDAToolkit) + if(CUDAToolkit_FOUND) + if(NOT DEFINED CMAKE_CUDA_COMPILER) + set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE}) + endif() + target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) else() message(FATAL_ERROR "No CUDA support found!") endif() # We need to explicitly state that we need all CUDA files in the particle # library to be built with -dc as the member functions could be called by # other libraries and executables. + set(CMAKE_CUDA_ARCHITECTURES "70;75;86") set_target_properties(polysolve PROPERTIES CUDA_SEPARABLE_COMPILATION ON) # Nvidia RTX8000 -> compute_75 @@ -167,11 +168,11 @@ if(POLYSOLVE_WITH_CUDA) PROPERTY BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) endif() - find_package(CUDAToolkit) polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cudart) polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cusparse) endif() + # Hypre (GNU Lesser General Public License) if(POLYSOLVE_WITH_HYPRE) include(hypre) From 64594096beec01e75ef96900cd2bbe1377287288 Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Tue, 26 Jul 2022 15:36:52 -0400 Subject: [PATCH 18/29] Solve the problem of undefined reference --- CMakeLists.txt | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9a0a6ec..4e6818d 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,11 +137,14 @@ target_link_libraries(polysolve PUBLIC Eigen3::Eigen) # CUDA if(POLYSOLVE_WITH_CUDA) + include(CheckLanguage) + check_language(CUDA) + set(CMAKE_CUDA_ARCHITECTURES "61") find_package(CUDAToolkit) if(CUDAToolkit_FOUND) - if(NOT DEFINED CMAKE_CUDA_COMPILER) - set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE}) - endif() + message(STATUS "Found CUDATooklkit") + set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE}) + enable_language(CUDA) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) else() message(FATAL_ERROR "No CUDA support found!") @@ -149,7 +152,6 @@ if(POLYSOLVE_WITH_CUDA) # We need to explicitly state that we need all CUDA files in the particle # library to be built with -dc as the member functions could be called by # other libraries and executables. - set(CMAKE_CUDA_ARCHITECTURES "70;75;86") set_target_properties(polysolve PROPERTIES CUDA_SEPARABLE_COMPILATION ON) # Nvidia RTX8000 -> compute_75 From 2137b93d3d1a31ac8057609955075921414ae03d Mon Sep 17 00:00:00 2001 From: njsyw1997 Date: Tue, 26 Jul 2022 15:40:36 -0400 Subject: [PATCH 19/29] Call cuda solver in cpp file, try it in polyfem --- cudatest/CMakeLists.txt | 6 ++++-- cudatest/{cudatest.cu => cudatest.cpp} | 1 + 2 files changed, 5 insertions(+), 2 deletions(-) rename cudatest/{cudatest.cu => cudatest.cpp} (99%) diff --git a/cudatest/CMakeLists.txt b/cudatest/CMakeLists.txt index 0ad97d0..f8f97be 100644 --- a/cudatest/CMakeLists.txt +++ b/cudatest/CMakeLists.txt @@ -3,7 +3,7 @@ ################################################################################ set(test_sources main.cpp - cudatest.cu + cudatest.cpp ) add_executable(cuda_test ${test_sources}) @@ -32,4 +32,6 @@ catch_discover_tests(cuda_test) # Data ################################################################################ set(DATA_DIR "${CMAKE_SOURCE_DIR}/tests/data") -target_compile_definitions(cuda_test PUBLIC -DPOLYSOLVE_DATA_DIR=\"${DATA_DIR}\") \ No newline at end of file +target_compile_definitions(cuda_test PUBLIC -DPOLYSOLVE_DATA_DIR=\"${DATA_DIR}\") +polysolve_target_link_system_libraries(cuda_test PRIVATE CUDA::cudart) +polysolve_target_link_system_libraries(cuda_test PRIVATE CUDA::cusparse) \ No newline at end of file diff --git a/cudatest/cudatest.cu b/cudatest/cudatest.cpp similarity index 99% rename from cudatest/cudatest.cu rename to cudatest/cudatest.cpp index e354725..d970177 100644 --- a/cudatest/cudatest.cu +++ b/cudatest/cudatest.cpp @@ -5,6 +5,7 @@ #include #include #include +#include ////////////////////////////////////////////////////////////////////////// #include From 34561712a7d78fc679f10bca5809f27764b5fe0c Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Thu, 27 Oct 2022 17:06:08 -0400 Subject: [PATCH 20/29] Hypre OFF and relaxation as preconditioner by default --- CMakeLists.txt | 113 +++++++++++++----------- src/polysolve/LinearSolverAMGCL_cuda.cu | 25 +----- 2 files changed, 63 insertions(+), 75 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4e6818d..8fa86cb 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,6 @@ # Detects whether this is a top-level project get_directory_property(HAS_PARENT PARENT_DIRECTORY) + if(HAS_PARENT) set(POLYSOLVE_TOPLEVEL_PROJECT OFF) else() @@ -8,6 +9,7 @@ endif() # Check required CMake version set(REQUIRED_CMAKE_VERSION "3.14.0") + if(POLYSOLVE_TOPLEVEL_PROJECT) cmake_minimum_required(VERSION ${REQUIRED_CMAKE_VERSION}) else() @@ -24,34 +26,35 @@ if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/PolySolveOptions.cmake) include(${CMAKE_CURRENT_SOURCE_DIR}/PolySolveOptions.cmake) endif() - -################################################################################ - +# ############################################################################### project(PolySolve - DESCRIPTION "Easy-to-use wrapper for linear solver" - LANGUAGES CXX) + DESCRIPTION "Easy-to-use wrapper for linear solver" + LANGUAGES CXX) # Polysolve options -option(POLYSOLVE_WITH_SANITIZERS "Enable sanitizers in compilation targets" OFF) +option(POLYSOLVE_WITH_SANITIZERS "Enable sanitizers in compilation targets" OFF) + # Polysolve options for enabling/disabling optional libraries -option(POLYSOLVE_WITH_CHOLMOD "Enable Cholmod library" ON) -option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" ON) -option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) -option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) - -option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) -option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) -option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) -option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) -option(POLYSOLVE_WITH_SPECTRA "Enable computing spectrum" ON) +option(POLYSOLVE_WITH_CHOLMOD "Enable Cholmod library" ON) +option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" ON) +option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) +option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) + +option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) +option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) +option(POLYSOLVE_WITH_HYPRE "Enable hypre" OFF) +option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) +option(POLYSOLVE_WITH_SPECTRA "Enable computing spectrum" ON) + # Sanitizer options -option(POLYSOLVE_SANITIZE_ADDRESS "Sanitize Address" OFF) -option(POLYSOLVE_SANITIZE_MEMORY "Sanitize Memory" OFF) -option(POLYSOLVE_SANITIZE_THREAD "Sanitize Thread" OFF) -option(POLYSOLVE_SANITIZE_UNDEFINED "Sanitize Undefined" OFF) +option(POLYSOLVE_SANITIZE_ADDRESS "Sanitize Address" OFF) +option(POLYSOLVE_SANITIZE_MEMORY "Sanitize Memory" OFF) +option(POLYSOLVE_SANITIZE_THREAD "Sanitize Thread" OFF) +option(POLYSOLVE_SANITIZE_UNDEFINED "Sanitize Undefined" OFF) + # Misc. -option(POLYSOLVE_LARGE_INDEX "Build for large indices" OFF) -option(POLYSOLVE_WITH_TESTS "Build unit-tests" ${POLYSOLVE_TOPLEVEL_PROJECT}) +option(POLYSOLVE_LARGE_INDEX "Build for large indices" OFF) +option(POLYSOLVE_WITH_TESTS "Build unit-tests" ${POLYSOLVE_TOPLEVEL_PROJECT}) include(CMakeDependentOption) cmake_dependent_option(SUITE_SPARSE_WITH_MKL "Build SuiteSparse using MKL" ON "POLYSOLVE_WITH_MKL" OFF) @@ -63,16 +66,17 @@ if(POLYSOLVE_TOPLEVEL_PROJECT) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) endif() + if(NOT DEFINED CMAKE_CUDA_STANDARD) set(CMAKE_CUDA_STANDARD 14) set(CMAKE_CUDA_STANDARD_REQUIRED ON) endif() -if (MSVC) +if(MSVC) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /bigobj") endif() -### Configuration +# ## Configuration list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/polysolve/") list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/recipes/") list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/find/") @@ -91,9 +95,9 @@ set_property(GLOBAL PROPERTY USE_FOLDERS ON) # Generate position independent code by default set(CMAKE_POSITION_INDEPENDENT_CODE ON) -################################################################################ +# ############################################################################### # PolySolve Library -################################################################################ +# ############################################################################### # Add an empty library and fill in the list of sources in `src/CMakeLists.txt`. add_library(polysolve) @@ -104,10 +108,9 @@ add_subdirectory(src) # Public include directory for Polysolve target_include_directories(polysolve PUBLIC ${PROJECT_SOURCE_DIR}/src) -################################################################################ +# ############################################################################### # Definitions -################################################################################ - +# ############################################################################### if(POLYSOLVE_LARGE_INDEX) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_LARGE_INDEX) endif() @@ -118,9 +121,9 @@ target_compile_definitions(polysolve PUBLIC -DEIGEN_STACK_ALLOCATION_LIMIT=0) # 8MB # target_compile_definitions(polysolve PUBLIC -DEIGEN_STACK_ALLOCATION_LIMIT=8388608) -################################################################################ +# ############################################################################### # Dependencies -################################################################################ +# ############################################################################### # Extra warnings include(polysolve_warnings) @@ -140,46 +143,50 @@ if(POLYSOLVE_WITH_CUDA) include(CheckLanguage) check_language(CUDA) set(CMAKE_CUDA_ARCHITECTURES "61") - find_package(CUDAToolkit) + find_package(CUDAToolkit) + if(CUDAToolkit_FOUND) message(STATUS "Found CUDATooklkit") set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE}) enable_language(CUDA) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) + target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) else() message(FATAL_ERROR "No CUDA support found!") endif() - # We need to explicitly state that we need all CUDA files in the particle - # library to be built with -dc as the member functions could be called by - # other libraries and executables. - set_target_properties(polysolve PROPERTIES CUDA_SEPARABLE_COMPILATION ON) - - # Nvidia RTX8000 -> compute_75 - # Nvidia V100 -> compute_70 - # Nvidia 1080/1080Ti -> compute_61 - # Nvidia 3080Ti -> compute_86 + + # We need to explicitly state that we need all CUDA files in the particle + # library to be built with -dc as the member functions could be called by + # other libraries and executables. + set_target_properties(polysolve PROPERTIES CUDA_SEPARABLE_COMPILATION ON) + + # Nvidia RTX8000 -> compute_75 + # Nvidia V100 -> compute_70 + # Nvidia 1080/1080Ti -> compute_61 + # Nvidia 3080Ti -> compute_86 if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) set(CMAKE_CUDA_ARCHITECTURES 70 75 86) endif() + set_target_properties(polysolve PROPERTIES CUDA_ARCHITECTURES "70;75;86") if(APPLE) # We need to add the path to the driver (libcuda.dylib) as an rpath, # so that the static cuda runtime can find it at runtime. set_property(TARGET polysolve - PROPERTY - BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) + PROPERTY + BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) endif() + polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cudart) polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cusparse) endif() - # 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) + if(HYPRE_WITH_MPI) target_compile_definitions(polysolve PUBLIC HYPRE_WITH_MPI) endif() @@ -205,6 +212,7 @@ endif() # Pardiso solver if(POLYSOLVE_WITH_PARDISO) include(pardiso) + if(TARGET Pardiso::Pardiso) target_link_libraries(polysolve PUBLIC Pardiso::Pardiso) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PARDISO) @@ -216,6 +224,7 @@ endif() # UmfPack solver if(POLYSOLVE_WITH_UMFPACK) include(umfpack) + if(TARGET UMFPACK::UMFPACK) target_link_libraries(polysolve PUBLIC UMFPACK::UMFPACK) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_UMFPACK) @@ -227,6 +236,7 @@ endif() # SuperLU solver if(POLYSOLVE_WITH_SUPERLU) include(superlu) + if(TARGET SuperLU::SuperLU) target_link_libraries(polysolve PUBLIC SuperLU::SuperLU) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_SUPERLU) @@ -235,7 +245,6 @@ if(POLYSOLVE_WITH_SUPERLU) endif() endif() - # AMGCL solver if(POLYSOLVE_WITH_AMGCL) include(amgcl) @@ -250,16 +259,16 @@ if(POLYSOLVE_WITH_SPECTRA) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_SPECTRA) endif() -################################################################################ +# ############################################################################### # Compiler options -################################################################################ +# ############################################################################### # Use C++14 target_compile_features(polysolve PUBLIC cxx_std_14) -################################################################################ +# ############################################################################### # Tests -################################################################################ +# ############################################################################### # Compile extras only if this is a top-level project if(POLYSOLVE_WITH_TESTS) @@ -273,11 +282,9 @@ if(POLYSOLVE_WITH_TESTS) include("${catch2_SOURCE_DIR}/contrib/Catch.cmake") add_subdirectory(tests) - + # Cuda test if(POLYSOLVE_WITH_CUDA) add_subdirectory(cudatest) endif() endif() - - diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu index 5b09387..40e1374 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.cu +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -34,26 +34,7 @@ namespace polysolve { json params = R"({ "precond": { - "relax": { - "degree": 16, - "type": "chebyshev", - "power_iters": 100, - "higher": 2, - "lower": 0.008333333333, - "scale": true - }, - "class": "amg", - "max_levels": 6, - "direct_coarse": false, - "ncycle": 2, - "coarsening": { - "type": "smoothed_aggregation", - "estimate_spectral_radius": true, - "relax": 1, - "aggr": { - "eps_strong": 0 - } - } + "class": "relaxation" }, "solver": { "tol": 1e-10, @@ -144,7 +125,7 @@ namespace polysolve boost::property_tree::ptree pt_params; boost::property_tree::read_json(ss_params, pt_params); auto A = std::tie(numRows, ia, ja, a); - solver_ = std::make_unique(A, pt_params,backend_params_); + solver_ = std::make_unique(A, pt_params, backend_params_); iterations_ = 0; residual_error_ = 0; } @@ -162,7 +143,7 @@ namespace polysolve auto x_b = Backend::copy_vector(x, backend_params_); std::tie(iterations_, residual_error_) = (*solver_)(*rhs_b, *x_b); - thrust::copy(x_b->begin(), x_b->end(),result.data()); + thrust::copy(x_b->begin(), x_b->end(), result.data()); } //////////////////////////////////////////////////////////////////////////////// From ceb1c68834d4c0400d8a2abf74bb90cf9aeeb343 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Thu, 27 Oct 2022 17:56:18 -0400 Subject: [PATCH 21/29] fixing warnings only for CXX --- cmake/polysolve/polysolve_warnings.cmake | 94 ++++++++++++------------ 1 file changed, 49 insertions(+), 45 deletions(-) diff --git a/cmake/polysolve/polysolve_warnings.cmake b/cmake/polysolve/polysolve_warnings.cmake index a5b6dd3..7c35967 100644 --- a/cmake/polysolve/polysolve_warnings.cmake +++ b/cmake/polysolve/polysolve_warnings.cmake @@ -1,7 +1,7 @@ -################################################################################ +# ############################################################################### # See comments and discussions here: # http://stackoverflow.com/questions/5088460/flags-to-enable-thorough-and-verbose-g-warnings -################################################################################ +# ############################################################################### if(TARGET polysolve::warnings) return() @@ -13,7 +13,7 @@ set(POLYSOLVE_WARNING_FLAGS -pedantic # -Wconversion - #-Wunsafe-loop-optimizations # broken with C++11 loops + # -Wunsafe-loop-optimizations # broken with C++11 loops -Wunused -Wno-long-long @@ -43,11 +43,11 @@ set(POLYSOLVE_WARNING_FLAGS -Wunused-but-set-variable -Wno-unused-parameter - #-Weffc++ + # -Weffc++ -Wno-old-style-cast - # -Wno-sign-conversion - #-Wsign-conversion + # -Wno-sign-conversion + # -Wsign-conversion -Wshadow -Wstrict-null-sentinel @@ -62,8 +62,10 @@ set(POLYSOLVE_WARNING_FLAGS -Wcast-align -Wdisabled-optimization - #-Winline # produces warning on default implicit destructor + + # -Winline # produces warning on default implicit destructor -Winvalid-pch + # -Wmissing-include-dirs -Wpacked -Wno-padded @@ -74,8 +76,8 @@ set(POLYSOLVE_WARNING_FLAGS -Wlogical-op -Wnoexcept -Woverloaded-virtual - # -Wundef + # -Wundef -Wnon-virtual-dtor -Wdelete-non-virtual-dtor -Werror=non-virtual-dtor @@ -83,62 +85,61 @@ set(POLYSOLVE_WARNING_FLAGS -Wno-sign-compare - ########### + # ########## # GCC 6.1 # - ########### - + # ########## -Wnull-dereference -fdelete-null-pointer-checks -Wduplicated-cond -Wmisleading-indentation - #-Weverything + # -Weverything - ########################### + # ########################## # Enabled by -Weverything # - ########################### - - #-Wdocumentation - #-Wdocumentation-unknown-command - #-Wfloat-equal - #-Wcovered-switch-default - - #-Wglobal-constructors - #-Wexit-time-destructors - #-Wmissing-variable-declarations - #-Wextra-semi - #-Wweak-vtables - #-Wno-source-uses-openmp - #-Wdeprecated - #-Wnewline-eof - #-Wmissing-prototypes - - #-Wno-c++98-compat - #-Wno-c++98-compat-pedantic - - ########################### + # ########################## + + # -Wdocumentation + # -Wdocumentation-unknown-command + # -Wfloat-equal + # -Wcovered-switch-default + + # -Wglobal-constructors + # -Wexit-time-destructors + # -Wmissing-variable-declarations + # -Wextra-semi + # -Wweak-vtables + # -Wno-source-uses-openmp + # -Wdeprecated + # -Wnewline-eof + # -Wmissing-prototypes + + # -Wno-c++98-compat + # -Wno-c++98-compat-pedantic + + # ########################## # Need to check if those are still valid today - ########################### + # ########################## - #-Wimplicit-atomic-properties - #-Wmissing-declarations - #-Wmissing-prototypes - #-Wstrict-selector-match - #-Wundeclared-selector - #-Wunreachable-code + # -Wimplicit-atomic-properties + # -Wmissing-declarations + # -Wmissing-prototypes + # -Wstrict-selector-match + # -Wundeclared-selector + # -Wunreachable-code # Not a warning, but enable link-time-optimization # TODO: Check out modern CMake version of setting this flag # https://cmake.org/cmake/help/latest/module/CheckIPOSupported.html - #-flto + # -flto # Gives meaningful stack traces -fno-omit-frame-pointer -fno-optimize-sibling-calls - ##################### + # #################### # Disabled warnings # - ##################### + # #################### -Wno-missing-noreturn -Wno-shadow -Wno-switch-enum @@ -162,10 +163,13 @@ add_library(polysolve::warnings ALIAS polysolve_warnings) foreach(FLAG IN ITEMS ${POLYSOLVE_WARNING_FLAGS}) string(REPLACE "=" "-" FLAG_VAR "${FLAG}") + if(NOT DEFINED IS_SUPPORTED_${FLAG_VAR}) check_cxx_compiler_flag("${FLAG}" IS_SUPPORTED_${FLAG_VAR}) endif() + if(IS_SUPPORTED_${FLAG_VAR}) - target_compile_options(polysolve_warnings INTERFACE ${FLAG}) + target_compile_options(polysolve_warnings INTERFACE $<$: ${FLAG} + >) endif() endforeach() From 3095b349962206f88a02112f490ca39428073ab3 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Fri, 28 Oct 2022 19:28:46 -0400 Subject: [PATCH 22/29] setting cusparse ilu0 and runtime type for relaxation --- src/polysolve/LinearSolverAMGCL_cuda.cu | 3 ++- src/polysolve/LinearSolverAMGCL_cuda.hpp | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu index 40e1374..d9f824d 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.cu +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -34,7 +34,8 @@ namespace polysolve { json params = R"({ "precond": { - "class": "relaxation" + "class": "relaxation", + "type": "ilu0" }, "solver": { "tol": 1e-10, diff --git a/src/polysolve/LinearSolverAMGCL_cuda.hpp b/src/polysolve/LinearSolverAMGCL_cuda.hpp index b3d0dce..ce27c5d 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.hpp +++ b/src/polysolve/LinearSolverAMGCL_cuda.hpp @@ -23,6 +23,10 @@ #include #include #include + +// SET THIS AS AN OPTIONAL HEADER +#include + #include #include #include From 803de42b475165177ebfb5bfa59d7e08ddc7a61d Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Fri, 28 Oct 2022 20:35:32 -0400 Subject: [PATCH 23/29] fixing AMGCL_cuda interface, Cusparseilu0 optional --- src/polysolve/LinearSolverAMGCL_cuda.cu | 12 ++++++------ src/polysolve/LinearSolverAMGCL_cuda.hpp | 2 ++ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu index d9f824d..ec8b4bc 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.cu +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -49,13 +49,13 @@ namespace polysolve void set_params(const json ¶ms, json &out) { - if (params.contains("AMGCL")) + if (params.contains("AMGCL_cuda")) { // Patch the stored params with input ones - if (params["AMGCL"].contains("precond")) - out["precond"].merge_patch(params["AMGCL"]["precond"]); - if (params["AMGCL"].contains("solver")) - out["solver"].merge_patch(params["AMGCL"]["solver"]); + if (params["AMGCL_cuda"].contains("precond")) + out["precond"].merge_patch(params["AMGCL_cuda"]["precond"]); + if (params["AMGCL_cuda"].contains("solver")) + out["solver"].merge_patch(params["AMGCL_cuda"]["solver"]); if (out["precond"]["class"] == "schur_pressure_correction") { @@ -89,7 +89,7 @@ namespace polysolve // Set solver parameters void LinearSolverAMGCL_cuda::setParameters(const json ¶ms) { - if (params.contains("AMGCL")) + if (params.contains("AMGCL_cuda")) { set_params(params, params_); } diff --git a/src/polysolve/LinearSolverAMGCL_cuda.hpp b/src/polysolve/LinearSolverAMGCL_cuda.hpp index ce27c5d..a71bf78 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.hpp +++ b/src/polysolve/LinearSolverAMGCL_cuda.hpp @@ -25,7 +25,9 @@ #include // SET THIS AS AN OPTIONAL HEADER +#ifdef POLYSOLVE_WITH_CUSPARSEILU0 #include +#endif #include #include From 288b9bda62af8176c9fdbc926021a10a696d271d Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Fri, 11 Nov 2022 10:50:49 -0500 Subject: [PATCH 24/29] AMG preconditioner as default for AMGCL_CUDA, Setting CUDA and HYPRE off as default --- CMakeLists.txt | 8 ++++---- src/polysolve/LinearSolver.cpp | 6 +++++- src/polysolve/LinearSolverAMGCL_cuda.cu | 18 ++++++++++++++++-- 3 files changed, 25 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8629535..56215d1 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,8 +47,8 @@ option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" ON) option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) -option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) -option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) +option(POLYSOLVE_WITH_CUDA "Enable CUDA" OFF) +option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" ON) option(POLYSOLVE_WITH_HYPRE "Enable hypre" OFF) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) option(POLYSOLVE_WITH_SPECTRA "Enable computing spectrum" ON) @@ -171,10 +171,10 @@ if(POLYSOLVE_WITH_CUDA) # Nvidia 1080/1080Ti -> compute_61 # Nvidia 3080Ti -> compute_86 if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) - set(CMAKE_CUDA_ARCHITECTURES 70 75 86) + set(CMAKE_CUDA_ARCHITECTURES 60 70 75 86) endif() - set_target_properties(polysolve PROPERTIES CUDA_ARCHITECTURES "70;75;86") + set_target_properties(polysolve PROPERTIES CUDA_ARCHITECTURES "60;70;75;86") if(APPLE) # We need to add the path to the driver (libcuda.dylib) as an rpath, diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp index b4bd18b..d10348f 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/LinearSolver.cpp @@ -228,10 +228,12 @@ namespace polysolve else if (solver == "AMGCL") { return std::make_unique(); +#ifdef POLYSOLVE_WITH_CUDA } else if (solver == "AMGCL_cuda") { - return std::make_unique(); + return std::make_unique(); +#endif #endif #if EIGEN_VERSION_AT_LEAST(3, 3, 0) @@ -302,8 +304,10 @@ namespace polysolve #endif #ifdef POLYSOLVE_WITH_AMGCL "AMGCL", +#ifdef POLYSOLVE_WITH_CUDA "AMGCL_cuda", #endif +#endif #if EIGEN_VERSION_AT_LEAST(3, 3, 0) #ifndef POLYSOLVE_LARGE_INDEX "Eigen::LeastSquaresConjugateGradient", diff --git a/src/polysolve/LinearSolverAMGCL_cuda.cu b/src/polysolve/LinearSolverAMGCL_cuda.cu index ec8b4bc..5a5246d 100644 --- a/src/polysolve/LinearSolverAMGCL_cuda.cu +++ b/src/polysolve/LinearSolverAMGCL_cuda.cu @@ -34,8 +34,21 @@ namespace polysolve { json params = R"({ "precond": { - "class": "relaxation", - "type": "ilu0" + "relax": { + "type": "spai0" + }, + "class": "amg", + "max_levels": 6, + "direct_coarse": false, + "ncycle": 2, + "coarsening": { + "type": "smoothed_aggregation", + "estimate_spectral_radius": true, + "relax": 1, + "aggr": { + "eps_strong": 0 + } + } }, "solver": { "tol": 1e-10, @@ -127,6 +140,7 @@ namespace polysolve boost::property_tree::read_json(ss_params, pt_params); auto A = std::tie(numRows, ia, ja, a); solver_ = std::make_unique(A, pt_params, backend_params_); + // std::cout << *solver_.get() << std::endl; iterations_ = 0; residual_error_ = 0; } From 1351f5829cf25dba28e73d8026c125d30f4f3b58 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Sat, 10 Dec 2022 19:01:37 -0500 Subject: [PATCH 25/29] added PETSC MUMPS solver --- .gitignore | 1 + CMakeLists.txt | 32 +++++-- cudatest/cudatest.cpp | 141 ++++++++++++++++++++++++++-- src/CMakeLists.txt | 3 +- src/polysolve/LinearSolver.cpp | 12 +++ src/polysolve/LinearSolver.hpp | 13 ++- src/polysolve/LinearSolverPETSC.cpp | 126 +++++++++++++++++++++++++ src/polysolve/LinearSolverPETSC.hpp | 74 +++++++++++++++ 8 files changed, 384 insertions(+), 18 deletions(-) create mode 100644 src/polysolve/LinearSolverPETSC.cpp create mode 100644 src/polysolve/LinearSolverPETSC.hpp diff --git a/.gitignore b/.gitignore index f4b3db8..3f2ca06 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ build* +*build* bin* 3rdparty* tests/data diff --git a/CMakeLists.txt b/CMakeLists.txt index 19c186d..e31ce4a 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,17 +42,18 @@ project(PolySolve option(POLYSOLVE_WITH_SANITIZERS "Enable sanitizers in compilation targets" OFF) # Polysolve options for enabling/disabling optional libraries -option(POLYSOLVE_WITH_CUDA "Enable CUDA" OFF) - option(POLYSOLVE_WITH_CHOLMOD "Enable Cholmod library" ON) option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" ON) option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) -option(POLYSOLVE_WITH_CUSOLVER "Enable cuSOLVER library" OFF) -option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) +option(POLYSOLVE_WITH_CUSOLVER "Enable cuSOLVER library" ON) + +option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) +option(POLYSOLVE_WITH_PETSC "Enable PETSC" ON) +option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" ON) option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) -option(HYPRE_WITH_MPI "Enable hypre" OFF) +option(HYPRE_WITH_MPI "Enable hypre MPI" OFF) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) option(POLYSOLVE_WITH_SPECTRA "Enable computing spectrum" ON) @@ -156,7 +157,7 @@ if(POLYSOLVE_WITH_CUDA) find_package(CUDAToolkit) if(CUDAToolkit_FOUND) - message(STATUS "Found CUDATooklkit") + message(STATUS "Found CUDAToolkit") set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE}) enable_language(CUDA) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUDA) @@ -199,6 +200,7 @@ if(POLYSOLVE_WITH_HYPRE) if(HYPRE_WITH_MPI) target_compile_definitions(polysolve PUBLIC HYPRE_WITH_MPI) + message(STATUS "HYPRE WITH MPI.") endif() endif() @@ -278,6 +280,7 @@ if(POLYSOLVE_WITH_CUSOLVER) if(TARGET CUDA::cusolver) target_link_libraries(polysolve PUBLIC CUDA::cusolver) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_CUSOLVER) + message(STATUS "cuSOLVER found!") else() message(WARNING "cuSOLVER not found, solver will not be available.") endif() @@ -285,6 +288,23 @@ endif() # ############################################################################### +# PETSC solvers +if(POLYSOLVE_WITH_PETSC) + # include(PETSC PETSc) + find_package(MPI) + + # if(TARGET PETSC::PETSc) + target_link_libraries(polysolve PUBLIC -lpetsc MPI::MPI_CXX) + target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PETSC) + message(STATUS "PETSc found!") + + # else() + # message(WARNING "PETSc not found, solver will not be available.") + # endif() +endif() + +# ############################################################################### + # Compiler options # ############################################################################### diff --git a/cudatest/cudatest.cpp b/cudatest/cudatest.cpp index d970177..8511133 100644 --- a/cudatest/cudatest.cpp +++ b/cudatest/cudatest.cpp @@ -38,7 +38,8 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) A.setFromTriplets(triple.begin(), triple.end()); }; -TEST_CASE("amgcl_crystm03_cg", "[solver]"){ +TEST_CASE("amgcl_crystm03_cg", "[solver]") +{ const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "crystm03.mtx"; Eigen::SparseMatrix A; @@ -68,12 +69,13 @@ TEST_CASE("amgcl_crystm03_cg", "[solver]"){ prof.toc("solve"); solver->getInfo(solver_info); REQUIRE(solver_info["num_iterations"] > 0); - std::cout<< prof << std::endl; + std::cout << prof << std::endl; } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); } -TEST_CASE("amgcl_crystm03_bicgstab", "[solver]"){ +TEST_CASE("amgcl_crystm03_bicgstab", "[solver]") +{ const std::string path = POLYSOLVE_DATA_DIR; std::string MatrixName = "crystm03.mtx"; Eigen::SparseMatrix A; @@ -103,7 +105,134 @@ TEST_CASE("amgcl_crystm03_bicgstab", "[solver]"){ prof.toc("solve"); solver->getInfo(solver_info); REQUIRE(solver_info["num_iterations"] > 0); - std::cout<< prof << std::endl; + std::cout << prof << std::endl; + } + REQUIRE((A * x - b).norm() / b.norm() < 1e-7); +} + +TEST_CASE("PETSC_MUMPS", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solver = LinearSolver::create("PETSC_Solver", ""); + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A, 0); + solver->solve(b, x); + + // std::cout<<"Solver error: "< A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solver = LinearSolver::create("cuSolverDN", ""); + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // std::cout<<"Solver error: "<setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // std::cout<<"Solver error: "<> A(m, n); + } + } + + Eigen::VectorXd b(A.rows()); + std::string gradient_path = path + "/matrixdata-5cubes/gradient" + std::to_string(i) + ".txt"; + std::ifstream gradient_file(gradient_path); + for (int m = 0; m < 120; m++) + { + gradient_file >> b(m); + } + + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + + // std::chrono::steady_clock::time_point beginf = std::chrono::steady_clock::now(); + solver->factorize(A); + // std::chrono::steady_clock::time_point endf = std::chrono::steady_clock::now(); + // std::cout << "time to factorize: " << std::chrono::duration_cast(endf-beginf).count() << std::endl; + // factorize_times_file << std::chrono::duration_cast(endf-beginf).count() << " "; + + // std::chrono::steady_clock::time_point begins = std::chrono::steady_clock::now(); + solver->solve(b, x); + // std::chrono::steady_clock::time_point ends = std::chrono::steady_clock::now(); + // std::cout << "time to solve: " << std::chrono::duration_cast(ends-begins).count() << std::endl; + // solve_times_file << std::chrono::duration_cast(ends-begins).count() << " "; + + // std::cout << "Ax norm: " << (A*x).norm() << std::endl; + // std::cout << "b norm: " << b.norm() << std::endl; + + const double err = (A * x - b).norm(); + REQUIRE(err < 1e-8); } - REQUIRE((A * x - b).norm() / b.norm() < 1e-7); } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 67133ea..9e8fec9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,7 +7,9 @@ set(SOURCES polysolve/LinearSolverAMGCL.cpp polysolve/LinearSolverAMGCL.hpp polysolve/LinearSolverCuSolverDN.cu + polysolve/LinearSolverPETSC.cpp polysolve/LinearSolverCuSolverDN.cuh + polysolve/LinearSolverPETSC.hpp polysolve/LinearSolverEigen.hpp polysolve/LinearSolverEigen.tpp polysolve/LinearSolverHypre.cpp @@ -23,5 +25,4 @@ set(SOURCES polysolve_prepend_current_path(SOURCES) polysolve_set_source_group(${SOURCES}) - target_sources(polysolve PRIVATE ${SOURCES}) diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp index 3f71ce7..3822522 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/LinearSolver.cpp @@ -32,6 +32,9 @@ #ifdef POLYSOLVE_WITH_CUSOLVER #include #endif +#ifdef POLYSOLVE_WITH_PETSC +#include +#endif #include //////////////////////////////////////////////////////////////////////////////// @@ -230,6 +233,12 @@ namespace polysolve { return std::make_unique(); #endif +#ifdef POLYSOLVE_WITH_PETSC + } + else if (solver == "PETSC_Solver") + { + return std::make_unique(); +#endif #ifdef POLYSOLVE_WITH_HYPRE } else if (solver == "Hypre") @@ -316,6 +325,9 @@ namespace polysolve #ifdef POLYSOLVE_WITH_CUSOLVER "cuSolverDN", #endif +#ifdef POLYSOLVE_WITH_PETSC + "PETSC_Solver", +#endif #ifdef POLYSOLVE_WITH_HYPRE "Hypre", #endif diff --git a/src/polysolve/LinearSolver.hpp b/src/polysolve/LinearSolver.hpp index 9654542..77b4d36 100644 --- a/src/polysolve/LinearSolver.hpp +++ b/src/polysolve/LinearSolver.hpp @@ -9,9 +9,9 @@ using json = nlohmann::json; #include #define POLYSOLVE_DELETE_MOVE_COPY(Base) \ - Base(Base &&) = delete; \ - Base &operator=(Base &&) = delete; \ - Base(const Base &) = delete; \ + Base(Base &&) = delete; \ + Base &operator=(Base &&) = delete; \ + Base(const Base &) = delete; \ Base &operator=(const Base &) = delete; //////////////////////////////////////////////////////////////////////////////// @@ -30,8 +30,8 @@ namespace polysolve typedef Eigen::SparseMatrix StiffnessMatrix; #endif /** - * @brief Base class for linear solver. - */ + * @brief Base class for linear solver. + */ class LinearSolver { @@ -85,6 +85,9 @@ namespace polysolve // Factorize system matrix virtual void factorize(const StiffnessMatrix &A) {} + // Factorize system matrix + virtual void factorize(StiffnessMatrix &A, int dummy) {} + // Analyze sparsity pattern of a dense matrix virtual void analyzePattern(const Eigen::MatrixXd &A, const int precond_num) {} diff --git a/src/polysolve/LinearSolverPETSC.cpp b/src/polysolve/LinearSolverPETSC.cpp new file mode 100644 index 0000000..6ccc22c --- /dev/null +++ b/src/polysolve/LinearSolverPETSC.cpp @@ -0,0 +1,126 @@ +#ifdef POLYSOLVE_WITH_PETSC + +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include + +namespace polysolve +{ + LinearSolverPETSC::LinearSolverPETSC() + { + init(); + } + + void LinearSolverPETSC::init() + { + // PetscCall(PetscInitialize(NULL, NULL, NULL, NULL)); + // PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); + PetscInitialize(NULL, NULL, NULL, NULL); + PetscDeviceInitialize(PETSC_DEVICE_CUDA); + return; + } + + void LinearSolverPETSC::setParameters(const json ¶ms) + { + } + + //////////////////////////////////////////////////////////////////////////////// + + void LinearSolverPETSC::getInfo(json ¶ms) const + { + } + + void LinearSolverPETSC::analyzePattern(const StiffnessMatrix &A, const int precond_num) + { + } + + void LinearSolverPETSC::analyzePattern(const Eigen::MatrixXd &A, const int precond_num) + { + } + + void LinearSolverPETSC::factorize(StiffnessMatrix &A, int dummy) + { + // PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(), &A_petsc)); + // PetscCall(MatConvert(A_petsc, MATAIJCUSPARSE, MAT_INPLACE_MATRIX, &A_petsc)); + // PetscCall(MatTranspose(A_petsc, MAT_INPLACE_MATRIX, &A_petsc)); + // PetscCall(MatCreateVecs(A_petsc, &x_petsc, NULL)); + MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(), &A_petsc); + // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A_outer, A_inner, A_value, &A_petsc); + MatConvert(A_petsc, MATAIJCUSPARSE, MAT_INPLACE_MATRIX, &A_petsc); + MatTranspose(A_petsc, MAT_INPLACE_MATRIX, &A_petsc); + MatCreateVecs(A_petsc, &x_petsc, NULL); + } + + void LinearSolverPETSC::solve(const Ref b, Ref x) + { + // copy b to device + // PetscCall(VecCreateSeqCUDAWithArrays(PETSC_COMM_SELF, 1, b.rows() * b.cols(), b.data(), NULL, &b_petsc)); + // + // PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); + // PetscCall(KSPSetOperators(ksp, A_petsc, A_petsc)); + // PetscCall(KSPSetTolerances(ksp, PETSC_DEFAULT, 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT)); + // + // KSPSetType(ksp, KSPPREONLY); + // PetscCall(KSPGetPC(ksp, &pc)); + // PCSetType(pc, PCLU); + // PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); + // + // PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */ + // PCFactorGetMatrix(pc, &F); + // MatMumpsSetCntl(F, 3, 1.e-6); + // + // PetscCall(KSPSetFromOptions(ksp)); + // KSPSetUp(ksp); + // PetscCall(KSPSolve(ksp, b_petsc, x_petsc)); + VecCreateSeqCUDAWithArrays(PETSC_COMM_SELF, 1, b.rows() * b.cols(), b.data(), NULL, &b_petsc); + + KSPCreate(PETSC_COMM_SELF, &ksp); + KSPSetOperators(ksp, A_petsc, A_petsc); + KSPSetTolerances(ksp, PETSC_DEFAULT, 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT); + // + KSPSetType(ksp, KSPPREONLY); + KSPGetPC(ksp, &pc); + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); + // + PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */ + PCFactorGetMatrix(pc, &F); + MatMumpsSetCntl(F, 3, 1.e-6); + // + KSPSetFromOptions(ksp); + KSPSetUp(ksp); + KSPSolve(ksp, b_petsc, x_petsc); + + x.resize(b.rows() * b.cols(), 1); + int pidx = 0; + for (int i = 0; i < b.rows() * b.cols(); ++i) + { + PetscInt ix[] = {pidx}; + PetscScalar y[] = {0}; + VecGetValues(x_petsc, 1, ix, y); + x(i, 0) = y[0]; + pidx++; + } + } + + //////////////////////////////////////////////////////////////////////////////// + + LinearSolverPETSC::~LinearSolverPETSC() + { + // PetscCall(KSPDestroy(&ksp)); + // PetscCall(MatDestroy(&A_petsc)); + // PetscCall(VecDestroy(&b_petsc)); + // PetscCall(VecDestroy(&x_petsc)); + // PetscCall(PetscFinalize()); + KSPDestroy(&ksp); + MatDestroy(&A_petsc); + VecDestroy(&b_petsc); + VecDestroy(&x_petsc); + PetscFinalize(); + } + +} // namespace polysolve + +#endif \ No newline at end of file diff --git a/src/polysolve/LinearSolverPETSC.hpp b/src/polysolve/LinearSolverPETSC.hpp new file mode 100644 index 0000000..5af76bf --- /dev/null +++ b/src/polysolve/LinearSolverPETSC.hpp @@ -0,0 +1,74 @@ +#pragma once + +#ifdef POLYSOLVE_WITH_PETSC + +//////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////// +// +// https://docs.nvidia.com/cuda/cusolver/index.html#cuSolverDN-function-reference +// + +namespace polysolve +{ + class LinearSolverPETSC : public LinearSolver + { + + public: + LinearSolverPETSC(); + ~LinearSolverPETSC(); + + private: + POLYSOLVE_DELETE_MOVE_COPY(LinearSolverPETSC) + + public: + ////////////////////// + // Public interface // + ////////////////////// + + // Set solver parameters + virtual void setParameters(const json ¶ms) override; + + // Retrieve memory information from cuSolverDN + virtual void getInfo(json ¶ms) const override; + + // Analyze sparsity pattern (sparse) + virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override; + + // Factorize system matrix (sparse) + virtual void factorize(StiffnessMatrix &A, int dummy) override; + + // Analyze sparsity pattern (dense, preferred) + virtual void analyzePattern(const Eigen::MatrixXd &A, const int precond_num) override; + + // Solve the linear system Ax = b + virtual void solve(const Ref b, Ref x) override; + + // Name of the solver type (for debugging purposes) + virtual std::string name() const override { return "PETSC_Solver"; } + + protected: + void init(); + + protected: + // PETSC variables + Vec b_petsc, x_petsc, y_petsc; + Mat A_petsc, F; + KSP ksp; + PC pc; + PetscReal norm; + PetscInt its; + + // Eigen variables + Eigen::SparseMatrix A; + }; + +} // namespace polysolve + +#endif From a8beae0ca6bc9dc184423216465d1d18eb30788b Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Wed, 14 Dec 2022 17:00:05 -0500 Subject: [PATCH 26/29] Adding multiple external sparse solvers with PETSC backend, setting GPU error check suitable for polyfem --- CMakeLists.txt | 3 +- cudatest/cudatest.cpp | 48 +++---- src/polysolve/LinearSolver.hpp | 4 +- src/polysolve/LinearSolverCuSolverDN.cu | 54 ++++---- src/polysolve/LinearSolverPETSC.cpp | 173 ++++++++++++++++-------- src/polysolve/LinearSolverPETSC.hpp | 17 ++- 6 files changed, 190 insertions(+), 109 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e31ce4a..c8a448e 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,7 +52,7 @@ option(POLYSOLVE_WITH_CUSOLVER "Enable cuSOLVER library" ON) option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) option(POLYSOLVE_WITH_PETSC "Enable PETSC" ON) option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" ON) -option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) +option(POLYSOLVE_WITH_HYPRE "Enable hypre" OFF) option(HYPRE_WITH_MPI "Enable hypre MPI" OFF) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) option(POLYSOLVE_WITH_SPECTRA "Enable computing spectrum" ON) @@ -293,6 +293,7 @@ if(POLYSOLVE_WITH_PETSC) # include(PETSC PETSc) find_package(MPI) + # TODO: CHECK IF PETSC IS AVAILABLE # if(TARGET PETSC::PETSc) target_link_libraries(polysolve PUBLIC -lpetsc MPI::MPI_CXX) target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PETSC) diff --git a/cudatest/cudatest.cpp b/cudatest/cudatest.cpp index 8511133..c364792 100644 --- a/cudatest/cudatest.cpp +++ b/cudatest/cudatest.cpp @@ -38,6 +38,29 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) A.setFromTriplets(triple.begin(), triple.end()); }; +TEST_CASE("PETSC-STRUMPACK", "[solver]") +{ + const std::string path = POLYSOLVE_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/crystm03.mtx"); + REQUIRE(ok); + + auto solver = LinearSolver::create("PETSC_Solver", ""); + // solver->setParameters(params); + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A, 0, 5); + solver->solve(b, x); + + // std::cout<<"Solver error: "< A; - const bool ok = loadMarket(A, path + "/A_2.mat"); - REQUIRE(ok); - - auto solver = LinearSolver::create("PETSC_Solver", ""); - // solver->setParameters(params); - Eigen::VectorXd b(A.rows()); - b.setRandom(); - Eigen::VectorXd x(b.size()); - x.setZero(); - - solver->analyzePattern(A, A.rows()); - solver->factorize(A, 0); - solver->solve(b, x); - - // std::cout<<"Solver error: "< #include -inline void gpuErrchk(cudaError_t code) { - if (code != cudaSuccess) { - throw cudaGetErrorString(code); - } +#define gpuErrchk(ans) \ + { \ + gpuAssert((ans), __FILE__, __LINE__); \ + } +inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) +{ + if (code != cudaSuccess) + { + fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); + if (abort) + exit(code); + } } //////////////////////////////////////////////////////////////////////////////// @@ -30,24 +38,20 @@ namespace polysolve void LinearSolverCuSolverDN::setParameters(const json ¶ms) { - } //////////////////////////////////////////////////////////////////////////////// void LinearSolverCuSolverDN::getInfo(json ¶ms) const { - } void LinearSolverCuSolverDN::analyzePattern(const StiffnessMatrix &A, const int precond_num) { - } void LinearSolverCuSolverDN::analyzePattern(const Eigen::MatrixXd &A, const int precond_num) { - } void LinearSolverCuSolverDN::factorize(const StiffnessMatrix &A) @@ -61,10 +65,10 @@ namespace polysolve { numrows = (int)A.rows(); - //copy A to device + // copy A to device gpuErrchk(cudaMalloc(reinterpret_cast(&d_A), sizeof(double) * A.size())); gpuErrchk(cudaMemcpy(d_A, (const void *)A.data(), sizeof(double) * A.size(), cudaMemcpyHostToDevice)); - + cusolverDnXgetrf_bufferSize(cuHandle, cuParams, numrows, numrows, CUDA_R_64F, d_A, numrows, CUDA_R_64F, &d_lwork, &h_lwork); gpuErrchk(cudaMalloc(reinterpret_cast(&d_work), sizeof(double) * d_lwork)); @@ -72,12 +76,13 @@ namespace polysolve gpuErrchk(cudaMalloc(reinterpret_cast(&d_info), sizeof(int))); gpuErrchk(cudaMalloc(reinterpret_cast(&d_Ipiv), sizeof(int64_t) * numrows)); int info = 0; - - //factorize - cusolverStatus_t solvererr = cusolverDnXgetrf(cuHandle, cuParams, numrows, numrows, CUDA_R_64F, d_A, - numrows, d_Ipiv, CUDA_R_64F, d_work, d_lwork, h_work, h_lwork, d_info); - if(solvererr == CUSOLVER_STATUS_INVALID_VALUE){ + // factorize + cusolverStatus_t solvererr = cusolverDnXgetrf(cuHandle, cuParams, numrows, numrows, CUDA_R_64F, d_A, + numrows, d_Ipiv, CUDA_R_64F, d_work, d_lwork, h_work, h_lwork, d_info); + + if (solvererr == CUSOLVER_STATUS_INVALID_VALUE) + { throw std::invalid_argument("CUDA returned invalid value"); } @@ -86,21 +91,22 @@ namespace polysolve void LinearSolverCuSolverDN::solve(const Ref b, Ref x) { - //copy b to device + // copy b to device gpuErrchk(cudaMalloc(reinterpret_cast(&d_b), sizeof(double) * b.size())); gpuErrchk(cudaMemcpy(d_b, (const void *)b.data(), sizeof(double) * b.size(), cudaMemcpyHostToDevice)); - //solve + // solve cusolverStatus_t solvererr = cusolverDnXgetrs(cuHandle, cuParams, CUBLAS_OP_N, numrows, 1, - CUDA_R_64F, d_A, numrows, d_Ipiv, - CUDA_R_64F, d_b, numrows, d_info); - if(solvererr == CUSOLVER_STATUS_INVALID_VALUE){ + CUDA_R_64F, d_A, numrows, d_Ipiv, + CUDA_R_64F, d_b, numrows, d_info); + if (solvererr == CUSOLVER_STATUS_INVALID_VALUE) + { throw std::invalid_argument("CUDA returned invalid value"); } int info = 0; gpuErrchk(cudaMemcpyAsync(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost, stream)); - //copy result to x + // copy result to x gpuErrchk(cudaMemcpy(x.data(), d_b, sizeof(double) * x.size(), cudaMemcpyDeviceToHost)); } @@ -117,9 +123,9 @@ namespace polysolve cusolverDnDestroyParams(cuParams); cusolverDnDestroy(cuHandle); cudaStreamDestroy(stream); - cudaDeviceReset(); + // cudaDeviceReset(); //not suitable for polyfem } - -}//namespace polysolve + +} // namespace polysolve #endif \ No newline at end of file diff --git a/src/polysolve/LinearSolverPETSC.cpp b/src/polysolve/LinearSolverPETSC.cpp index 6ccc22c..8a8b8c8 100644 --- a/src/polysolve/LinearSolverPETSC.cpp +++ b/src/polysolve/LinearSolverPETSC.cpp @@ -13,13 +13,11 @@ namespace polysolve init(); } - void LinearSolverPETSC::init() + int LinearSolverPETSC::init() { - // PetscCall(PetscInitialize(NULL, NULL, NULL, NULL)); - // PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); - PetscInitialize(NULL, NULL, NULL, NULL); - PetscDeviceInitialize(PETSC_DEVICE_CUDA); - return; + PetscCall(PetscInitialize(NULL, NULL, NULL, NULL)); + PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); + return 0; } void LinearSolverPETSC::setParameters(const json ¶ms) @@ -40,60 +38,131 @@ namespace polysolve { } - void LinearSolverPETSC::factorize(StiffnessMatrix &A, int dummy) + void LinearSolverPETSC::factorize(StiffnessMatrix &A, int AIJ_CUSPARSE, int SOLVER_INDEX) { - // PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(), &A_petsc)); - // PetscCall(MatConvert(A_petsc, MATAIJCUSPARSE, MAT_INPLACE_MATRIX, &A_petsc)); - // PetscCall(MatTranspose(A_petsc, MAT_INPLACE_MATRIX, &A_petsc)); - // PetscCall(MatCreateVecs(A_petsc, &x_petsc, NULL)); - MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(), &A_petsc); - // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.rows(), A.cols(), A_outer, A_inner, A_value, &A_petsc); + GPU_vec = AIJ_CUSPARSE; + /*CHOLMOD requires 64-bit int indices for GPU backend support*/ +#ifdef CHOLMOD_WITH_GPU + std::vector outer_(A.outerSize() + 1); + std::vector inner_(A.nonZeros()); + for (int k = 0; k < A.outerSize() + 1; ++k) + { + outer_[k] = A.outerIndexPtr()[k]; + } + for (int j = 0; j < A.nonZeros(); ++j) + { + inner_[j] = A.innerIndexPtr()[j]; + } + MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, A.rows(), A.cols(), outer_.data(), inner_.data(), A.valuePtr(), &A_petsc); MatConvert(A_petsc, MATAIJCUSPARSE, MAT_INPLACE_MATRIX, &A_petsc); +#else + MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(), &A_petsc); + if (AIJ_CUSPARSE) + MatConvert(A_petsc, MATAIJCUSPARSE, MAT_INPLACE_MATRIX, &A_petsc); +#endif + // AS EIGEN IS COLUMN MAJOR WE DO A TRANSPOSE MatTranspose(A_petsc, MAT_INPLACE_MATRIX, &A_petsc); + MatCreateVecs(A_petsc, &x_petsc, NULL); - } - void LinearSolverPETSC::solve(const Ref b, Ref x) - { - // copy b to device - // PetscCall(VecCreateSeqCUDAWithArrays(PETSC_COMM_SELF, 1, b.rows() * b.cols(), b.data(), NULL, &b_petsc)); - // - // PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); - // PetscCall(KSPSetOperators(ksp, A_petsc, A_petsc)); - // PetscCall(KSPSetTolerances(ksp, PETSC_DEFAULT, 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT)); - // - // KSPSetType(ksp, KSPPREONLY); - // PetscCall(KSPGetPC(ksp, &pc)); - // PCSetType(pc, PCLU); - // PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); - // - // PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */ - // PCFactorGetMatrix(pc, &F); - // MatMumpsSetCntl(F, 3, 1.e-6); - // - // PetscCall(KSPSetFromOptions(ksp)); - // KSPSetUp(ksp); - // PetscCall(KSPSolve(ksp, b_petsc, x_petsc)); - VecCreateSeqCUDAWithArrays(PETSC_COMM_SELF, 1, b.rows() * b.cols(), b.data(), NULL, &b_petsc); - - KSPCreate(PETSC_COMM_SELF, &ksp); + KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, A_petsc, A_petsc); - KSPSetTolerances(ksp, PETSC_DEFAULT, 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT); - // - KSPSetType(ksp, KSPPREONLY); + KSPSetTolerances(ksp, PETSC_DEFAULT, 1e-50, PETSC_DEFAULT, PETSC_DEFAULT); + if (SOLVER_INDEX == 5) + KSPSetType(ksp, KSPGMRES); + else + KSPSetType(ksp, KSPPREONLY); KSPGetPC(ksp, &pc); - PCSetType(pc, PCLU); - PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); - // + + switch (SOLVER_INDEX) + { + case 0: + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO); + break; + case 1: + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU_DIST); + break; + case 2: + PCSetType(pc, PCCHOLESKY); + PCFactorSetMatSolverType(pc, MATSOLVERCHOLMOD); + break; + case 3: + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); + break; + case 4: + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERCUSPARSE); + break; + case 5: + PCSetType(pc, PCLU); + PCFactorSetMatSolverType(pc, MATSOLVERSTRUMPACK); + break; + case 6: + // TODO + PCSetType(pc, PCHYPRE); + PCHYPRESetType(pc, "boomeramg"); + // PCFactorSetMatSolverType(pc, MATSOLVERCUSPARSE); + break; + default: + // TODO + break; + } + PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */ PCFactorGetMatrix(pc, &F); - MatMumpsSetCntl(F, 3, 1.e-6); - // - KSPSetFromOptions(ksp); - KSPSetUp(ksp); + + if (SOLVER_INDEX == 3) + { + MatMumpsSetIcntl(F, 7, 2); + /* threshold for row pivot detection */ + MatMumpsSetIcntl(F, 24, 1); + MatMumpsSetCntl(F, 3, 1.e-6); + } + if (SOLVER_INDEX == 1) + { + MatSuperluSetILUDropTol(F, 1.e-8); + } + if (SOLVER_INDEX == 5) + { + /* Set the fill-reducing reordering. */ + MatSTRUMPACKSetReordering(F, MAT_STRUMPACK_METIS); + /* Since this is a simple discretization, the diagonal is always */ + /* nonzero, and there is no need for the extra MC64 permutation. */ + MatSTRUMPACKSetColPerm(F, PETSC_FALSE); + /* The compression tolerance used when doing low-rank compression */ + /* in the preconditioner. This is problem specific! */ + MatSTRUMPACKSetHSSRelTol(F, 1.e-3); + /* Set minimum matrix size for HSS compression to 15 in order to */ + /* demonstrate preconditioner on small problems. For performance */ + /* a value of say 500 is better. */ + MatSTRUMPACKSetHSSMinSepSize(F, 500); + /* You can further limit the fill in the preconditioner by */ + /* setting a maximum rank */ + MatSTRUMPACKSetHSSMaxRank(F, 100); + /* Set the size of the diagonal blocks (the leafs) in the HSS */ + /* approximation. The default value should be better for real */ + /* problems. This is mostly for illustration on a small problem. */ + // MatSTRUMPACKSetHSSLeafSize(F, 4); + } + return; + } + + void LinearSolverPETSC::solve(const Ref b, Ref x) + { + if (GPU_vec) + VecCreateSeqCUDAWithArrays(PETSC_COMM_WORLD, 1, b.rows() * b.cols(), b.data(), NULL, &b_petsc); + else + VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, b.rows() * b.cols(), b.data(), &b_petsc); + + // KSPSetFromOptions(ksp); + // KSPSetUp(ksp); KSPSolve(ksp, b_petsc, x_petsc); x.resize(b.rows() * b.cols(), 1); + int pidx = 0; for (int i = 0; i < b.rows() * b.cols(); ++i) { @@ -109,16 +178,12 @@ namespace polysolve LinearSolverPETSC::~LinearSolverPETSC() { - // PetscCall(KSPDestroy(&ksp)); - // PetscCall(MatDestroy(&A_petsc)); - // PetscCall(VecDestroy(&b_petsc)); - // PetscCall(VecDestroy(&x_petsc)); - // PetscCall(PetscFinalize()); KSPDestroy(&ksp); MatDestroy(&A_petsc); VecDestroy(&b_petsc); VecDestroy(&x_petsc); - PetscFinalize(); + // PetscFinalize(); + // MISSING: SET A EXTERNAL FINALIZE } } // namespace polysolve diff --git a/src/polysolve/LinearSolverPETSC.hpp b/src/polysolve/LinearSolverPETSC.hpp index 5af76bf..24e366a 100644 --- a/src/polysolve/LinearSolverPETSC.hpp +++ b/src/polysolve/LinearSolverPETSC.hpp @@ -41,8 +41,17 @@ namespace polysolve // Analyze sparsity pattern (sparse) virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override; - // Factorize system matrix (sparse) - virtual void factorize(StiffnessMatrix &A, int dummy) override; + // Factorize system matrix for PETSC(AIJ_CUSPARSE: TRUE OR FALSE) + /*SOLVER INDEX + 0 = PARDISO + 1 = SUPERLU_DIST + 2 = CHOLMOD + 3 = MUMPS + 4 = CUSPARSE + 5 = STRUMPACK + 6 = HYPRE // NOT FULLY IMPLEMENTED YET + */ + virtual void factorize(StiffnessMatrix &A, int AIJ_CUSPARSE, int SOLVER_INDEX) override; // Analyze sparsity pattern (dense, preferred) virtual void analyzePattern(const Eigen::MatrixXd &A, const int precond_num) override; @@ -54,7 +63,7 @@ namespace polysolve virtual std::string name() const override { return "PETSC_Solver"; } protected: - void init(); + int init(); protected: // PETSC variables @@ -63,7 +72,7 @@ namespace polysolve KSP ksp; PC pc; PetscReal norm; - PetscInt its; + PetscInt its, GPU_vec; // Eigen variables Eigen::SparseMatrix A; From 96bffb16be52da331e3a773d7c5235eaa29ff8e8 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Thu, 15 Dec 2022 20:12:44 -0500 Subject: [PATCH 27/29] creating index for solver and adding comments, removed deprecated linking interface --- CMakeLists.txt | 29 +++++++++---------- ...lysolve_target_link_system_libraries.cmake | 28 ------------------ cudatest/cudatest.cpp | 9 +++--- src/polysolve/LinearSolverPETSC.cpp | 22 +++++++++----- src/polysolve/LinearSolverPETSC.hpp | 2 +- tests/test_solver.cpp | 28 ++++++++++++++++++ 6 files changed, 61 insertions(+), 57 deletions(-) delete mode 100644 cmake/polysolve/polysolve_target_link_system_libraries.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index c8a448e..f8acb79 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,11 +47,11 @@ option(POLYSOLVE_WITH_UMFPACK "Enable UmfPack library" ON) option(POLYSOLVE_WITH_SUPERLU "Enable SuperLU library" ON) option(POLYSOLVE_WITH_MKL "Enable MKL library" ON) -option(POLYSOLVE_WITH_CUSOLVER "Enable cuSOLVER library" ON) +option(POLYSOLVE_WITH_CUSOLVER "Enable cuSOLVER library" OFF) -option(POLYSOLVE_WITH_CUDA "Enable CUDA" ON) -option(POLYSOLVE_WITH_PETSC "Enable PETSC" ON) -option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" ON) +option(POLYSOLVE_WITH_CUDA "Enable CUDA" OFF) +option(POLYSOLVE_WITH_PETSC "Enable PETSC" OFF) +option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" OFF) option(POLYSOLVE_WITH_HYPRE "Enable hypre" OFF) option(HYPRE_WITH_MPI "Enable hypre MPI" OFF) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) @@ -98,7 +98,6 @@ include(polysolve_use_colors) # IPC Toolkit utils include(polysolve_prepend_current_path) include(polysolve_set_source_group) -include(polysolve_target_link_system_libraries) # Sort projects inside the solution set_property(GLOBAL PROPERTY USE_FOLDERS ON) @@ -288,20 +287,18 @@ endif() # ############################################################################### -# PETSC solvers +# PETSC solvers (REQUIRES MPI) if(POLYSOLVE_WITH_PETSC) - # include(PETSC PETSc) find_package(MPI) - # TODO: CHECK IF PETSC IS AVAILABLE - # if(TARGET PETSC::PETSc) - target_link_libraries(polysolve PUBLIC -lpetsc MPI::MPI_CXX) - target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PETSC) - message(STATUS "PETSc found!") - - # else() - # message(WARNING "PETSc not found, solver will not be available.") - # endif() + # TODO: CHECK IF PETSC IS PRESENT BEFOREHAND + if(MPI_FOUND) + target_link_libraries(polysolve PUBLIC -lpetsc MPI::MPI_CXX) + target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_PETSC) + message(STATUS "PETSc found!") + else() + message(FATAL_ERROR "PETSc requires MPI!") + endif() endif() # ############################################################################### diff --git a/cmake/polysolve/polysolve_target_link_system_libraries.cmake b/cmake/polysolve/polysolve_target_link_system_libraries.cmake deleted file mode 100644 index 43f8a14..0000000 --- a/cmake/polysolve/polysolve_target_link_system_libraries.cmake +++ /dev/null @@ -1,28 +0,0 @@ -function(polysolve_target_link_system_libraries target) - set(options PRIVATE PUBLIC INTERFACE) - cmake_parse_arguments(TLLS "${options}" "" "" ${ARGN}) - foreach(op ${options}) - if(TLLS_${op}) - set(scope ${op}) - endif() - endforeach(op) - set(libs ${TLLS_UNPARSED_ARGUMENTS}) - - foreach(lib ${libs}) - get_target_property(lib_include_dirs ${lib} INTERFACE_INCLUDE_DIRECTORIES) - if(lib_include_dirs) - if(scope) - target_include_directories(${target} SYSTEM ${scope} ${lib_include_dirs}) - else() - target_include_directories(${target} SYSTEM PRIVATE ${lib_include_dirs}) - endif() - else() - message(WARNING "${lib} doesn't set INTERFACE_INCLUDE_DIRECTORIES. No include_directories set.") - endif() - if(scope) - target_link_libraries(${target} ${scope} ${lib}) - else() - target_link_libraries(${target} ${lib}) - endif() - endforeach() -endfunction(polysolve_target_link_system_libraries) \ No newline at end of file diff --git a/cudatest/cudatest.cpp b/cudatest/cudatest.cpp index c364792..f59cc00 100644 --- a/cudatest/cudatest.cpp +++ b/cudatest/cudatest.cpp @@ -37,12 +37,12 @@ void loadSymmetric(Eigen::SparseMatrix &A, std::string PATH) fin.close(); A.setFromTriplets(triple.begin(), triple.end()); }; - -TEST_CASE("PETSC-STRUMPACK", "[solver]") +#ifdef POLYSOLVE_WITH_PETSC +TEST_CASE("PETSC-DEFAULTWITHCUDA", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; Eigen::SparseMatrix A; - const bool ok = loadMarket(A, path + "/crystm03.mtx"); + const bool ok = loadMarket(A, path + "/A_2.mat"); REQUIRE(ok); auto solver = LinearSolver::create("PETSC_Solver", ""); @@ -53,13 +53,14 @@ TEST_CASE("PETSC-STRUMPACK", "[solver]") x.setZero(); solver->analyzePattern(A, A.rows()); - solver->factorize(A, 0, 5); + solver->factorize(A, 1, 99); solver->solve(b, x); // std::cout<<"Solver error: "< A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + auto solver = LinearSolver::create("PETSC_Solver", ""); + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyzePattern(A, A.rows()); + solver->factorize(A, 1, 99); + solver->solve(b, x); + + const double err = (A * x - b).norm(); + REQUIRE(err < 1e-8); +} #endif \ No newline at end of file From fe05427cd799a746f4cc01caf8fadbe6977323f5 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Sat, 28 Jan 2023 16:48:25 -0500 Subject: [PATCH 28/29] fixed target linking library call --- CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f8acb79..b33c061 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -187,8 +187,7 @@ if(POLYSOLVE_WITH_CUDA) BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) endif() - polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cudart) - polysolve_target_link_system_libraries(polysolve PRIVATE CUDA::cusparse) + target_link_libraries(polysolve PUBLIC -lcudart ${CUDA_cusparse_LIBRARY}) endif() # Hypre (GNU Lesser General Public License) From 5ac5e5602e0d910159c9a1549c388b873edcb660 Mon Sep 17 00:00:00 2001 From: AlexTru96 Date: Tue, 28 Mar 2023 09:53:35 -0400 Subject: [PATCH 29/29] working version for sonic --- CMakeLists.txt | 2 +- cudatest/CMakeLists.txt | 22 +++++++++++----------- src/polysolve/LinearSolverPETSC.cpp | 8 ++++---- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b33c061..83134a0 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -187,7 +187,7 @@ if(POLYSOLVE_WITH_CUDA) BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) endif() - target_link_libraries(polysolve PUBLIC -lcudart ${CUDA_cusparse_LIBRARY}) + target_link_libraries(polysolve PUBLIC -lcudart -lcusparse) endif() # Hypre (GNU Lesser General Public License) diff --git a/cudatest/CMakeLists.txt b/cudatest/CMakeLists.txt index f8f97be..8b0febb 100644 --- a/cudatest/CMakeLists.txt +++ b/cudatest/CMakeLists.txt @@ -1,15 +1,15 @@ -################################################################################ +# ############################################################################### # Tests -################################################################################ +# ############################################################################### set(test_sources main.cpp cudatest.cpp ) add_executable(cuda_test ${test_sources}) -################################################################################ +# ############################################################################### # Required Libraries -################################################################################ +# ############################################################################### include(catch2) target_link_libraries(cuda_test PUBLIC Catch2::Catch2) @@ -17,21 +17,21 @@ target_link_libraries(cuda_test PUBLIC polysolve::polysolve) include(polysolve_warnings) target_link_libraries(cuda_test PRIVATE polysolve::warnings) -################################################################################ + +# ############################################################################### # Register tests -################################################################################ +# ############################################################################### foreach(source IN ITEMS ${test_sources}) - source_group("tests" FILES "${source}") + source_group("tests" FILES "${source}") endforeach() # Register tests set(PARSE_CATCH_TESTS_ADD_TO_CONFIGURE_DEPENDS ON) catch_discover_tests(cuda_test) -################################################################################ +# ############################################################################### # Data -################################################################################ +# ############################################################################### set(DATA_DIR "${CMAKE_SOURCE_DIR}/tests/data") target_compile_definitions(cuda_test PUBLIC -DPOLYSOLVE_DATA_DIR=\"${DATA_DIR}\") -polysolve_target_link_system_libraries(cuda_test PRIVATE CUDA::cudart) -polysolve_target_link_system_libraries(cuda_test PRIVATE CUDA::cusparse) \ No newline at end of file +target_link_libraries(cuda_test PUBLIC -lcudart -lcusparse) \ No newline at end of file diff --git a/src/polysolve/LinearSolverPETSC.cpp b/src/polysolve/LinearSolverPETSC.cpp index c0a85da..1334c4d 100644 --- a/src/polysolve/LinearSolverPETSC.cpp +++ b/src/polysolve/LinearSolverPETSC.cpp @@ -104,11 +104,11 @@ namespace polysolve PCSetType(pc, PCLU); PCFactorSetMatSolverType(pc, MATSOLVERSTRUMPACK); break; - case 6: + // case 6: // TODO : FIX HYPRE PARAMETERS - PCSetType(pc, PCHYPRE); - PCHYPRESetType(pc, "boomeramg"); - break; + // PCSetType(pc, PCHYPRE); + // PCHYPRESetType(pc, "boomeramg"); + // break; default: PCSetType(pc, PCLU); PCFactorSetMatSolverType(pc, MATSOLVERPETSC);