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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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