diff --git a/.gitignore b/.gitignore index 254fae3f..68ce9b4c 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,6 @@ tests/data *.code-workspace PolySolveOptions.cmake .vs -__cmake_systeminformation \ No newline at end of file +__cmake_systeminformation +scaling-experiment +restart-experiments \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 78e1bb93..fc7ecffc 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -96,6 +96,7 @@ option(POLYSOLVE_WITH_PARDISO "Enable Pardiso library" option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) option(POLYSOLVE_WITH_SPECTRA "Enable Spectra library" ON) +option(POLYSOLVE_WITH_STRUMPACK "Enable STRUMPACK (experimental CPUHybrid subdomain factorization)" OFF) # Sanitizer options option(POLYSOLVE_SANITIZE_ADDRESS "Sanitize Address" OFF) @@ -304,11 +305,12 @@ if(POLYSOLVE_WITH_CUDA) target_link_libraries(polysolve PUBLIC CCCL::CCCL) endif() -# Cuda runtime and cuSparse +# Cuda runtime, cuSparse and cuBLAS (cuBLAS is used by GPUHybridSolver's +# contact_patch_schwarz dense batched subdomain solve) if (POLYSOLVE_WITH_CUDA) find_package(CUDAToolkit REQUIRED) - target_link_libraries(polysolve PUBLIC CUDA::cudart CUDA::cusparse) - target_link_libraries(polysolve_linear PUBLIC CUDA::cudart CUDA::cusparse) + target_link_libraries(polysolve PUBLIC CUDA::cudart CUDA::cusparse CUDA::cublas) + target_link_libraries(polysolve_linear PUBLIC CUDA::cudart CUDA::cusparse CUDA::cublas) endif() # ------ @@ -432,6 +434,13 @@ if(POLYSOLVE_WITH_SPECTRA) target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_SPECTRA) endif() +# STRUMPACK (BSD 3-Clause), used only by CPUHybrid's experimental use_strumpack +if(POLYSOLVE_WITH_STRUMPACK) + include(strumpack) + target_link_libraries(polysolve_linear PRIVATE STRUMPACK::strumpack) + target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_STRUMPACK) +endif() + # cuSolver solvers if(POLYSOLVE_WITH_CUSOLVER) if (NOT POLYSOLVE_WITH_CUDA) diff --git a/PolySolveOptions.cmake.sample b/PolySolveOptions.cmake.sample index c23eca08..6993251f 100644 --- a/PolySolveOptions.cmake.sample +++ b/PolySolveOptions.cmake.sample @@ -34,6 +34,7 @@ # option(POLYSOLVE_WITH_HYPRE "Enable hypre" ON) # option(POLYSOLVE_WITH_AMGCL "Use AMGCL" ON) # option(POLYSOLVE_WITH_SPECTRA "Enable Spectra library" ON) +# option(POLYSOLVE_WITH_STRUMPACK "Enable STRUMPACK (experimental CPUHybrid subdomain factorization)" OFF) # Options for third-party libraries # option(EIGEN_WITH_MKL "Use Eigen with MKL" ON) diff --git a/cmake/recipes/amgcl.cmake b/cmake/recipes/amgcl.cmake index 53849d5f..8a26ca05 100644 --- a/cmake/recipes/amgcl.cmake +++ b/cmake/recipes/amgcl.cmake @@ -39,6 +39,22 @@ function(amgcl_import_target) # Prefer Config mode before Module mode to prevent lib from loading its own FindXXX.cmake set(CMAKE_FIND_PACKAGE_PREFER_CONFIG TRUE) + # AMGCL's own build optionally calls find_package(MPI) to enable its + # distributed-memory backend. When nano-mpi is in play, CMAKE_FIND_PACKAGE_PREFER_CONFIG + # above makes that call skip Module mode entirely -- bypassing + # CMAKE_MODULE_PATH and nano-mpi's FindMPI shim -- and resolve via Config + # mode straight to a real system MPI instead (Homebrew's Open MPI ships + # an MPIConfig.cmake). Since MPI::MPI_C is a single GLOBAL imported + # target, that unconditionally overwrites its INTERFACE_INCLUDE_DIRECTORIES + # for everyone, including hypre, silently pointing their #include + # at a real MPI nano-mpi cannot run against. AMGCL's distributed backend + # would need a real launcher anyway, which nano-mpi does not provide, so + # disable that find_package(MPI) call here rather than let it corrupt the + # shared target. + if(TARGET nanompi::nanompi) + set(CMAKE_DISABLE_FIND_PACKAGE_MPI TRUE) + endif() + # Ready to include third-party lib include(CPM) CPMAddPackage( diff --git a/cmake/recipes/argparse.cmake b/cmake/recipes/argparse.cmake new file mode 100644 index 00000000..0e116e1b --- /dev/null +++ b/cmake/recipes/argparse.cmake @@ -0,0 +1,16 @@ +if(TARGET argparse::argparse) + return() +endif() + +message(STATUS "Third-party: creating target 'argparse::argparse'") + +include(CPM) +CPMAddPackage( + NAME argparse + GITHUB_REPOSITORY p-ranav/argparse + GIT_TAG v3.2 +) + +if(TARGET argparse) + set_target_properties(argparse PROPERTIES SYSTEM ON) +endif() diff --git a/cmake/recipes/metis.cmake b/cmake/recipes/metis.cmake new file mode 100644 index 00000000..bd557f7e --- /dev/null +++ b/cmake/recipes/metis.cmake @@ -0,0 +1,69 @@ +# METIS (https://github.com/KarypisLab/METIS), with GKlib (https://github.com/KarypisLab/GKlib) +# License: Apache 2.0 + +if(TARGET METIS::metis) + return() +endif() + +message(STATUS "Third-party: creating target 'METIS::metis'") + +# Neither project can be add_subdirectory()'d: METIS's CMake expects GKlib to +# be installed under CMAKE_INSTALL_PREFIX already, and a metis.h that its +# `make config` step fills in. So both are fetched as plain sources and +# compiled here into a single static library. +include(CPM) +CPMAddPackage( + NAME gklib + GITHUB_REPOSITORY KarypisLab/GKlib + # Pinned: GKlib's master has since moved its sources into src/ and include/. + GIT_TAG 8bd6bad750b2b0d90800c632cf18e8ee93ad72d7 + DOWNLOAD_ONLY YES +) +CPMAddPackage( + NAME metis + GITHUB_REPOSITORY KarypisLab/METIS + GIT_TAG v5.2.1 + DOWNLOAD_ONLY YES +) + +enable_language(C) + +# metis.h ships with its index and real widths commented out, for `make config` +# to fill in. Fill in 32 bits, matching the int indices of the Eigen matrices +# every caller here passes in. file(CONFIGURE) rewrites the header only when +# its content changes, so a reconfigure does not trigger a rebuild. +file(READ "${metis_SOURCE_DIR}/include/metis.h" METIS_HEADER) +string(REPLACE "//#define IDXTYPEWIDTH 32" "#define IDXTYPEWIDTH 32" METIS_HEADER "${METIS_HEADER}") +string(REPLACE "//#define REALTYPEWIDTH 32" "#define REALTYPEWIDTH 32" METIS_HEADER "${METIS_HEADER}") +file(CONFIGURE OUTPUT "${metis_BINARY_DIR}/include/metis.h" CONTENT "${METIS_HEADER}" @ONLY) + +file(GLOB GKLIB_SOURCES "${gklib_SOURCE_DIR}/*.c") +file(GLOB METIS_SOURCES "${metis_SOURCE_DIR}/libmetis/*.c") + +add_library(METIS_metis STATIC ${GKLIB_SOURCES} ${METIS_SOURCES}) +add_library(METIS::metis ALIAS METIS_metis) +set_target_properties(METIS_metis PROPERTIES FOLDER third_party/metis) + +target_include_directories(METIS_metis + PUBLIC + "${metis_BINARY_DIR}/include" + PRIVATE + "${gklib_SOURCE_DIR}" + "${metis_SOURCE_DIR}/libmetis" +) + +# The definitions GKlib's own build passes (see its GKlibSystem.cmake), minus +# -march=native and -Werror. +if(MSVC) + target_include_directories(METIS_metis PRIVATE "${gklib_SOURCE_DIR}/win32") + target_compile_definitions(METIS_metis PRIVATE WIN32 MSC _CRT_SECURE_NO_DEPRECATE USE_GKREGEX "__thread=__declspec(thread)") +elseif(MINGW) + target_compile_definitions(METIS_metis PRIVATE USE_GKREGEX) +else() + target_compile_definitions(METIS_metis PRIVATE LINUX _FILE_OFFSET_BITS=64) +endif() +target_compile_definitions(METIS_metis PRIVATE NDEBUG NDEBUG2) + +if(UNIX) + target_link_libraries(METIS_metis PRIVATE m) +endif() diff --git a/cmake/recipes/nanompi.cmake b/cmake/recipes/nanompi.cmake index 204f8c86..7c75ba38 100644 --- a/cmake/recipes/nanompi.cmake +++ b/cmake/recipes/nanompi.cmake @@ -14,6 +14,17 @@ list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/nanompi-as-mpi/") list(REMOVE_DUPLICATES CMAKE_MODULE_PATH) +# Promote to CACHE so the shim is visible in every directory scope, not just +# ones that inherit this list(APPEND) through add_subdirectory(). A dependency +# added from a scope that never saw it (e.g. AMGCL's own optional +# find_package(MPI) call) would otherwise fall through to CMake's real +# FindMPI, which unconditionally overwrites MPI::MPI_C's +# INTERFACE_INCLUDE_DIRECTORIES on every invocation (see _MPI_create_imported_target +# in CMake's own FindMPI.cmake) -- silently re-pointing everyone's +# #include , hypre included, at a real system MPI nano-mpi cannot +# actually run against. +set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" CACHE STRING "" FORCE) + if(TARGET nanompi::nanompi) return() endif() diff --git a/cmake/recipes/strumpack.cmake b/cmake/recipes/strumpack.cmake new file mode 100644 index 00000000..8c5c41a9 --- /dev/null +++ b/cmake/recipes/strumpack.cmake @@ -0,0 +1,134 @@ +# STRUMPACK (https://github.com/pghysels/STRUMPACK) +# License: BSD 3-Clause (LBNL) +# +# Sparse multifrontal LU, with optional low-rank (BLR, HSS) compression of its +# frontal matrices. CPUHybridSolver can use it to factorize its subdomains +# (the experimental use_strumpack option). + +if(TARGET STRUMPACK::strumpack) + return() +endif() + +message(STATUS "Third-party: creating target 'STRUMPACK::strumpack'") + +# STRUMPACK's project() declares Fortran (for BLAS name mangling and its +# Fortran interface), so it cannot even configure without a Fortran compiler. +include(CheckLanguage) +check_language(Fortran) +if(NOT CMAKE_Fortran_COMPILER) + message(FATAL_ERROR + "STRUMPACK needs a Fortran compiler, and none was found. Install one " + "(e.g. gfortran) or turn POLYSOLVE_WITH_STRUMPACK off.") +endif() + +include(metis) +include(blas) +include(lapack) + +include(CPM) +CPMAddPackage( + NAME strumpack + GITHUB_REPOSITORY pghysels/STRUMPACK + GIT_TAG v8.0.0 + DOWNLOAD_ONLY YES +) + +# Unlike the other recipes, STRUMPACK is not add_subdirectory()'d: its CMake +# only works as the top-level project. It finds its own modules and test +# sources through CMAKE_SOURCE_DIR, which would point at polysolve instead, and +# unconditionally adds its ~200 tests to CTest, which would land in polysolve's. +# So it is configured and built in a CMake run of its own, and only the +# resulting static library is imported below. +set(STRUMPACK_INSTALL_DIR "${FETCHCONTENT_BASE_DIR}/strumpack-install") +set(STRUMPACK_LIBRARY "${STRUMPACK_INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}strumpack${CMAKE_STATIC_LIBRARY_SUFFIX}") + +# Follow polysolve's build type; with none set, STRUMPACK would be built +# unoptimized. +if(CMAKE_BUILD_TYPE) + set(STRUMPACK_BUILD_TYPE ${CMAKE_BUILD_TYPE}) +else() + set(STRUMPACK_BUILD_TYPE Release) +endif() + +get_target_property(STRUMPACK_METIS_INCLUDE_DIR METIS::metis INTERFACE_INCLUDE_DIRECTORIES) + +set(STRUMPACK_CMAKE_ARGS + -DCMAKE_BUILD_TYPE=${STRUMPACK_BUILD_TYPE} + -DCMAKE_INSTALL_PREFIX= + -DCMAKE_INSTALL_LIBDIR=lib + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} + -DCMAKE_Fortran_COMPILER=${CMAKE_Fortran_COMPILER} + -DCMAKE_POSITION_INDEPENDENT_CODE=ON + -DBUILD_SHARED_LIBS=OFF + + # Each subdomain is factorized by a single nanompi rank -- a thread of + # this process -- so STRUMPACK runs sequentially inside it, with neither + # MPI nor OpenMP of its own. + -DSTRUMPACK_USE_MPI=OFF + -DSTRUMPACK_USE_OPENMP=OFF + -DSTRUMPACK_USE_CUDA=OFF + -DSTRUMPACK_USE_HIP=OFF + -DSTRUMPACK_USE_SYCL=OFF + + # No optional third-party libraries: BLR and HSS compression are built in. + -DTPL_ENABLE_SLATE=OFF + -DTPL_ENABLE_PARMETIS=OFF + -DTPL_ENABLE_SCOTCH=OFF + -DTPL_ENABLE_PTSCOTCH=OFF + -DTPL_ENABLE_BPACK=OFF + -DTPL_ENABLE_ZFP=OFF + -DTPL_ENABLE_SZ3=OFF + -DTPL_ENABLE_MAGMA=OFF + -DTPL_ENABLE_KBLAS=OFF + -DTPL_ENABLE_COMBBLAS=OFF + -DTPL_ENABLE_PAPI=OFF + -DTPL_ENABLE_MATLAB=OFF + + # METIS, the one required dependency, is the one built by the metis + # recipe. Setting METIS_INCLUDE_DIR outright keeps STRUMPACK's FindMETIS + # from settling on a system metis.h instead. The library itself is only + # needed once polysolve links, so it does not have to be built yet. + -DMETIS_INCLUDE_DIR=${STRUMPACK_METIS_INCLUDE_DIR} + -DTPL_METIS_LIBRARIES=$ + + # BLAS and LAPACK are left for STRUMPACK's configure to find: it uses + # them only to check that a test program links. What STRUMPACK actually + # calls is whatever polysolve links it with, below. +) + +# Build the library only: STRUMPACK's default target also builds its test +# programs. With a Makefile generator, $(MAKE) lets the sub-build share the +# outer make's job slots instead of running serially. +if(CMAKE_GENERATOR MATCHES "Makefiles") + set(STRUMPACK_BUILD_COMMAND $(MAKE) strumpack) +else() + set(STRUMPACK_BUILD_COMMAND ${CMAKE_COMMAND} --build --target strumpack --config ${STRUMPACK_BUILD_TYPE}) +endif() + +include(ExternalProject) +ExternalProject_Add(strumpack_build + SOURCE_DIR "${strumpack_SOURCE_DIR}" + BINARY_DIR "${strumpack_BINARY_DIR}" + INSTALL_DIR "${STRUMPACK_INSTALL_DIR}" + CMAKE_ARGS ${STRUMPACK_CMAKE_ARGS} + BUILD_COMMAND ${STRUMPACK_BUILD_COMMAND} + # `cmake --install`, unlike the install target, does not first build + # everything else. + INSTALL_COMMAND ${CMAKE_COMMAND} --install --config ${STRUMPACK_BUILD_TYPE} + BUILD_BYPRODUCTS "${STRUMPACK_LIBRARY}" +) + +# An imported target's include directory has to exist at generate time, well +# before the ExternalProject has installed anything into it. +file(MAKE_DIRECTORY "${STRUMPACK_INSTALL_DIR}/include") + +add_library(STRUMPACK::strumpack STATIC IMPORTED GLOBAL) +set_target_properties(STRUMPACK::strumpack PROPERTIES + IMPORTED_LOCATION "${STRUMPACK_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${STRUMPACK_INSTALL_DIR}/include" +) +# Link polysolve's BLAS and LAPACK (MKL, Accelerate, ...), not whichever ones +# STRUMPACK's configure happened to find. +target_link_libraries(STRUMPACK::strumpack INTERFACE METIS::metis LAPACK::LAPACK BLAS::BLAS) +add_dependencies(STRUMPACK::strumpack strumpack_build) diff --git a/linear-solver-spec.json b/linear-solver-spec.json index fbf698c4..23309d44 100644 --- a/linear-solver-spec.json +++ b/linear-solver-spec.json @@ -18,7 +18,8 @@ "AMGCL", "MAS", "CPUHybrid", - "GPUHybrid" + "GPUHybrid", + "RestartedAMGPCG" ], "doc": "Settings for the linear solver." }, @@ -53,7 +54,8 @@ "Eigen::MINRES", "MAS", "CPUHybrid", - "GPUHybrid" + "GPUHybrid", + "RestartedAMGPCG" ] }, { @@ -148,8 +150,7 @@ "tolerance", "theta", "nodal_coarsening", - "interp_rbms", - "dimension" + "interp_rbms" ], "doc": "Settings for the Hypre solver." }, @@ -288,12 +289,6 @@ "type": "bool", "doc": "Whether or not to include nodal coarsening options." }, - { - "pointer": "/Hypre/dimension", - "default": 1, - "type": "int", - "doc": "Dimension of problem." - }, { "pointer": "/AMGCL/solver", "default": null, @@ -461,7 +456,6 @@ "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -470,17 +464,6 @@ ], "doc": "Settings for the MAS solver." }, - { - "pointer": "/MAS/block_dim", - "default": 1, - "type": "int", - "options": [ - 1, - 2, - 3 - ], - "doc": "Block size of the BSR matrix PCG sovler internally use. Choose 3 for 3d problem, 2 for 2d, etc." - }, { "pointer": "/MAS/max_iter", "default": 10000, @@ -516,7 +499,6 @@ "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -529,16 +511,19 @@ "gmm_tol", "max_gmm_iterations", "conditioning_threshold", - "additive_mode" + "additive_mode", + "select_bad_dofs_from_l1_norm", + "subdomain_selection_strategy", + "contact_patch_schwarz", + "contact_patch_remove_covered", + "ras_correction", + "use_strumpack", + "strumpack_compression", + "strumpack_compression_rel_tol", + "detailed_log" ], "doc": "Settings for the CPUHybrid solver." }, - { - "pointer": "/CPUHybrid/block_dim", - "default": 1, - "type": "int", - "doc": "Block size of the system matrix. " - }, { "pointer": "/CPUHybrid/max_iter", "default": 10000, @@ -617,12 +602,77 @@ "type": "bool", "doc": "Apply subspace correction additively. " }, + { + "pointer": "/CPUHybrid/select_bad_dofs_from_l1_norm", + "default": true, + "type": "bool", + "doc": "Select the problematic subspace from the L1 row norm heuristic. If false (AMGF mode), the problematic subspace is taken only from set_problematic_dofs (and, if enabled, contact_patch_schwarz) instead." + }, + { + "pointer": "/CPUHybrid/subdomain_selection_strategy", + "default": "gmm", + "type": "string", + "options": [ + "knee", + "gmm", + "fd", + "cost", + "aposteriori" + ], + "doc": "Strategy used to threshold row norms into the problematic subspace: gmm (Gaussian mixture model), knee (max-distance-to-chord knee point), fd (largest finite-difference jump), cost (minimize estimated solve cost), or aposteriori (not yet implemented)." + }, + { + "pointer": "/CPUHybrid/contact_patch_schwarz", + "default": false, + "type": "bool", + "doc": "Additionally give each contact patch (see Solver::contact_patches) its own overlapping subdomain, combined additively (Schwarz) with the rest of the problematic subspace." + }, + { + "pointer": "/CPUHybrid/contact_patch_remove_covered", + "default": false, + "type": "bool", + "doc": "Only used when contact_patch_schwarz is enabled. Before factorizing, drop any contact patch that is fully contained in another patch, since it would only ever add redundant correction; logs the number of patches before and after." + }, + { + "pointer": "/CPUHybrid/ras_correction", + "default": false, + "type": "bool", + "doc": "Restricted Additive Schwarz correction: use the same problematic dof set as AMGF/select_bad_dofs, but build one subdomain per dof (every dof with a nonzero entry in that dof's row) and let each subdomain correct only the single dof it was built around, instead of merging dofs into shared clusters. Mutually exclusive with contact_patch_schwarz." + }, + { + "pointer": "/CPUHybrid/use_strumpack", + "default": false, + "type": "bool", + "doc": "Experimental: factorize every subdomain larger than 1000 dofs with STRUMPACK's sparse multifrontal LU instead of Pardiso/Accelerate (smaller subdomains keep using SimplicialLDLT). Requires polysolve to be built with POLYSOLVE_WITH_STRUMPACK." + }, + { + "pointer": "/CPUHybrid/strumpack_compression", + "default": "none", + "type": "string", + "options": [ + "none", + "blr", + "hss" + ], + "doc": "Only used when use_strumpack is enabled. Low-rank compression of STRUMPACK's frontal matrices: none (exact LU), blr (block low-rank) or hss (hierarchically semi-separable). STRUMPACK only compresses fronts whose separator has at least 512 (blr) or 1000 (hss) dofs, so small subdomains are factorized exactly either way. A compressed factorization is not exactly symmetric, so PCG's convergence guarantees no longer strictly hold." + }, + { + "pointer": "/CPUHybrid/strumpack_compression_rel_tol", + "default": 1e-4, + "type": "float", + "doc": "Only used when strumpack_compression is not none. Relative tolerance of the low-rank compression." + }, + { + "pointer": "/CPUHybrid/detailed_log", + "default": true, + "type": "bool", + "doc": "If false, skip the per-PCG-iteration timing logs (matmul, amg_v_cycle, subdomain_solve, pcg_iter) so a long solve doesn't flood the log; every other log line is unaffected." + }, { "pointer": "/GPUHybrid", "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -635,16 +685,14 @@ "gmm_tol", "max_gmm_iterations", "conditioning_threshold", - "additive_mode" + "additive_mode", + "select_bad_dofs_from_l1_norm", + "subdomain_selection_strategy", + "contact_patch_schwarz", + "detailed_log" ], "doc": "Settings for the GPUHybrid solver." }, - { - "pointer": "/GPUHybrid/block_dim", - "default": 1, - "type": "int", - "doc": "Block size of the system matrix. " - }, { "pointer": "/GPUHybrid/max_iter", "default": 10000, @@ -722,5 +770,87 @@ "default": false, "type": "bool", "doc": "Apply subspace correction additively. " + }, + { + "pointer": "/GPUHybrid/select_bad_dofs_from_l1_norm", + "default": true, + "type": "bool", + "doc": "Select the problematic subspace from the L1 row norm heuristic. If false (AMGF mode), the problematic subspace is taken only from set_problematic_dofs (and, if enabled, contact_patch_schwarz) instead." + }, + { + "pointer": "/GPUHybrid/subdomain_selection_strategy", + "default": "gmm", + "type": "string", + "options": [ + "knee", + "gmm", + "fd", + "cost", + "aposteriori" + ], + "doc": "Strategy used to threshold row norms into the problematic subspace: gmm (Gaussian mixture model), knee (max-distance-to-chord knee point), fd (largest finite-difference jump), cost (minimize estimated solve cost), or aposteriori (not yet implemented)." + }, + { + "pointer": "/GPUHybrid/contact_patch_schwarz", + "default": false, + "type": "bool", + "doc": "Additionally give each contact patch (see Solver::contact_patches) its own overlapping subdomain, combined additively (Schwarz) with the rest of the problematic subspace. Unlike the generic problematic subspace, contact patches are factorized and solved with a fixed-size (64 dof) dense batched GPU solve; a contact patch with more than 64 dofs throws." + }, + { + "pointer": "/GPUHybrid/detailed_log", + "default": true, + "type": "bool", + "doc": "If false, skip the per-PCG-iteration timing logs (matmul, amg_v_cycle, subdomain_solve, contact_patch_dense_subdomain_solve, pcg_iter) so a long solve doesn't flood the log; every other log line is unaffected." + }, + { + "pointer": "/RestartedAMGPCG", + "default": null, + "type": "object", + "optional": [ + "max_iter", + "relative_tolerance", + "absolute_tolerance", + "theta", + "restart_interval", + "detailed_log" + ], + "doc": "Settings for the RestartedAMGPCG solver: GPU BoomerAMG-preconditioned CG (AMG set up as in GPUHybrid, without its subspace correction), periodically restarted from the true residual recomputed in double-double precision." + }, + { + "pointer": "/RestartedAMGPCG/max_iter", + "default": 10000, + "type": "int", + "doc": "Maximum number of PCG iterations, summed over all restarts." + }, + { + "pointer": "/RestartedAMGPCG/relative_tolerance", + "default": 1e-10, + "type": "float", + "doc": "Relative tolerance for solver convergence, checked against the double-double precision residual." + }, + { + "pointer": "/RestartedAMGPCG/absolute_tolerance", + "default": 0, + "type": "float", + "doc": "Absolute tolerance for solver convergence, checked against the double-double precision residual." + }, + { + "pointer": "/RestartedAMGPCG/theta", + "default": 0.5, + "type": "float", + "doc": "Strong threshold for AMG." + }, + { + "pointer": "/RestartedAMGPCG/restart_interval", + "default": 50, + "type": "int", + "min": 1, + "doc": "Every restart_interval PCG iterations (or sooner, whenever the recursively-updated residual reaches the tolerance), recompute the true residual b - Ax in double-double precision and restart PCG from it. Convergence is only declared on this recomputed residual." + }, + { + "pointer": "/RestartedAMGPCG/detailed_log", + "default": true, + "type": "bool", + "doc": "If false, skip the per-PCG-iteration timing logs (matmul, amg_v_cycle, pcg_iter) so a long solve doesn't flood the log; every other log line, including the per-restart residual log, is unaffected." } ] diff --git a/nonlinear-solver-spec.json b/nonlinear-solver-spec.json index 8981feaa..9e056f84 100644 --- a/nonlinear-solver-spec.json +++ b/nonlinear-solver-spec.json @@ -23,7 +23,9 @@ "ADAM", "StochasticADAM", "StochasticGradientDescent", + "NonlinearCG", "box_constraints", + "simple_linear", "advanced" ], "doc": "Settings for nonlinear solver. Interior-loop linear solver settings are defined in the solver/linear section." @@ -35,7 +37,9 @@ "options": [ "Newton", "DenseNewton", + "NewtonCG", "GradientDescent", + "NonlinearCG", "ADAM", "StochasticADAM", "StochasticGradientDescent", @@ -314,6 +318,34 @@ "type": "float", "doc": "Probability of erasing a component on the gradient for StochasticGradientDescent." }, + { + "pointer": "/NonlinearCG", + "default": null, + "type": "object", + "optional": [ + "formula", + "restart_frequency" + ], + "doc": "Options for (unpreconditioned) Nonlinear Conjugate Gradient." + }, + { + "pointer": "/NonlinearCG/formula", + "default": "PolakRibiere", + "type": "string", + "options": [ + "FletcherReeves", + "PolakRibiere", + "HestenesStiefel" + ], + "doc": "Formula used to compute the beta coefficient combining the previous and current search directions. A negative beta always triggers a restart to steepest descent." + }, + { + "pointer": "/NonlinearCG/restart_frequency", + "default": 0, + "type": "int", + "min": 0, + "doc": "Number of iterations between forced restarts to steepest descent. A value of 0 defaults to the number of degrees of freedom of the problem." + }, { "pointer": "/solver", "type": "list", @@ -373,6 +405,18 @@ ], "doc": "Options for regularized projected Newton." }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "NewtonCG", + "required": [ + "type" + ], + "optional": [ + "force_psd_projection" + ], + "doc": "Inexact Newton's method that requires the linear solver to be CPUHybrid or GPUHybrid. The linear solver's relative tolerance is set to min(0.5, sqrt(||grad||)) at each iteration and its residual tolerance is treated as infinite, so the (possibly inexact) CG direction is always accepted." + }, { "pointer": "/solver/*", "type": "object", @@ -448,6 +492,19 @@ ], "doc": "Options for Stochastic Gradient Descent." }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "NonlinearCG", + "required": [ + "type" + ], + "optional": [ + "formula", + "restart_frequency" + ], + "doc": "Options for (unpreconditioned) Nonlinear Conjugate Gradient." + }, { "pointer": "/solver/*", "type": "object", @@ -512,8 +569,10 @@ "DenseRegularizedNewton", "RegularizedProjectedNewton", "DenseRegularizedProjectedNewton", + "NewtonCG", "GradientDescent", "StochasticGradientDescent", + "NonlinearCG", "ADAM", "StochasticADAM", "L-BFGS", @@ -545,12 +604,36 @@ "type": "float", "doc": "Regulariztion weight increment." }, + { + "pointer": "/solver/*/force_psd_projection", + "default": false, + "type": "bool", + "doc": "NewtonCG only: force the Hessian to be projected to PSD before each (inexact, CG-based) linear solve." + }, { "pointer": "/solver/*/erase_component_probability", "default": 0.3, "type": "float", "doc": "Probability of erasing a component on the gradient for stochastic solvers." }, + { + "pointer": "/solver/*/formula", + "default": "PolakRibiere", + "type": "string", + "options": [ + "FletcherReeves", + "PolakRibiere", + "HestenesStiefel" + ], + "doc": "Formula used to compute the beta coefficient combining the previous and current search directions. A negative beta always triggers a restart to steepest descent." + }, + { + "pointer": "/solver/*/restart_frequency", + "default": 0, + "type": "int", + "min": 0, + "doc": "Number of iterations between forced restarts to steepest descent. A value of 0 defaults to the number of degrees of freedom of the problem." + }, { "pointer": "/solver/*/history_size", "default": 6, @@ -686,6 +769,12 @@ "min_value": 0, "doc": "Relative tolerance on E to switch to approximate." }, + { + "pointer": "/simple_linear", + "default": "Eigen::SimplicialLDLT", + "type": "string", + "doc": "Linear solver used internally by callers for small, auxiliary linear systems unrelated to the main Hessian/energy solve (e.g. a constraint-elimination projection), as opposed to the (possibly heavyweight, e.g. Hybrid/AMGF) solver configured in solver/linear for the main solve. Not read by polysolve::nonlinear::Solver itself; exposed here so callers building such auxiliary solves have a single, user-configurable place to look." + }, { "pointer": "/box_constraints", "type": "object", diff --git a/src/polysolve/linear/AMGCL.cpp b/src/polysolve/linear/AMGCL.cpp index cc638648..36704e70 100644 --- a/src/polysolve/linear/AMGCL.cpp +++ b/src/polysolve/linear/AMGCL.cpp @@ -107,11 +107,6 @@ namespace polysolve::linear { if (params.contains("AMGCL")) { - // Specially named parameters to match other solvers - if (params["AMGCL"].contains("block_size")) - { - block_size_ = params["AMGCL"]["block_size"]; - } if (block_size_ == 2) { block2_solver_.set_parameters(params); @@ -127,6 +122,11 @@ namespace polysolve::linear } } + void AMGCL::set_block_size(int block_size) + { + block_size_ = block_size; + } + void AMGCL::get_info(json ¶ms) const { if (block_size_ == 2) diff --git a/src/polysolve/linear/AMGCL.hpp b/src/polysolve/linear/AMGCL.hpp index 2518358f..9485d6c2 100644 --- a/src/polysolve/linear/AMGCL.hpp +++ b/src/polysolve/linear/AMGCL.hpp @@ -115,6 +115,9 @@ namespace polysolve::linear // Set solver parameters virtual void set_parameters(const json ¶ms) override; + // Set block size for multigrid solvers + virtual void set_block_size(int block_size) override; + // Retrieve information virtual void get_info(json ¶ms) const override; diff --git a/src/polysolve/linear/CMakeLists.txt b/src/polysolve/linear/CMakeLists.txt index 03916f9d..d240576a 100644 --- a/src/polysolve/linear/CMakeLists.txt +++ b/src/polysolve/linear/CMakeLists.txt @@ -23,6 +23,8 @@ if(POLYSOLVE_WITH_CUDA) list(APPEND SOURCES GPUHybridSolver.cu GPUHybridSolver.hpp + RestartedAMGPCG.cu + RestartedAMGPCG.hpp cuDSS.hpp cuDSS.cpp ) diff --git a/src/polysolve/linear/CPUHybridSolver.cpp b/src/polysolve/linear/CPUHybridSolver.cpp index 39cf9d6e..6a8820f6 100644 --- a/src/polysolve/linear/CPUHybridSolver.cpp +++ b/src/polysolve/linear/CPUHybridSolver.cpp @@ -5,7 +5,12 @@ #include "hybrid_utils/DisjointSet.hpp" +#include +#include #include +#include +#include +#include #include #include @@ -18,6 +23,10 @@ #include #endif +#ifdef POLYSOLVE_WITH_STRUMPACK +#include +#endif + #include #include #include @@ -37,8 +46,116 @@ namespace polysolve::linear { return std::chrono::duration(clock::now() - begin).count(); } + +#ifdef POLYSOLVE_WITH_STRUMPACK + // A subdomain factorized by STRUMPACK, under use_strumpack. Only the + // rank that owns the subdomain -- a thread of this process -- ever + // touches it, and STRUMPACK is built without MPI or OpenMP, so it runs + // sequentially on that thread. + // + // Note that STRUMPACK's solver installs a process-wide new-handler + // for its lifetime, which aborts on out-of-memory rather than + // throwing std::bad_alloc. + class StrumpackWrapper : public AbstractSolver + { + strumpack::StrumpackSparseSolver solver{/*verbose=*/false}; + + public: + StrumpackWrapper(CPUHybridSolver::StrumpackCompression compression, double compression_rel_tol) + { + auto &opts = solver.options(); + + // Apply the factorization exactly once per solve. The default + // (AUTO) wraps it in iterative refinement -- or, once + // compressed, in GMRES -- which would make this preconditioner + // nonlinear under PCG. + opts.set_Krylov_solver(strumpack::KrylovSolver::DIRECT); + // Subdomains are principal submatrices of an SPD system, so + // their diagonal needs no help from MC64's row permutation, + // which would only throw away their symmetry. + opts.set_matching(strumpack::MatchingJob::NONE); + + switch (compression) + { + case CPUHybridSolver::StrumpackCompression::NONE: + opts.set_compression(strumpack::CompressionType::NONE); + break; + case CPUHybridSolver::StrumpackCompression::BLR: + opts.set_compression(strumpack::CompressionType::BLR); + break; + case CPUHybridSolver::StrumpackCompression::HSS: + opts.set_compression(strumpack::CompressionType::HSS); + break; + } + opts.set_compression_rel_tol(compression_rel_tol); + } + + void compute(const Eigen::SparseMatrix &A) override + { + // STRUMPACK takes CSR, and keeps a copy of its own. + Eigen::SparseMatrix A_csr = A; + A_csr.makeCompressed(); + + // Don't promise a symmetric pattern: an assembled matrix can + // store an explicit zero on one side of the diagonal only, + // and STRUMPACK would trust the promise. + solver.set_csr_matrix(A_csr.rows(), A_csr.outerIndexPtr(), A_csr.innerIndexPtr(), A_csr.valuePtr(), /*symmetric_pattern=*/false); + + // The reordering is METIS nested dissection, and METIS keeps + // its random number generator in a process-wide static that + // the other ranks' subdomains, reordered at the same time, + // would race on. So reorder one subdomain at a time; METIS + // reseeds on every call, which also makes the orderings + // deterministic. The numeric factorization stays concurrent. + strumpack::ReturnCode status; + { + static std::mutex metis_mutex; + std::lock_guard lock(metis_mutex); + status = solver.reorder(); + } + if (status == strumpack::ReturnCode::SUCCESS) + { + status = solver.factor(); + } + if (status != strumpack::ReturnCode::SUCCESS) + { + SPDLOG_ERROR("[CPUHybrid] STRUMPACK failed to factorize a subdomain of size {} (return code {})", A.rows(), static_cast(status)); + } + } + + Eigen::VectorXd solve(const Eigen::VectorXd &b) override + { + Eigen::VectorXd x(b.size()); + solver.solve(b.data(), x.data()); + return x; + } + }; +#endif } // namespace + // Workers (nanompi ranks 1..N-1, each a thread of this process) run the + // very same collective code paths as rank 0 -- see run_worker_loop -- + // so an unguarded SPDLOG_INFO below would print once per rank. Route + // rank-scoped logging through this macro instead of muting the logger + // globally: spdlog::set_level() mutates the single process-wide default + // logger, so calling it from a worker thread silences rank 0 as well. +#define CPUHYBRID_LOG_INFO(...) \ + do \ + { \ + if (myid == 0) \ + SPDLOG_INFO(__VA_ARGS__); \ + } while (0) + + // Per-PCG-iteration timing logs (matmul, amg_v_cycle, subdomain_solve, + // pcg_iter) go through this instead of CPUHYBRID_LOG_INFO, so they can be + // silenced with detailed_log=false without affecting any other log line. +#define CPUHYBRID_LOG_ITER(...) \ + do \ + { \ + if (myid == 0 && detailed_log) \ + SPDLOG_INFO(__VA_ARGS__); \ + } while (0) + CPUHybridSolver::CPUHybridSolver() { // Ranks are threads of this process. The first solver constructed @@ -63,11 +180,7 @@ namespace polysolve::linear HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST); HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST); - if (myid != 0) - { - spdlog::set_level(spdlog::level::off); - } - else + if (myid == 0) { spdlog::flush_on(spdlog::level::info); } @@ -131,10 +244,6 @@ namespace polysolve::linear { theta = shared_params["CPUHybrid"]["theta"]; } - if (shared_params["CPUHybrid"].contains("block_dim")) - { - dimension_ = shared_params["CPUHybrid"]["block_dim"]; - } if (shared_params["CPUHybrid"].contains("decompose_subdomains")) { decompose_subdomains = shared_params["CPUHybrid"]["decompose_subdomains"]; @@ -171,7 +280,124 @@ namespace polysolve::linear { additive_mode = shared_params["CPUHybrid"]["additive_mode"]; } + if (shared_params["CPUHybrid"].contains("select_bad_dofs_from_l1_norm")) + { + select_bad_dofs_from_l1_row_norm = shared_params["CPUHybrid"]["select_bad_dofs_from_l1_norm"]; + } + if (shared_params["CPUHybrid"].contains("subdomain_selection_strategy")) + { + const std::string strategy_str = shared_params["CPUHybrid"]["subdomain_selection_strategy"]; + if (strategy_str == "knee") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::KNEE; + } + else if (strategy_str == "gmm") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::GMM; + } + else if (strategy_str == "fd") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::FD; + } + else if (strategy_str == "cost") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::COST; + } + else if (strategy_str == "aposteriori") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::APOSTERIORI; + } + else + { + throw std::runtime_error("Unknown subdomain selection strategy: " + strategy_str); + } + } + if (shared_params["CPUHybrid"].contains("contact_patch_schwarz")) + { + contact_patch_schwarz = shared_params["CPUHybrid"]["contact_patch_schwarz"]; + } + if (shared_params["CPUHybrid"].contains("contact_patch_remove_covered")) + { + contact_patch_remove_covered = shared_params["CPUHybrid"]["contact_patch_remove_covered"]; + } + if (shared_params["CPUHybrid"].contains("ras_correction")) + { + ras_correction = shared_params["CPUHybrid"]["ras_correction"]; + } + if (shared_params["CPUHybrid"].contains("use_strumpack")) + { + use_strumpack = shared_params["CPUHybrid"]["use_strumpack"]; + } + if (shared_params["CPUHybrid"].contains("strumpack_compression")) + { + const std::string compression_str = shared_params["CPUHybrid"]["strumpack_compression"]; + if (compression_str == "none") + { + strumpack_compression = StrumpackCompression::NONE; + } + else if (compression_str == "blr") + { + strumpack_compression = StrumpackCompression::BLR; + } + else if (compression_str == "hss") + { + strumpack_compression = StrumpackCompression::HSS; + } + else + { + throw std::runtime_error("Unknown STRUMPACK compression: " + compression_str); + } + } + if (shared_params["CPUHybrid"].contains("strumpack_compression_rel_tol")) + { + strumpack_compression_rel_tol = shared_params["CPUHybrid"]["strumpack_compression_rel_tol"]; + } + if (shared_params["CPUHybrid"].contains("detailed_log")) + { + detailed_log = shared_params["CPUHybrid"]["detailed_log"]; + } + } + } + + // Set block size for multigrid solvers + void CPUHybridSolver::set_block_size(int block_size) + { + if (myid == 0) + { + SolverCmd cmd = CMD_SET_BLOCK_SIZE; + MPI_Bcast(&cmd, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&solver_id, 1, MPI_INT, 0, MPI_COMM_WORLD); + } + + // Broadcast the block size itself, so worker ranks (which pass in a + // dummy value from run_worker_loop) end up with the root's value. + MPI_Bcast(&block_size, 1, MPI_INT, 0, MPI_COMM_WORLD); + dimension_ = block_size; + } + + // Set the function (block) assigned to each row for multigrid solvers + void CPUHybridSolver::set_block_mapping(const Eigen::VectorXi &block_mapping) + { + if (myid == 0) + { + SolverCmd cmd = CMD_SET_BLOCK_MAPPING; + MPI_Bcast(&cmd, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&solver_id, 1, MPI_INT, 0, MPI_COMM_WORLD); + + block_mapping_ = block_mapping; } + + // Broadcast the full (global) mapping to every rank; each rank keeps the + // whole thing since the local row range ([starts[myid], ends[myid]]) is + // only known once factorize() has partitioned the rows. + int size = myid == 0 ? (int)block_mapping_.size() : 0; + MPI_Bcast(&size, 1, MPI_INT, 0, MPI_COMM_WORLD); + + if (myid != 0) + block_mapping_.resize(size); + + if (size > 0) + MPI_Bcast(block_mapping_.data(), size, MPI_INT, 0, MPI_COMM_WORLD); } void CPUHybridSolver::check_settings() const @@ -180,6 +406,23 @@ namespace polysolve::linear { return; } + + if (subdomain_selection_strategy == SubdomainSelectionStrategy::APOSTERIORI) + { + throw std::runtime_error("A Posteriori subdomain selection strategy is not yet implemented!"); + } + + if (contact_patch_schwarz && ras_correction) + { + throw std::runtime_error("contact_patch_schwarz and ras_correction cannot both be enabled!"); + } + +#ifndef POLYSOLVE_WITH_STRUMPACK + if (use_strumpack) + { + throw std::runtime_error("use_strumpack requires polysolve to be built with POLYSOLVE_WITH_STRUMPACK=ON!"); + } +#endif } void CPUHybridSolver::get_info(json ¶ms) const @@ -204,7 +447,7 @@ namespace polysolve::linear int rows, cols, nnz; if (myid == 0) { - SPDLOG_TRACE("[{}] [start_solve] [0.000000] [num_procs={}] [problem_size={}]", name(), num_procs, Ain.rows()); + CPUHYBRID_LOG_INFO("[{}] [start_solve] [0.000000] [num_procs={}] [problem_size={}]", name(), num_procs, Ain.rows()); rows = Ain.rows(); cols = Ain.cols(); nnz = Ain.nonZeros(); @@ -254,7 +497,7 @@ namespace polysolve::linear MPI_Win_fence(0, A_win); - SPDLOG_TRACE("[{}] [create_shared_matrix_window] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [create_shared_matrix_window] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } SharedSparseMatrix shared_A(rows, cols, nnz, outer_pointers, inner_indices, values); @@ -263,25 +506,49 @@ namespace polysolve::linear bad_indices_sets.clear(); bad_indices_arrays.clear(); + ras_target_dofs.clear(); select_bad_dofs(shared_A); + if (myid == 0 && contact_patch_schwarz) + { + for (const auto &patch : contact_patches) + { + bad_indices_sets.emplace_back(patch.begin(), patch.end()); + bad_indices_arrays.emplace_back(patch.begin(), patch.end()); + } + + if (contact_patch_remove_covered) + { + remove_covered_contact_patches(); + } + } + + if (myid == 0 && ras_correction) + { + // Snapshot all_bad_dofs into per-dof RAS subdomains now, before + // the filter/expand/decompose block below can mutate or replace + // all_bad_dofs -- the same reason contact_patch_schwarz builds + // its subdomains up front instead of going through that block. + build_ras_subdomains(shared_A); + } + if (myid == 0) { - if (decompose_subdomains) + if (decompose_subdomains && !ras_correction) { filter_subdomains(shared_A); } - if (expand_subdomains) + if (expand_subdomains && !ras_correction) { expand_subdomains_to_strongly_connected(shared_A); } - if (decompose_subdomains) + if (decompose_subdomains && !ras_correction) { decompose_subdomains_to_disjoint_subsets(shared_A); } - else + else if (!contact_patch_schwarz && !ras_correction) { bad_indices_sets.emplace_back(all_bad_dofs.begin(), all_bad_dofs.end()); } @@ -290,9 +557,15 @@ namespace polysolve::linear } share_bad_subdomains(); + + if (contact_patch_schwarz) + { + compute_overlap_scale(); + } + factorize_submatrix(shared_A); - SPDLOG_TRACE("[{}] [setup_problematic_dof_precond] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [setup_problematic_dof_precond] [{:.6f}]", name(), elapsed_seconds(phase_begin)); if (has_matrix_) { @@ -305,7 +578,7 @@ namespace polysolve::linear auto phase_begin = clock::now(); copy_matrix_to_hypre(shared_A); has_matrix_ = true; - SPDLOG_TRACE("[{}] [copy_matrix_to_hypre] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [copy_matrix_to_hypre] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } MPI_Win_free(&A_win); @@ -483,6 +756,21 @@ namespace polysolve::linear precond, dimension_, theta); + + if (block_mapping_.size() > 0) + { + assert(block_mapping_.size() == problem_size); + + // HYPRE_IJMatrixCreate for this rank used the global row range + // [starts[myid], ends[myid]], so dof_func must be the matching + // *local* slice of the (globally-indexed) mapping -- each rank + // only reports the function numbers for the rows it owns. + const int n_local = my_size(); + HYPRE_Int *dof_func = hypre_CTAlloc(HYPRE_Int, n_local, HYPRE_MEMORY_HOST); + for (int i = 0; i < n_local; ++i) + dof_func[i] = block_mapping_[starts[myid] + i]; + HYPRE_BoomerAMGSetDofFunc(precond, dof_func); + } } MPI_Barrier(MPI_COMM_WORLD); @@ -496,7 +784,7 @@ namespace polysolve::linear HYPRE_IJVectorGetObject(ij_x, (void **)&par_x); HYPRE_BoomerAMGSetup(precond, parcsr_A, par_b, par_x); - SPDLOG_TRACE("[{}] [amg_setup] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [amg_setup] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } /* Now setup and solve! */ @@ -520,7 +808,7 @@ namespace polysolve::linear matmul(shared_result, buffer); buffer = shared_rhs - buffer; final_res_norm = sqrt(dot(buffer, buffer)); - SPDLOG_TRACE("[{}] [pcg_solve] [{:.6f}] [pcg_iters={}] [residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, final_res_norm); + CPUHYBRID_LOG_INFO("[{}] [pcg_solve] [{:.6f}] [pcg_iters={}] [residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, final_res_norm); } /* Destroy preconditioner */ @@ -539,6 +827,8 @@ namespace polysolve::linear double pre_loop_time; double bi_prod, rel_eps, abs_eps, gamma, old_gamma; + num_iterations = 0; + r.setZero(); { auto phase_begin = clock::now(); @@ -570,7 +860,7 @@ namespace polysolve::linear gamma = dot(r, z1); old_gamma = gamma; - SPDLOG_TRACE("[{}] [pre_loop] [{:.6f}] [rhs_norm={}]", name(), elapsed_seconds(phase_begin), sqrt(bi_prod)); + CPUHYBRID_LOG_INFO("[{}] [pre_loop] [{:.6f}] [rhs_norm={}]", name(), elapsed_seconds(phase_begin), sqrt(bi_prod)); } for (int k = 0; k < max_iter_; ++k) @@ -583,7 +873,7 @@ namespace polysolve::linear if (sdotp == 0.0) { - SPDLOG_TRACE("[{}] [err_zero_sdotp] [0.000000]", name()); + CPUHYBRID_LOG_INFO("[{}] [err_zero_sdotp] [0.000000]", name()); break; } @@ -591,12 +881,12 @@ namespace polysolve::linear if (alpha <= 0.0) { - SPDLOG_TRACE("[{}] [err_negative_alpha] [0.000000]", name()); + CPUHYBRID_LOG_INFO("[{}] [err_negative_alpha] [0.000000]", name()); break; } else if (alpha < __DBL_MIN__) { - SPDLOG_TRACE("[{}] [err_subnormal_alpha] [0.000000]", name()); + CPUHYBRID_LOG_INFO("[{}] [err_subnormal_alpha] [0.000000]", name()); break; } @@ -606,13 +896,13 @@ namespace polysolve::linear double i_prod = dot(r, r); if (rel_eps > 0 && (i_prod / bi_prod) < rel_eps) { - SPDLOG_TRACE("[{}] [converged_rel] [0.000000]", name()); + CPUHYBRID_LOG_INFO("[{}] [converged_rel] [0.000000]", name()); break; } if (abs_eps > 0 && i_prod < abs_eps) { - SPDLOG_TRACE("[{}] [converged_abs] [0.000000]", name()); + CPUHYBRID_LOG_INFO("[{}] [converged_abs] [0.000000]", name()); break; } @@ -625,7 +915,7 @@ namespace polysolve::linear old_gamma = gamma; p = z1 + beta * p; - SPDLOG_TRACE("[{}] [pcg_iter] [{:.6f}] [iter={}] [residual={}]", name(), elapsed_seconds(phase_begin), k, sqrt(i_prod)); + CPUHYBRID_LOG_ITER("[{}] [pcg_iter] [{:.6f}] [iter={}] [residual={}]", name(), elapsed_seconds(phase_begin), k, sqrt(i_prod)); } MPI_Barrier(MPI_COMM_WORLD); } @@ -677,7 +967,7 @@ namespace polysolve::linear HYPRE_BoomerAMGSolve(precond, parcsr_A, par_b, par_x); MPI_Barrier(MPI_COMM_WORLD); - SPDLOG_TRACE("[{}] [amg_v_cycle] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [amg_v_cycle] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } void CPUHybridSolver::dss_precond_iter(Eigen::VectorXd &z, Eigen::VectorXd &r, Eigen::VectorXd &next_z, SharedVector &vec, MPI_Win &vec_win) @@ -694,6 +984,18 @@ namespace polysolve::linear } MPI_Win_fence(0, vec_win); + // Contact patches may overlap, including patches assigned to + // different ranks. Rather than have every rank atomically update a + // shared correction vector, each rank accumulates the subdomains it + // owns into a private buffer (plain += is safe -- nothing else + // touches it) and the per-dof sums are then combined once, across + // all ranks, with a single sum reduction. + Eigen::VectorXd local_correction; + if (contact_patch_schwarz) + { + local_correction = Eigen::VectorXd::Zero(problem_size); + } + int index_counter = 0; for (int index : bad_subdomain_assignments[myid]) { @@ -710,28 +1012,65 @@ namespace polysolve::linear sub_result = D_solvers[index_counter]->solve(sub_rhs); } - for (int i = 0; i < subdomain.size(); ++i) + if (ras_correction) { - vec(2 * problem_size + subdomain[i]) = sub_result(index_mappings[index_counter][subdomain[i]]); + // Restricted Additive Schwarz: this subdomain exists only to + // correct the one dof it was built around, even though the + // solve above covered its whole row neighborhood -- so + // overlapping subdomains never double-correct a shared dof. + const int target = ras_target_dofs[index]; + vec(2 * problem_size + target) = sub_result(index_mappings[index_counter][target]); + } + else + { + for (int i = 0; i < subdomain.size(); ++i) + { + if (contact_patch_schwarz) + { + local_correction(subdomain[i]) += sub_result(index_mappings[index_counter][subdomain[i]]); + } + else + { + vec(2 * problem_size + subdomain[i]) = sub_result(index_mappings[index_counter][subdomain[i]]); + } + } } ++index_counter; } + + if (contact_patch_schwarz) + { + // Sum each rank's private contribution into the true overall + // correction, then damp it by the (global, not per-dof) average + // amount of patch overlap so overlapping patches don't + // over-correct on average. + MPI_Allreduce(MPI_IN_PLACE, local_correction.data(), problem_size, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + local_correction *= overlap_scale_; + } + MPI_Barrier(MPI_COMM_WORLD); MPI_Win_fence(0, vec_win); for (int i = 0; i < my_size(); ++i) { - next_z(i) = vec(starts[myid] + i) + vec(2 * problem_size + starts[myid] + i); + next_z(i) = vec(starts[myid] + i) + (contact_patch_schwarz ? local_correction(starts[myid] + i) : vec(2 * problem_size + starts[myid] + i)); } MPI_Win_fence(0, vec_win); - SPDLOG_TRACE("[{}] [subdomain_solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_ITER("[{}] [subdomain_solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } void CPUHybridSolver::select_bad_dofs(SharedSparseMatrix &sparse_A) { auto phase_begin = clock::now(); + if (!select_bad_dofs_from_l1_row_norm) + { + // AMGF mode: bad DOFs come only from set_problematic_dofs (and, + // if enabled, contact_patch_schwarz), not from this heuristic. + return; + } + MPI_Win row_norm_win; int local_alloc_size = (myid == 0) ? (sparse_A.rows() * sizeof(double)) : 0; void *row_norm_ptr; @@ -775,96 +1114,214 @@ namespace polysolve::linear double min_cost = -1.0; int split_idx = 0; - MPI_Win gamma_win; - void *gamma_ptr; - MPI_Win_allocate_shared(2 * local_alloc_size, 1, MPI_INFO_NULL, MPI_COMM_WORLD, &gamma_ptr, &gamma_win); - if (myid != 0) + switch (subdomain_selection_strategy) { - int disp_unit; - MPI_Aint sz; - MPI_Win_shared_query(gamma_win, 0, &sz, &disp_unit, &gamma_ptr); - } + case SubdomainSelectionStrategy::KNEE: + case SubdomainSelectionStrategy::FD: + case SubdomainSelectionStrategy::COST: + { + if (myid == 0) + { + std::vector> sorted_norms; + sorted_norms.reserve(row_norms.size()); + for (int i = 0; i < row_norms.size(); ++i) + { + sorted_norms.push_back({row_norms(i), i}); + } + + std::sort(sorted_norms.begin(), sorted_norms.end()); + int N = sorted_norms.size(); + + if (N >= 2) + { + if (subdomain_selection_strategy == SubdomainSelectionStrategy::KNEE) + { + double x1 = 0.0; + double y1 = std::log(sorted_norms[0].first + 1e-12); + double x2 = N - 1.0; + double y2 = std::log(sorted_norms[N - 1].first + 1e-12); + + for (int i = 1; i < N - 1; ++i) + { + double x0 = (double)i; + double y0 = std::log(sorted_norms[i].first + 1e-12); + + double dist = std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1); + + if (dist > max_dist) + { + max_dist = dist; + split_idx = i; + } + } + } + else if (subdomain_selection_strategy == SubdomainSelectionStrategy::FD) + { + for (int i = 1; i < N; ++i) + { + double jump = sorted_norms[i].first - sorted_norms[i - 1].first; + if (jump > max_jump) + { + max_jump = jump; + split_idx = i; + } + } + } + else if (subdomain_selection_strategy == SubdomainSelectionStrategy::COST) + { + // TODO: Migrate these to class variables later + const double c_pcg_conv = 1.0; + const double c_pcg = 1.0; + const double c_fac = 1.0; + const double c_solve = 1.0; - SharedVector gamma((double *)gamma_ptr, 2 * sparse_A.rows()); + min_cost = std::numeric_limits::max(); - mean_0 = row_norms.minCoeff(); - var_0 = global_var; - mean_1 = row_norms.maxCoeff(); - var_1 = global_var; - double w0 = 0.5; - double w1 = 0.5; - double var_reg = 1e-6; + double lambda_min = sorted_norms[0].first; - for (gmm_iter = 0; gmm_iter < max_gmm_iterations; ++gmm_iter) + for (int i = 0; i < N; ++i) + { + double p = static_cast(N - i); + double lambda_max = sorted_norms[i].first; + + double m = c_pcg_conv * std::sqrt(lambda_max / lambda_min); + double cost = (m * c_pcg * N) + + (c_fac * p * p) + + (m * c_solve * std::pow(p, 4.0 / 3.0)); + + if (cost < min_cost) + { + min_cost = cost; + split_idx = i; + } + } + } + + for (int i = split_idx; i < N; ++i) + { + all_bad_dofs.insert(sorted_norms[i].second); + } + } + } + break; + } + case SubdomainSelectionStrategy::GMM: { - double log_w0 = std::log(w0); - double log_w1 = std::log(w1); - double log_norm_const_0 = -0.5 * std::log(2.0 * M_PI * var_0); - double log_norm_const_1 = -0.5 * std::log(2.0 * M_PI * var_1); + MPI_Win gamma_win; + void *gamma_ptr; + MPI_Win_allocate_shared(2 * local_alloc_size, 1, MPI_INFO_NULL, MPI_COMM_WORLD, &gamma_ptr, &gamma_win); + if (myid != 0) + { + int disp_unit; + MPI_Aint sz; + MPI_Win_shared_query(gamma_win, 0, &sz, &disp_unit, &gamma_ptr); + } - auto x = row_norms.segment(starts[myid], my_size()).array(); + SharedVector gamma((double *)gamma_ptr, 2 * sparse_A.rows()); - Eigen::ArrayXd log_g0 = log_w0 + log_norm_const_0 - 0.5 * (x - mean_0).square() / var_0; - Eigen::ArrayXd log_g1 = log_w1 + log_norm_const_1 - 0.5 * (x - mean_1).square() / var_1; - Eigen::ArrayXd max_log_g = log_g0.cwiseMax(log_g1); - Eigen::ArrayXd log_total = max_log_g + ((log_g0 - max_log_g).exp() + (log_g1 - max_log_g).exp()).log(); + mean_0 = row_norms.minCoeff(); + var_0 = global_var; + mean_1 = row_norms.maxCoeff(); + var_1 = global_var; + double w0 = 0.5; + double w1 = 0.5; + double var_reg = 1e-6; - MPI_Win_fence(0, gamma_win); - gamma.segment(starts[myid], my_size()).array() = (log_g0 - log_total).exp(); - gamma.segment(starts[myid] + row_norms.size(), my_size()).array() = (log_g1 - log_total).exp(); - MPI_Win_fence(0, gamma_win); + for (gmm_iter = 0; gmm_iter < max_gmm_iterations; ++gmm_iter) + { + double log_w0 = std::log(w0); + double log_w1 = std::log(w1); + double log_norm_const_0 = -0.5 * std::log(2.0 * M_PI * var_0); + double log_norm_const_1 = -0.5 * std::log(2.0 * M_PI * var_1); - w0 = 1.0 / row_norms.size() * gamma.segment(starts[myid], my_size()).sum(); - w1 = 1.0 / row_norms.size() * gamma.segment(starts[myid] + row_norms.size(), my_size()).sum(); + auto x = row_norms.segment(starts[myid], my_size()).array(); - MPI_Allreduce(MPI_IN_PLACE, &w0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &w1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + Eigen::ArrayXd log_g0 = log_w0 + log_norm_const_0 - 0.5 * (x - mean_0).square() / var_0; + Eigen::ArrayXd log_g1 = log_w1 + log_norm_const_1 - 0.5 * (x - mean_1).square() / var_1; + Eigen::ArrayXd max_log_g = log_g0.cwiseMax(log_g1); + Eigen::ArrayXd log_total = max_log_g + ((log_g0 - max_log_g).exp() + (log_g1 - max_log_g).exp()).log(); - double old_mean_0 = mean_0; - double old_mean_1 = mean_1; - double old_var_0 = var_0; - double old_var_1 = var_1; + MPI_Win_fence(0, gamma_win); + gamma.segment(starts[myid], my_size()).array() = (log_g0 - log_total).exp(); + gamma.segment(starts[myid] + row_norms.size(), my_size()).array() = (log_g1 - log_total).exp(); + MPI_Win_fence(0, gamma_win); - mean_0 = (row_norms.segment(starts[myid], my_size()).array() * gamma.segment(starts[myid], my_size()).array()).sum() / (w0 * row_norms.size()); - mean_1 = (row_norms.segment(starts[myid], my_size()).array() * gamma.segment(starts[myid] + row_norms.size(), my_size()).array()).sum() / (w1 * row_norms.size()); + w0 = 1.0 / row_norms.size() * gamma.segment(starts[myid], my_size()).sum(); + w1 = 1.0 / row_norms.size() * gamma.segment(starts[myid] + row_norms.size(), my_size()).sum(); - MPI_Allreduce(MPI_IN_PLACE, &mean_0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &mean_1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &w0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &w1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - var_0 = (gamma.segment(starts[myid], my_size()).array() * (row_norms.segment(starts[myid], my_size()).array() - mean_0).square()).sum() / (w0 * row_norms.size()); - var_1 = (gamma.segment(starts[myid] + row_norms.size(), my_size()).array() * (row_norms.segment(starts[myid], my_size()).array() - mean_1).square()).sum() / (w1 * row_norms.size()); + double old_mean_0 = mean_0; + double old_mean_1 = mean_1; + double old_var_0 = var_0; + double old_var_1 = var_1; - MPI_Allreduce(MPI_IN_PLACE, &var_0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &var_1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + mean_0 = (row_norms.segment(starts[myid], my_size()).array() * gamma.segment(starts[myid], my_size()).array()).sum() / (w0 * row_norms.size()); + mean_1 = (row_norms.segment(starts[myid], my_size()).array() * gamma.segment(starts[myid] + row_norms.size(), my_size()).array()).sum() / (w1 * row_norms.size()); - var_0 += var_reg; - var_1 += var_reg; + MPI_Allreduce(MPI_IN_PLACE, &mean_0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &mean_1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - if (std::abs(mean_0 - old_mean_0) / std::abs(old_mean_0) < gmm_tol && std::abs(mean_1 - old_mean_1) / std::abs(old_mean_1) < gmm_tol && std::abs(var_0 - old_var_0) / std::abs(old_var_0) < gmm_tol && std::abs(var_1 - old_var_1) / std::abs(old_var_1) < gmm_tol) - { - break; + var_0 = (gamma.segment(starts[myid], my_size()).array() * (row_norms.segment(starts[myid], my_size()).array() - mean_0).square()).sum() / (w0 * row_norms.size()); + var_1 = (gamma.segment(starts[myid] + row_norms.size(), my_size()).array() * (row_norms.segment(starts[myid], my_size()).array() - mean_1).square()).sum() / (w1 * row_norms.size()); + + MPI_Allreduce(MPI_IN_PLACE, &var_0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &var_1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + var_0 += var_reg; + var_1 += var_reg; + + if (std::abs(mean_0 - old_mean_0) / std::abs(old_mean_0) < gmm_tol && std::abs(mean_1 - old_mean_1) / std::abs(old_mean_1) < gmm_tol && std::abs(var_0 - old_var_0) / std::abs(old_var_0) < gmm_tol && std::abs(var_1 - old_var_1) / std::abs(old_var_1) < gmm_tol) + { + break; + } } - } - if (myid == 0) - { - if (std::abs(mean_1) / std::abs(mean_0) > gmm_jump_threshold) + if (myid == 0) { - for (int i = 0; i < row_norms.size(); ++i) + if (std::abs(mean_1) / std::abs(mean_0) > gmm_jump_threshold) { - if (gamma(i) < gamma(i + row_norms.size())) + for (int i = 0; i < row_norms.size(); ++i) { - all_bad_dofs.insert(i); + if (gamma(i) < gamma(i + row_norms.size())) + { + all_bad_dofs.insert(i); + } } } } + + MPI_Win_free(&gamma_win); + break; + } + case SubdomainSelectionStrategy::APOSTERIORI: + break; } - MPI_Win_free(&gamma_win); MPI_Win_free(&row_norm_win); - SPDLOG_TRACE("[{}] [bad_dof_selection] [{}] [strategy=GMM] [global_mean={}] [global_var={}] [mean_0={}] [mean_1={}] [var_0={}] [var_1={}] [gmm_iters={}] [num_bad_dofs={}]", - name(), elapsed_seconds(phase_begin), global_mean, global_var, mean_0, mean_1, var_0, var_1, gmm_iter, all_bad_dofs.size()); + switch (subdomain_selection_strategy) + { + case SubdomainSelectionStrategy::KNEE: + CPUHYBRID_LOG_INFO("[{}] [bad_dof_selection] [{}] [strategy=KNEE] [global_mean={}] [global_var={}] [split_idx={}] [max_dist={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, max_dist, all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::FD: + CPUHYBRID_LOG_INFO("[{}] [bad_dof_selection] [{}] [strategy=FD] [global_mean={}] [global_var={}] [split_idx={}] [max_jump={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, max_jump, all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::COST: + CPUHYBRID_LOG_INFO("[{}] [bad_dof_selection] [{}] [strategy=COST] [global_mean={}] [global_var={}] [split_idx={}] [min_cost={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, min_cost, all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::GMM: + CPUHYBRID_LOG_INFO("[{}] [bad_dof_selection] [{}] [strategy=GMM] [global_mean={}] [global_var={}] [mean_0={}] [mean_1={}] [var_0={}] [var_1={}] [gmm_iters={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, mean_0, mean_1, var_0, var_1, gmm_iter, all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::APOSTERIORI: + break; + } } void CPUHybridSolver::factorize_submatrix(SharedSparseMatrix &sparse_A) @@ -876,6 +1333,13 @@ namespace polysolve::linear { if (bad_indices_sets[i].size() > 1000) { +#ifdef POLYSOLVE_WITH_STRUMPACK + if (use_strumpack) + { + D_solvers.push_back(std::make_unique(strumpack_compression, strumpack_compression_rel_tol)); + continue; + } +#endif #if POLYSOLVE_WITH_MKL D_solvers.push_back(std::make_unique>>>()); #elif POLYSOLVE_WITH_ACCELERATE @@ -905,7 +1369,7 @@ namespace polysolve::linear } MPI_Barrier(MPI_COMM_WORLD); - SPDLOG_TRACE("[{}] [factorize_submatrix] [{}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_INFO("[{}] [factorize_submatrix] [{}]", name(), elapsed_seconds(phase_begin)); } void CPUHybridSolver::matmul(Eigen::VectorXd &x, Eigen::VectorXd &result) @@ -922,7 +1386,7 @@ namespace polysolve::linear HYPRE_IJVectorAssemble(ij_b); HYPRE_IJVectorGetObject(ij_b, (void **)&par_result); HYPRE_ParCSRMatrixMatvec(1.0, parcsr_A, par_x, 0.0, par_result); - SPDLOG_TRACE("[{}] [matmul] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + CPUHYBRID_LOG_ITER("[{}] [matmul] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } double CPUHybridSolver::dot(Eigen::VectorXd &a, Eigen::VectorXd &b) @@ -1177,7 +1641,7 @@ namespace polysolve::linear ++num_not_poorly_conditioned; } - SPDLOG_TRACE("[{}] [subdomain_filtering] [{}] [total_dofs_before={}] [total_dofs_after={}] [num_too_small={}] [num_too_large={}] [num_not_poorly_conditioned={}]", + CPUHYBRID_LOG_INFO("[{}] [subdomain_filtering] [{}] [total_dofs_before={}] [total_dofs_after={}] [num_too_small={}] [num_too_large={}] [num_not_poorly_conditioned={}]", name(), elapsed_seconds(phase_begin), original_num_bad_dofs, all_bad_dofs.size(), num_too_small, num_too_large, num_not_poorly_conditioned); } @@ -1187,7 +1651,6 @@ namespace polysolve::linear int num_bad_dofs_before = all_bad_dofs.size(); std::set new_bad_dofs; - ; for (int k : all_bad_dofs) { @@ -1198,10 +1661,109 @@ namespace polysolve::linear } all_bad_dofs = std::move(new_bad_dofs); - SPDLOG_TRACE("[{}] [subdomain_expansion] [{}] [num_dofs_before={}] [num_dofs_after={}]", + CPUHYBRID_LOG_INFO("[{}] [subdomain_expansion] [{}] [num_dofs_before={}] [num_dofs_after={}]", name(), elapsed_seconds(phase_begin), num_bad_dofs_before, all_bad_dofs.size()); } + // A patch fully contained in another contributes nothing that the larger + // patch doesn't already correct, so it is dropped rather than solved and + // added redundantly into the overlap sum in dss_precond_iter. Patches are + // processed smallest-first so each is only ever tested against patches + // that could actually contain it; ties (identical patches) keep exactly + // one copy rather than deleting each other. + void CPUHybridSolver::remove_covered_contact_patches() + { + auto phase_begin = clock::now(); + const int num_before = static_cast(bad_indices_sets.size()); + + std::vector order(bad_indices_sets.size()); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](int a, int b) { + if (bad_indices_sets[a].size() != bad_indices_sets[b].size()) + return bad_indices_sets[a].size() < bad_indices_sets[b].size(); + return a < b; + }); + + std::vector covered(bad_indices_sets.size(), false); + for (size_t oi = 0; oi < order.size(); ++oi) + { + const int i = order[oi]; + if (covered[i]) + continue; + + for (size_t oj = oi + 1; oj < order.size(); ++oj) + { + const int j = order[oj]; + if (covered[j] || bad_indices_sets[i].size() > bad_indices_sets[j].size()) + continue; + + if (std::includes(bad_indices_sets[j].begin(), bad_indices_sets[j].end(), + bad_indices_sets[i].begin(), bad_indices_sets[i].end())) + { + if (bad_indices_sets[i].size() == bad_indices_sets[j].size()) + { + // Identical patches: keep i (processed first), drop j. + covered[j] = true; + } + else + { + covered[i] = true; + break; + } + } + } + } + + std::vector> kept_sets; + std::vector> kept_arrays; + kept_sets.reserve(bad_indices_sets.size()); + kept_arrays.reserve(bad_indices_arrays.size()); + for (int i = 0; i < static_cast(bad_indices_sets.size()); ++i) + { + if (!covered[i]) + { + kept_sets.push_back(std::move(bad_indices_sets[i])); + kept_arrays.push_back(std::move(bad_indices_arrays[i])); + } + } + bad_indices_sets = std::move(kept_sets); + bad_indices_arrays = std::move(kept_arrays); + + CPUHYBRID_LOG_INFO("[{}] [remove_covered_contact_patches] [{:.6f}] [num_patches_before={}] [num_patches_after={}] [num_removed={}]", + name(), elapsed_seconds(phase_begin), num_before, bad_indices_sets.size(), num_before - static_cast(bad_indices_sets.size())); + } + + // RAS: one subdomain per bad dof, made of every dof with a nonzero entry + // in that dof's row. sparse_A is column-major, so column k is used as a + // stand-in for row k -- the same structural-symmetry assumption already + // relied on elsewhere in this file (e.g. assemble_D's global_to_row). + // Unlike the other subdomain-construction paths, the subdomain a dof + // lands in here is never merged, filtered, or expanded: it is exactly the + // dof's own sparsity neighborhood, and dss_precond_iter restricts its + // correction to that one owning dof (ras_target_dofs), so overlapping + // subdomains never double-correct a dof the way additive Schwarz would. + void CPUHybridSolver::build_ras_subdomains(SharedSparseMatrix &sparse_A) + { + auto phase_begin = clock::now(); + + for (int d : all_bad_dofs) + { + std::set subdomain; + for (SharedSparseMatrix::InnerIterator it(sparse_A, d); it; ++it) + { + subdomain.insert(it.row()); + } + subdomain.insert(d); + + bad_indices_sets.push_back(subdomain); + bad_indices_arrays.emplace_back(subdomain.begin(), subdomain.end()); + ras_target_dofs.push_back(d); + } + + CPUHYBRID_LOG_INFO("[{}] [build_ras_subdomains] [{:.6f}] [num_bad_dofs={}] [num_subdomains={}]", + name(), elapsed_seconds(phase_begin), all_bad_dofs.size(), bad_indices_sets.size()); + } + void CPUHybridSolver::decompose_subdomains_to_disjoint_subsets(SharedSparseMatrix &sparse_A) { auto phase_begin = clock::now(); @@ -1242,7 +1804,7 @@ namespace polysolve::linear } bad_indices_sets.emplace_back(kv.second.begin(), kv.second.end()); } - SPDLOG_TRACE("[{}] [subdomain_decomposition] [{}] [num_subdomains={}] ", + CPUHYBRID_LOG_INFO("[{}] [subdomain_decomposition] [{}] [num_subdomains={}] ", name(), elapsed_seconds(phase_begin), bad_indices_sets.size()); } @@ -1304,6 +1866,20 @@ namespace polysolve::linear } } + if (ras_correction) + { + int ras_size = (myid == 0) ? static_cast(ras_target_dofs.size()) : 0; + MPI_Bcast(&ras_size, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (myid != 0) + { + ras_target_dofs.resize(ras_size); + } + if (ras_size > 0) + { + MPI_Bcast(ras_target_dofs.data(), ras_size, MPI_INT, 0, MPI_COMM_WORLD); + } + } + bad_indices_arrays.clear(); bad_indices_arrays.resize(bad_indices_sets.size()); for (int i = 0; i < bad_indices_sets.size(); ++i) @@ -1314,10 +1890,46 @@ namespace polysolve::linear bad_indices_arrays[i].push_back(index); } } - SPDLOG_TRACE("[{}] [share_bad_subdomains] [{}] ", + CPUHYBRID_LOG_INFO("[{}] [share_bad_subdomains] [{}] ", name(), elapsed_seconds(phase_begin)); } + // Under contact_patch_schwarz, bad_indices_sets (the contact patches) may + // overlap, so a covered dof can belong to more than one patch. Compute a + // single global scale factor -- 1 / (average number of patches per + // covered dof) -- used to damp the summed correction in dss_precond_iter + // so overlapping patches don't over-correct on average. bad_indices_sets + // is identical on every rank by this point (share_bad_subdomains already + // broadcast it), so no communication is needed here. + void CPUHybridSolver::compute_overlap_scale() + { + std::unordered_map dof_patch_count; + for (const auto &subdomain : bad_indices_sets) + { + for (int dof : subdomain) + { + ++dof_patch_count[dof]; + } + } + + if (dof_patch_count.empty()) + { + overlap_scale_ = 1.0; + return; + } + + long long total_incidences = 0; + for (const auto &kv : dof_patch_count) + { + total_incidences += kv.second; + } + + const double avg_overlap = static_cast(total_incidences) / static_cast(dof_patch_count.size()); + overlap_scale_ = 1.0 / avg_overlap; + + CPUHYBRID_LOG_INFO("[{}] [contact_patch_overlap] [avg_overlap={}] [overlap_scale={}]", name(), avg_overlap, overlap_scale_); + } + void CPUHybridSolver::load_balance_subdomains() { auto phase_begin = clock::now(); @@ -1356,7 +1968,7 @@ namespace polysolve::linear const int max_size = subdomain_sizes.size() > 0 ? subdomain_sizes.front().second : 0; const int min_size = subdomain_sizes.size() > 0 ? subdomain_sizes.back().second : 0; - SPDLOG_TRACE("[{}] [subdomain_load_balance] [{}] [max_size={}] [min_size={}] [total_dofs={}]", + CPUHYBRID_LOG_INFO("[{}] [subdomain_load_balance] [{}] [max_size={}] [min_size={}] [total_dofs={}]", name(), elapsed_seconds(phase_begin), max_size, min_size, total_bad_dofs); } @@ -1375,7 +1987,6 @@ namespace polysolve::linear Eigen::setNbThreads(1); HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST); HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST); - spdlog::set_level(spdlog::level::off); CPUHybridSolver::run_worker_loop(); return 0; } @@ -1450,6 +2061,18 @@ namespace polysolve::linear worker_registry[id]->set_parameters(dummy_params); break; } + case CMD_SET_BLOCK_SIZE: + { + int dummy_block_size = 0; + worker_registry[id]->set_block_size(dummy_block_size); + break; + } + case CMD_SET_BLOCK_MAPPING: + { + Eigen::VectorXi dummy_block_mapping; + worker_registry[id]->set_block_mapping(dummy_block_mapping); + break; + } case CMD_FACTORIZE: { StiffnessMatrix dummy_A; diff --git a/src/polysolve/linear/CPUHybridSolver.hpp b/src/polysolve/linear/CPUHybridSolver.hpp index ddc13d8a..52ca36c0 100644 --- a/src/polysolve/linear/CPUHybridSolver.hpp +++ b/src/polysolve/linear/CPUHybridSolver.hpp @@ -66,6 +66,8 @@ namespace polysolve::linear { CMD_CREATE, CMD_SET_PARAMETERS, + CMD_SET_BLOCK_SIZE, + CMD_SET_BLOCK_MAPPING, CMD_FACTORIZE, CMD_SOLVE, CMD_DESTROY, @@ -104,6 +106,12 @@ namespace polysolve::linear // Set solver parameters virtual void set_parameters(const json ¶ms) override; + // Set block size for multigrid solvers + virtual void set_block_size(int block_size) override; + + // Set the function (block) assigned to each row for multigrid solvers + virtual void set_block_mapping(const Eigen::VectorXi &block_mapping) override; + // Retrieve solve information virtual void get_info(json ¶ms) const override; @@ -118,14 +126,41 @@ namespace polysolve::linear // Name of the solver type (for debugging purposes) virtual std::string name() const override { - return "CPUHybrid"; + return select_bad_dofs_from_l1_row_norm ? "CPUHybrid" : "CPUAMGF"; } + virtual void set_problematic_dofs(const std::set &bad_dofs) override { all_bad_dofs = bad_dofs; } + + // Alternative strategies for picking which DOFs get the dense/direct + // subdomain treatment, in place of the default GMM row-norm clustering. + enum class SubdomainSelectionStrategy + { + KNEE, + GMM, + FD, + COST, + APOSTERIORI + }; + + // Low-rank compression of STRUMPACK's frontal matrices, under + // use_strumpack. NONE is an exact LU factorization. + enum class StrumpackCompression + { + NONE, + BLR, + HSS + }; + protected: // AMG settings double theta = 0.5; // Hybrid preconditioner settings + // When false, bad DOFs are taken only from whatever set_problematic_dofs + // (and, if enabled, contact_patch_schwarz) supplies instead of being + // inferred from the row-norm heuristic below (AMGF mode). + bool select_bad_dofs_from_l1_row_norm = true; + SubdomainSelectionStrategy subdomain_selection_strategy = SubdomainSelectionStrategy::GMM; bool decompose_subdomains = true; int min_subdomain_size = 3; int max_subdomain_size = 1e9; @@ -134,17 +169,52 @@ namespace polysolve::linear int max_gmm_iterations = 20; bool expand_subdomains = true; bool additive_mode = false; + // Give each contact patch its own (possibly overlapping) subdomain, + // solved additively -- an additive Schwarz treatment of contact DOFs. + bool contact_patch_schwarz = false; + // Under contact_patch_schwarz, drop any patch that is fully contained + // in another patch before factorizing -- it would only ever be + // redundant additive correction, since the larger patch already + // covers every one of its dofs. + bool contact_patch_remove_covered = false; + // Restricted Additive Schwarz: build one subdomain per bad dof (every + // dof with a nonzero entry in that dof's row/column, since the matrix + // is assumed structurally symmetric -- see assemble_D), but only let + // that subdomain correct the single dof it was built around. Uses the + // same bad-dof set as AMGF/select_bad_dofs. Mutually exclusive with + // contact_patch_schwarz. + bool ras_correction = false; + // Experimental: factorize the large subdomains -- the ones that would + // otherwise go to Pardiso/Accelerate -- with STRUMPACK's multifrontal + // LU instead, optionally with low-rank compression of its fronts. + // Needs polysolve built with POLYSOLVE_WITH_STRUMPACK. + bool use_strumpack = false; + StrumpackCompression strumpack_compression = StrumpackCompression::NONE; + double strumpack_compression_rel_tol = 1e-4; // General solver settings int dimension_ = 1; // 1 = scalar (Laplace), 2 or 3 = vector (Elasticity) + + // Optional per-row function (block) assignment for HYPRE_BoomerAMGSetDofFunc, + // indexed globally (i.e. by the row indices factorize() is given). Empty means + // the default interleaved mapping (row i belongs to function i % dimension_). + // Kept identically on every rank; each rank slices out its own local rows + // ([starts[myid], ends[myid]]) when handing HYPRE its local dof_func array. + Eigen::VectorXi block_mapping_; + int max_iter_ = 10000; double rel_conv_tol_ = 1e-10; double abs_conv_tol_ = 0.0; double conditioning_threshold = 100.0; + // When false, skips the per-PCG-iteration timing logs (matmul, + // amg_v_cycle, subdomain_solve, pcg_iter) so a long solve doesn't + // flood the log; every other log line is unaffected. + bool detailed_log = true; + // solve information - HYPRE_Int num_iterations; - HYPRE_Complex final_res_norm; + HYPRE_Int num_iterations = 0; + HYPRE_Complex final_res_norm = 0; private: bool has_matrix_ = false; @@ -172,6 +242,21 @@ namespace polysolve::linear std::vector> bad_subdomain_assignments; std::vector> index_mappings; + // Under ras_correction, bad_indices_sets/bad_indices_arrays[i] holds + // the RAS subdomain built around a single dof; ras_target_dofs[i] is + // that dof -- the only one whose correction the subdomain is allowed + // to write back. Shared to every rank in share_bad_subdomains, same + // as bad_indices_sets. + std::vector ras_target_dofs; + + // Under contact_patch_schwarz, contact patches (bad_indices_sets) may + // overlap, so a dof's additive-Schwarz correction is the sum of every + // patch that covers it. That sum is damped by this single global + // factor -- 1 / (average number of patches covering a covered dof) -- + // rather than dividing each dof by its own multiplicity, so overlap + // is corrected for on average without per-dof bookkeeping. + double overlap_scale_ = 1.0; + Eigen::VectorXd shared_rhs, shared_result; Eigen::VectorXd z1, z2, z3; Eigen::VectorXd r, p, buffer; @@ -190,10 +275,13 @@ namespace polysolve::linear void decompose_subdomains_to_disjoint_subsets(SharedSparseMatrix &sparse_A); void filter_subdomains(SharedSparseMatrix &sparse_A); void expand_subdomains_to_strongly_connected(SharedSparseMatrix &sparse_A); + void remove_covered_contact_patches(); + void build_ras_subdomains(SharedSparseMatrix &sparse_A); void share_bad_subdomains(); void load_balance_subdomains(); void select_bad_dofs(SharedSparseMatrix &sparse_A); void factorize_submatrix(SharedSparseMatrix &sparse_A); + void compute_overlap_scale(); // matrix multiplication void matmul(Eigen::VectorXd &x, Eigen::VectorXd &result); diff --git a/src/polysolve/linear/EigenSolver.tpp b/src/polysolve/linear/EigenSolver.tpp index 4e8d3852..00f7c653 100644 --- a/src/polysolve/linear/EigenSolver.tpp +++ b/src/polysolve/linear/EigenSolver.tpp @@ -2,7 +2,9 @@ //////////////////////////////////////////////////////////////////////////////// #include "EigenSolver.hpp" +#include #include +#include //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -10,6 +12,16 @@ //////////////////////////////////////////////////////////////////////////////// namespace polysolve::linear { + namespace + { + using clock = std::chrono::steady_clock; + + double elapsed_seconds(const std::chrono::time_point &begin) + { + return std::chrono::duration(clock::now() - begin).count(); + } + } // namespace + // Get info on the last solve step template void EigenDirect::get_info(json ¶ms) const @@ -37,18 +49,22 @@ namespace polysolve::linear template void EigenDirect::analyze_pattern(const StiffnessMatrix &A, const int precond_num) { + auto phase_begin = clock::now(); m_Solver.analyzePattern(A); + SPDLOG_INFO("[{}] [analyze_pattern] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } // Factorize system matrix template void EigenDirect::factorize(const StiffnessMatrix &A) { + auto phase_begin = clock::now(); m_Solver.factorize(A); if (m_Solver.info() == Eigen::NumericalIssue) { throw std::runtime_error("[EigenDirect] NumericalIssue encountered."); } + SPDLOG_INFO("[{}] [factorize] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } // Solve the linear system @@ -56,7 +72,9 @@ namespace polysolve::linear void EigenDirect::solve( const Ref b, Ref x) { + auto phase_begin = clock::now(); x = m_Solver.solve(b); + SPDLOG_INFO("[{}] [solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } //////////////////////////////////////////////////////////////////////////////// @@ -93,15 +111,19 @@ namespace polysolve::linear template void EigenIterative::analyze_pattern(const StiffnessMatrix &A, const int precond_num) { + auto phase_begin = clock::now(); m_Solver.analyzePattern(A); + SPDLOG_INFO("[{}] [analyze_pattern] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } // Factorize system matrix template void EigenIterative::factorize(const StiffnessMatrix &A) { + auto phase_begin = clock::now(); m_A = A; m_Solver.factorize(m_A); + SPDLOG_INFO("[{}] [factorize] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } // Solve the linear system @@ -110,7 +132,9 @@ namespace polysolve::linear const Ref b, Ref x) { assert(x.size() == b.size()); + auto phase_begin = clock::now(); x = m_Solver.solveWithGuess(b, x); + SPDLOG_INFO("[{}] [solve] [{:.6f}] [iters={}]", name(), elapsed_seconds(phase_begin), m_Solver.iterations()); } //////////////////////////////////////////////////////////////////////////////// @@ -134,7 +158,9 @@ namespace polysolve::linear template void EigenDenseSolver::factorize_dense(const Eigen::MatrixXd &A) { + auto phase_begin = clock::now(); m_Solver.compute(A); + SPDLOG_INFO("[{}] [factorize] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } // Solve the linear system @@ -142,6 +168,8 @@ namespace polysolve::linear void EigenDenseSolver::solve( const Ref b, Ref x) { + auto phase_begin = clock::now(); x = m_Solver.solve(b); + SPDLOG_INFO("[{}] [solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } } // namespace polysolve::linear diff --git a/src/polysolve/linear/GPUHybridSolver.cu b/src/polysolve/linear/GPUHybridSolver.cu index e34ded99..245d1b80 100644 --- a/src/polysolve/linear/GPUHybridSolver.cu +++ b/src/polysolve/linear/GPUHybridSolver.cu @@ -21,15 +21,25 @@ #include #include #include +#include +#include +#include #include #include +#include +#include #include #include #include +// Internal header, needed for hypre_CTAlloc/hypre_TFree/hypre_TMemcpy so the +// dof_func array handed to HYPRE_BoomerAMGSetDofFunc uses the allocator HYPRE +// itself frees it with (and lives in the device memory HYPRE expects here). +#include "_hypre_utilities.h" + #ifdef HYPRE_ENABLE_MPI #include #endif @@ -59,6 +69,29 @@ } \ } while (0) +#define CHECK_CUBLAS(call) \ + do \ + { \ + cublasStatus_t status = call; \ + if (status != CUBLAS_STATUS_SUCCESS) \ + { \ + std::cerr << "cuBLAS Error at " << __FILE__ << ":" << __LINE__ \ + << " code " << (int)status << std::endl; \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + + // Per-PCG-iteration timing logs (matmul, amg_v_cycle, subdomain_solve, + // contact_patch_dense_subdomain_solve, pcg_iter) go through this instead + // of SPDLOG_INFO directly, so they can be silenced with + // detailed_log=false without affecting any other log line. +#define GPUHYBRID_LOG_ITER(...) \ + do \ + { \ + if (detailed_log) \ + SPDLOG_INFO(__VA_ARGS__); \ + } while (0) + namespace polysolve::linear { @@ -94,6 +127,7 @@ namespace polysolve::linear HYPRE_SetUseGpuRand(true); CHECK_CUDSS(cudssCreate(&cudss_handle)); + CHECK_CUBLAS(cublasCreate(&cublas_handle)); } void GPUHybridSolver::set_parameters(const json ¶ms) @@ -116,10 +150,6 @@ namespace polysolve::linear { theta = params["GPUHybrid"]["theta"]; } - if (params["GPUHybrid"].contains("block_dim")) - { - dimension_ = params["GPUHybrid"]["block_dim"]; - } if (params["GPUHybrid"].contains("decompose_subdomains")) { decompose_subdomains = params["GPUHybrid"]["decompose_subdomains"]; @@ -156,9 +186,59 @@ namespace polysolve::linear { additive_mode = params["GPUHybrid"]["additive_mode"]; } + if (params["GPUHybrid"].contains("contact_patch_schwarz")) + { + contact_patch_schwarz = params["GPUHybrid"]["contact_patch_schwarz"]; + } + if (params["GPUHybrid"].contains("detailed_log")) + { + detailed_log = params["GPUHybrid"]["detailed_log"]; + } + if (params["GPUHybrid"].contains("select_bad_dofs_from_l1_norm")) + { + select_bad_dofs_from_l1_row_norm = params["GPUHybrid"]["select_bad_dofs_from_l1_norm"]; + } + if (params["GPUHybrid"].contains("subdomain_selection_strategy")) + { + const std::string strategy_str = params["GPUHybrid"]["subdomain_selection_strategy"]; + if (strategy_str == "knee") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::KNEE; + } + else if (strategy_str == "gmm") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::GMM; + } + else if (strategy_str == "fd") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::FD; + } + else if (strategy_str == "cost") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::COST; + } + else if (strategy_str == "aposteriori") + { + subdomain_selection_strategy = SubdomainSelectionStrategy::APOSTERIORI; + } + else + { + throw std::runtime_error("Unknown subdomain selection strategy: " + strategy_str); + } + } } } + void GPUHybridSolver::set_block_size(int block_size) + { + dimension_ = block_size; + } + + void GPUHybridSolver::set_block_mapping(const Eigen::VectorXi &block_mapping) + { + block_mapping_ = block_mapping; + } + void GPUHybridSolver::get_info(json ¶ms) const { params["num_iterations"] = num_iterations; @@ -167,12 +247,16 @@ namespace polysolve::linear void GPUHybridSolver::check_settings() const { + if (subdomain_selection_strategy == SubdomainSelectionStrategy::APOSTERIORI) + { + throw std::runtime_error("A Posteriori subdomain selection strategy is not yet implemented!"); + } } void GPUHybridSolver::factorize(const StiffnessMatrix &Ain) { check_settings(); - SPDLOG_TRACE("[{}] [start_solve] [0.000000] [problem_size={}]", name(), Ain.rows()); + SPDLOG_INFO("[{}] [start_solve] [0.000000] [problem_size={}]", name(), Ain.rows()); { auto phase_begin = clock::now(); @@ -185,7 +269,7 @@ namespace polysolve::linear thrust::copy(Ain.innerIndexPtr(), Ain.innerIndexPtr() + d_inner_indices.size(), d_inner_indices.begin()); thrust::copy(Ain.valuePtr(), Ain.valuePtr() + d_values.size(), d_values.begin()); - SPDLOG_TRACE("[{}] [copy_matrix_to_gpu] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + SPDLOG_INFO("[{}] [copy_matrix_to_gpu] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } { @@ -194,6 +278,14 @@ namespace polysolve::linear bad_indices_arrays.clear(); select_bad_dofs(); + if (contact_patch_schwarz) + { + for (const auto &patch : contact_patches) + { + bad_indices_arrays.emplace_back(patch.begin(), patch.end()); + } + } + if (decompose_subdomains) { filter_subdomains(Ain); @@ -208,26 +300,33 @@ namespace polysolve::linear { decompose_subdomains_to_disjoint_subsets(Ain); } - else + else if (!contact_patch_schwarz) { bad_indices_arrays.emplace_back(h_all_bad_dofs.begin(), h_all_bad_dofs.end()); } - d_all_bad_dofs.clear(); - h_subdomain_sizes.clear(); - d_subdomain_sizes.clear(); - - for (int i = 0; i < bad_indices_arrays.size(); ++i) + if (contact_patch_schwarz) { - d_all_bad_dofs.insert(d_all_bad_dofs.end(), bad_indices_arrays[i].begin(), bad_indices_arrays[i].end()); - h_subdomain_sizes.push_back(bad_indices_arrays[i].size()); + factorize_contact_patches_dense(Ain); } + else + { + d_all_bad_dofs.clear(); + h_subdomain_sizes.clear(); + d_subdomain_sizes.clear(); + + for (int i = 0; i < bad_indices_arrays.size(); ++i) + { + d_all_bad_dofs.insert(d_all_bad_dofs.end(), bad_indices_arrays[i].begin(), bad_indices_arrays[i].end()); + h_subdomain_sizes.push_back(bad_indices_arrays[i].size()); + } - d_subdomain_sizes.insert(d_subdomain_sizes.end(), h_subdomain_sizes.begin(), h_subdomain_sizes.end()); + d_subdomain_sizes.insert(d_subdomain_sizes.end(), h_subdomain_sizes.begin(), h_subdomain_sizes.end()); - factorize_submatrix(); + factorize_submatrix(); + } - SPDLOG_TRACE("[{}] [setup_problematic_dof_precond] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + SPDLOG_INFO("[{}] [setup_problematic_dof_precond] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } if (has_matrix_) @@ -296,8 +395,6 @@ namespace polysolve::linear // Make sure the systems AMG options are set HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim); - // HYPRE_BoomerAMGSetDofFunc(amg_precond, (HYPRE_Int*) dof_to_function.data()); - // More robust options with respect to convergence HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0); HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta); @@ -332,11 +429,29 @@ namespace polysolve::linear precond, dimension_, theta); + + if (block_mapping_.size() > 0) + { + assert(block_mapping_.size() == b.size()); + + // HYPRE_MEMORY_DEVICE is the active memory location here (see the + // constructor), so the array must live on the GPU. Build it on the + // host first, then copy it over with HYPRE's own allocator/copy so + // HYPRE can free it (via hypre_TFree) when the AMG precond is destroyed. + std::vector dof_func_host(b.size()); + for (int i = 0; i < b.size(); ++i) + dof_func_host[i] = block_mapping_[i]; + + HYPRE_Int *dof_func = hypre_CTAlloc(HYPRE_Int, b.size(), HYPRE_MEMORY_DEVICE); + hypre_TMemcpy(dof_func, dof_func_host.data(), HYPRE_Int, b.size(), + HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + HYPRE_BoomerAMGSetDofFunc(precond, dof_func); + } } HYPRE_BoomerAMGSetup(precond, parcsr_A, par_b, par_x); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [amg_setup] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + SPDLOG_INFO("[{}] [amg_setup] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } { @@ -352,7 +467,7 @@ namespace polysolve::linear thrust::copy(d_x.begin(), d_x.end(), x.data()); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [pcg_solve] [{:.6f}] [pcg_iters={}] [residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, final_res_norm); + SPDLOG_INFO("[{}] [pcg_solve] [{:.6f}] [pcg_iters={}] [residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, final_res_norm); } { @@ -405,7 +520,7 @@ namespace polysolve::linear HYPRE_IJMatrixGetObject(A, &temp_A); parcsr_A = static_cast(temp_A); - SPDLOG_TRACE("[{}] [copy_matrix_to_hypre] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + SPDLOG_INFO("[{}] [copy_matrix_to_hypre] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } void GPUHybridSolver::init_hypre_vectors(const int size) @@ -437,7 +552,7 @@ namespace polysolve::linear HYPRE_ParCSRMatrixMatvec(1.0, parcsr_A, par_x, 0.0, par_result); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [matmul] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + GPUHYBRID_LOG_ITER("[{}] [matmul] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } double GPUHybridSolver::dot(const thrust::device_vector &a, const thrust::device_vector &b) @@ -495,7 +610,7 @@ namespace polysolve::linear void GPUHybridSolver::custom_mixed_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &r, thrust::device_vector &z, thrust::device_vector &buffer, thrust::device_vector &z2) { - if (d_all_bad_dofs.size() == 0) + if (d_all_bad_dofs.size() == 0 && num_contact_patches_ == 0) { amg_precond_iter(precond, r, z); } @@ -524,6 +639,12 @@ namespace polysolve::linear void GPUHybridSolver::dss_precond_iter(thrust::device_vector &z, thrust::device_vector &r, thrust::device_vector &next_z) { + if (contact_patch_schwarz) + { + dense_contact_patch_precond_iter(z, r, next_z); + return; + } + auto phase_begin = clock::now(); matmul(z, next_z); @@ -557,7 +678,98 @@ namespace polysolve::linear vector_add(1.0, z, next_z); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [subdomain_solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + GPUHYBRID_LOG_ITER("[{}] [subdomain_solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + void GPUHybridSolver::dense_contact_patch_precond_iter(thrust::device_vector &z, thrust::device_vector &r, thrust::device_vector &next_z) + { + auto phase_begin = clock::now(); + + matmul(z, next_z); + vector_scale(-1.0, next_z); + vector_add(1.0, r, next_z); + + if (num_contact_patches_ > 0) + { + thrust::fill(d_dense_b.begin(), d_dense_b.end(), 0.0); + + const double *raw_next_z = thrust::raw_pointer_cast(next_z.data()); + double *raw_dense_b = thrust::raw_pointer_cast(d_dense_b.data()); + const int *raw_dof_map = thrust::raw_pointer_cast(d_contact_dof_map.data()); + const int *raw_slot = thrust::raw_pointer_cast(d_contact_real_to_dense_slot.data()); + + // Scatter each real dof's current rhs into its patch's padded + // slot; padding slots stay 0 from the fill above, so they solve + // to 0 and never affect the real dofs (see factorize_contact_patches_dense). + thrust::for_each(thrust::device, + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(total_contact_patch_real_dofs_), + [=] __device__(int i) { + raw_dense_b[raw_slot[i]] = raw_next_z[raw_dof_map[i]]; + }); + + // getrsBatched's info is a single HOST scalar (parameter validation + // only), unlike getrfBatched's per-matrix DEVICE info array used in + // factorize_contact_patches_dense -- these are not interchangeable. + int host_info = 0; + CHECK_CUBLAS(cublasDgetrsBatched( + cublas_handle, CUBLAS_OP_N, kContactPatchDenseSize, 1, + thrust::raw_pointer_cast(d_dense_A_ptrs.data()), kContactPatchDenseSize, + thrust::raw_pointer_cast(d_dense_pivot.data()), + thrust::raw_pointer_cast(d_dense_b_ptrs.data()), kContactPatchDenseSize, + &host_info, num_contact_patches_)); + if (host_info != 0) + { + throw std::runtime_error("GPUHybridSolver: cublasDgetrsBatched reported an invalid argument (info=" + std::to_string(host_info) + ")"); + } + } + + thrust::fill(next_z.begin(), next_z.end(), 0.0); + + if (num_contact_patches_ > 0) + { + // Contact patches may overlap, so a dof's solved value can come + // from more than one patch. Gather every real entry's solved + // value (still in the padded per-patch layout) into the order + // sorted by global dof -- precomputed once at factorize time -- + // then sum each dof's duplicate entries with reduce_by_key + // instead of racing to scatter overlapping writes into next_z. + const double *raw_dense_b = thrust::raw_pointer_cast(d_dense_b.data()); + const int *raw_slot = thrust::raw_pointer_cast(d_contact_real_to_dense_slot.data()); + const int *raw_perm = thrust::raw_pointer_cast(d_reduce_perm.data()); + double *raw_reduce_values = thrust::raw_pointer_cast(d_reduce_values.data()); + + thrust::for_each(thrust::device, + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(total_contact_patch_real_dofs_), + [=] __device__(int k) { + int i = raw_perm[k]; + raw_reduce_values[k] = raw_dense_b[raw_slot[i]]; + }); + + thrust::reduce_by_key( + thrust::device, + d_reduce_keys.begin(), d_reduce_keys.end(), + d_reduce_values.begin(), + d_reduce_unique_keys.begin(), d_reduce_unique_values.begin()); + + const double scale = overlap_scale_; + thrust::transform(thrust::device, + d_reduce_unique_values.begin(), d_reduce_unique_values.begin() + num_unique_overlap_dofs_, + d_reduce_unique_values.begin(), + [scale] __device__(double v) { return v * scale; }); + + thrust::scatter( + thrust::device, + d_reduce_unique_values.begin(), d_reduce_unique_values.begin() + num_unique_overlap_dofs_, + d_reduce_unique_keys.begin(), + next_z.begin()); + } + + vector_add(1.0, z, next_z); + + CHECK_CUDA(cudaDeviceSynchronize()); + GPUHYBRID_LOG_ITER("[{}] [contact_patch_dense_subdomain_solve] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } void GPUHybridSolver::amg_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &b, thrust::device_vector &x) @@ -571,7 +783,7 @@ namespace polysolve::linear HYPRE_BoomerAMGSolve(precond, parcsr_A, par_b, par_x); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [amg_v_cycle] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + GPUHYBRID_LOG_ITER("[{}] [amg_v_cycle] [{:.6f}]", name(), elapsed_seconds(phase_begin)); } void GPUHybridSolver::decompose_subdomains_to_disjoint_subsets(const Eigen::SparseMatrix &sparse_A) @@ -616,12 +828,19 @@ namespace polysolve::linear bad_indices_arrays.emplace_back(kv.second.begin(), kv.second.end()); } - SPDLOG_TRACE("[{}] [subdomain_decomposition] [{}] [num_subdomains={}] ", + SPDLOG_INFO("[{}] [subdomain_decomposition] [{}] [num_subdomains={}] ", name(), elapsed_seconds(phase_begin), bad_indices_arrays.size()); } void GPUHybridSolver::select_bad_dofs() { + if (!select_bad_dofs_from_l1_row_norm) + { + // AMGF mode: bad DOFs come only from whatever set_problematic_dofs + // supplied, not from this heuristic. + return; + } + auto phase_begin = clock::now(); const int num_rows = d_outer_indices.size() - 1; @@ -652,92 +871,223 @@ namespace polysolve::linear double var_sum = thrust::transform_reduce(thrust::device, d_row_norms.begin(), d_row_norms.end(), [global_mean] __device__(double x) -> double { return (x - global_mean) * (x - global_mean); }, 0.0, thrust::plus()); double global_var = var_sum / num_rows; - auto minmax = thrust::minmax_element(d_row_norms.begin(), d_row_norms.end()); - double mean_0 = *minmax.first; - double mean_1 = *minmax.second; - double var_0 = global_var; - double var_1 = global_var; - double w0 = 0.5; - double w1 = 0.5; - - double var_reg = 1e-6; + h_all_bad_dofs.clear(); - thrust::device_vector d_gamma0(num_rows); - thrust::device_vector d_gamma1(num_rows); - double *g0 = thrust::raw_pointer_cast(d_gamma0.data()); - double *g1 = thrust::raw_pointer_cast(d_gamma1.data()); + double mean_0 = 0.0, mean_1 = 0.0; + double var_0 = 0.0, var_1 = 0.0; + int gmm_iter = 0; - int gmm_iter; + double max_dist = -1.0; + double max_jump = -1.0; + double min_cost = -1.0; + int split_idx = 0; - for (gmm_iter = 0; gmm_iter < max_gmm_iterations; ++gmm_iter) + switch (subdomain_selection_strategy) + { + case SubdomainSelectionStrategy::KNEE: + case SubdomainSelectionStrategy::FD: + case SubdomainSelectionStrategy::COST: { + // Sort on GPU, preserving original row indices + thrust::device_vector d_sorted_norms = d_row_norms; + thrust::device_vector d_indices(num_rows); + thrust::sequence(d_indices.begin(), d_indices.end()); + + thrust::sort_by_key(d_sorted_norms.begin(), d_sorted_norms.end(), d_indices.begin()); + + // Bring sorted data back to host for O(N) linear scanning + std::vector h_sorted_norms(num_rows); + std::vector h_indices(num_rows); + + thrust::copy(d_sorted_norms.begin(), d_sorted_norms.end(), h_sorted_norms.begin()); + thrust::copy(d_indices.begin(), d_indices.end(), h_indices.begin()); + + int N = num_rows; + if (N >= 2) + { + if (subdomain_selection_strategy == SubdomainSelectionStrategy::KNEE) + { + double x1 = 0.0; + double y1 = std::log(h_sorted_norms[0] + 1e-12); + double x2 = N - 1.0; + double y2 = std::log(h_sorted_norms[N - 1] + 1e-12); + + for (int i = 1; i < N - 1; ++i) + { + double x0 = (double)i; + double y0 = std::log(h_sorted_norms[i] + 1e-12); - double log_likelihood = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { - double x = row_norms[i]; + double dist = std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1); - double log_w0 = log(w0); - double log_w1 = log(w1); - - double log_N0 = -0.5 * log(2.0 * M_PI * var_0) - 0.5 * (x - mean_0) * (x - mean_0) / var_0; - double log_N1 = -0.5 * log(2.0 * M_PI * var_1) - 0.5 * (x - mean_1) * (x - mean_1) / var_1; - - double log_g0 = log_w0 + log_N0; - double log_g1 = log_w1 + log_N1; - - double max_log_g = max(log_g0, log_g1); - double log_total = max_log_g + log(exp(log_g0 - max_log_g) + exp(log_g1 - max_log_g)); - - g0[i] = exp(log_g0 - log_total); - g1[i] = exp(log_g1 - log_total); + if (dist > max_dist) + { + max_dist = dist; + split_idx = i; + } + } + } + else if (subdomain_selection_strategy == SubdomainSelectionStrategy::FD) + { + for (int i = 1; i < N; ++i) + { + double jump = h_sorted_norms[i] - h_sorted_norms[i - 1]; + if (jump > max_jump) + { + max_jump = jump; + split_idx = i; + } + } + } + else if (subdomain_selection_strategy == SubdomainSelectionStrategy::COST) + { + const double c_pcg_conv = 1.0; + const double c_pcg = 1.0; + const double c_fac = 1.0; + const double c_solve = 1.0; - return log_total; }, 0.0, thrust::plus()); + min_cost = std::numeric_limits::max(); - double sum_g0 = thrust::reduce(d_gamma0.begin(), d_gamma0.end(), 0.0); - double sum_g1 = thrust::reduce(d_gamma1.begin(), d_gamma1.end(), 0.0); + double lambda_min = h_sorted_norms[0]; - w0 = sum_g0 / num_rows; - w1 = sum_g1 / num_rows; + for (int i = 0; i < N; ++i) + { + double p = static_cast(N - i); + double lambda_max = h_sorted_norms[i]; - double old_mean_0 = mean_0, old_mean_1 = mean_1; - double old_var_0 = var_0, old_var_1 = var_1; + double m = c_pcg_conv * std::sqrt(lambda_max / lambda_min); + double cost = (m * c_pcg * N) + + (c_fac * p * p) + + (m * c_solve * std::pow(p, 4.0 / 3.0)); - mean_0 = thrust::inner_product(d_gamma0.begin(), d_gamma0.end(), d_row_norms.begin(), 0.0) / sum_g0; - mean_1 = thrust::inner_product(d_gamma1.begin(), d_gamma1.end(), d_row_norms.begin(), 0.0) / sum_g1; + if (cost < min_cost) + { + min_cost = cost; + split_idx = i; + } + } + } - var_0 = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { return g0[i] * (row_norms[i] - mean_0) * (row_norms[i] - mean_0); }, 0.0, thrust::plus()) / sum_g0 + var_reg; + std::vector h_bad_dofs; + h_bad_dofs.reserve(N - split_idx); - var_1 = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { return g1[i] * (row_norms[i] - mean_1) * (row_norms[i] - mean_1); }, 0.0, thrust::plus()) / sum_g1 + var_reg; + for (int i = split_idx; i < N; ++i) + { + h_all_bad_dofs.insert(h_indices[i]); + h_bad_dofs.push_back(h_indices[i]); + } - // Check Convergence - if (abs(mean_0 - old_mean_0) / abs(old_mean_0) < gmm_tol && abs(mean_1 - old_mean_1) / abs(old_mean_1) < gmm_tol && abs(var_0 - old_var_0) / abs(old_var_0) < gmm_tol && abs(var_1 - old_var_1) / abs(old_var_1) < gmm_tol) - { - break; + d_all_bad_dofs = h_bad_dofs; } + break; } - - int num_bad_dofs = 0; - if (abs(mean_1) / abs(mean_0) > gmm_jump_threshold) + case SubdomainSelectionStrategy::GMM: { - d_all_bad_dofs.resize(num_rows); - auto end_it = thrust::copy_if(thrust::device, - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(num_rows), - d_all_bad_dofs.begin(), - [=] __device__(int i) { return g0[i] < g1[i]; }); - - num_bad_dofs = thrust::distance(d_all_bad_dofs.begin(), end_it); - } + auto minmax = thrust::minmax_element(d_row_norms.begin(), d_row_norms.end()); + mean_0 = *minmax.first; + mean_1 = *minmax.second; + var_0 = global_var; + var_1 = global_var; + double w0 = 0.5; + double w1 = 0.5; + + double var_reg = 1e-6; + + thrust::device_vector d_gamma0(num_rows); + thrust::device_vector d_gamma1(num_rows); + double *g0 = thrust::raw_pointer_cast(d_gamma0.data()); + double *g1 = thrust::raw_pointer_cast(d_gamma1.data()); + + for (gmm_iter = 0; gmm_iter < max_gmm_iterations; ++gmm_iter) + { + double log_likelihood = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { + double x = row_norms[i]; - d_all_bad_dofs.resize(num_bad_dofs); + double log_w0 = log(w0); + double log_w1 = log(w1); - std::vector h_bad_dofs(num_bad_dofs); - thrust::copy(d_all_bad_dofs.begin(), d_all_bad_dofs.end(), h_bad_dofs.begin()); + double log_N0 = -0.5 * log(2.0 * M_PI * var_0) - 0.5 * (x - mean_0) * (x - mean_0) / var_0; + double log_N1 = -0.5 * log(2.0 * M_PI * var_1) - 0.5 * (x - mean_1) * (x - mean_1) / var_1; - h_all_bad_dofs.clear(); - h_all_bad_dofs.insert(h_bad_dofs.begin(), h_bad_dofs.end()); + double log_g0 = log_w0 + log_N0; + double log_g1 = log_w1 + log_N1; + + double max_log_g = max(log_g0, log_g1); + double log_total = max_log_g + log(exp(log_g0 - max_log_g) + exp(log_g1 - max_log_g)); + + g0[i] = exp(log_g0 - log_total); + g1[i] = exp(log_g1 - log_total); + + return log_total; }, 0.0, thrust::plus()); + + double sum_g0 = thrust::reduce(d_gamma0.begin(), d_gamma0.end(), 0.0); + double sum_g1 = thrust::reduce(d_gamma1.begin(), d_gamma1.end(), 0.0); + + w0 = sum_g0 / num_rows; + w1 = sum_g1 / num_rows; - SPDLOG_TRACE("[{}] [bad_dof_selection] [{:.6f}] [global_mean={}] [global_var={}] [mean_0={}] [mean_1={}] [var_0={}] [var_1={}] [gmm_iters={}] [num_bad_dofs={}]", - name(), elapsed_seconds(phase_begin), global_mean, global_var, mean_0, mean_1, var_0, var_1, gmm_iter, num_bad_dofs); + double old_mean_0 = mean_0, old_mean_1 = mean_1; + double old_var_0 = var_0, old_var_1 = var_1; + + mean_0 = thrust::inner_product(d_gamma0.begin(), d_gamma0.end(), d_row_norms.begin(), 0.0) / sum_g0; + mean_1 = thrust::inner_product(d_gamma1.begin(), d_gamma1.end(), d_row_norms.begin(), 0.0) / sum_g1; + + var_0 = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { return g0[i] * (row_norms[i] - mean_0) * (row_norms[i] - mean_0); }, 0.0, thrust::plus()) / sum_g0 + var_reg; + + var_1 = thrust::transform_reduce(thrust::device, thrust::make_counting_iterator(0), thrust::make_counting_iterator(num_rows), [=] __device__(int i) -> double { return g1[i] * (row_norms[i] - mean_1) * (row_norms[i] - mean_1); }, 0.0, thrust::plus()) / sum_g1 + var_reg; + + // Check Convergence + if (abs(mean_0 - old_mean_0) / abs(old_mean_0) < gmm_tol && abs(mean_1 - old_mean_1) / abs(old_mean_1) < gmm_tol && abs(var_0 - old_var_0) / abs(old_var_0) < gmm_tol && abs(var_1 - old_var_1) / abs(old_var_1) < gmm_tol) + { + break; + } + } + + int num_bad_dofs = 0; + if (abs(mean_1) / abs(mean_0) > gmm_jump_threshold) + { + d_all_bad_dofs.resize(num_rows); + auto end_it = thrust::copy_if(thrust::device, + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(num_rows), + d_all_bad_dofs.begin(), + [=] __device__(int i) { return g0[i] < g1[i]; }); + + num_bad_dofs = thrust::distance(d_all_bad_dofs.begin(), end_it); + } + + d_all_bad_dofs.resize(num_bad_dofs); + + std::vector h_bad_dofs(num_bad_dofs); + thrust::copy(d_all_bad_dofs.begin(), d_all_bad_dofs.end(), h_bad_dofs.begin()); + + h_all_bad_dofs.insert(h_bad_dofs.begin(), h_bad_dofs.end()); + break; + } + case SubdomainSelectionStrategy::APOSTERIORI: + break; + } + + switch (subdomain_selection_strategy) + { + case SubdomainSelectionStrategy::KNEE: + SPDLOG_INFO("[{}] [bad_dof_selection] [{:.6f}] [strategy=KNEE] [global_mean={}] [global_var={}] [split_idx={}] [max_dist={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, max_dist, h_all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::FD: + SPDLOG_INFO("[{}] [bad_dof_selection] [{:.6f}] [strategy=FD] [global_mean={}] [global_var={}] [split_idx={}] [max_jump={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, max_jump, h_all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::COST: + SPDLOG_INFO("[{}] [bad_dof_selection] [{:.6f}] [strategy=COST] [global_mean={}] [global_var={}] [split_idx={}] [min_cost={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, split_idx, min_cost, h_all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::GMM: + SPDLOG_INFO("[{}] [bad_dof_selection] [{:.6f}] [strategy=GMM] [global_mean={}] [global_var={}] [mean_0={}] [mean_1={}] [var_0={}] [var_1={}] [gmm_iters={}] [num_bad_dofs={}]", + name(), elapsed_seconds(phase_begin), global_mean, global_var, mean_0, mean_1, var_0, var_1, gmm_iter, h_all_bad_dofs.size()); + break; + case SubdomainSelectionStrategy::APOSTERIORI: + break; + } } void GPUHybridSolver::filter_subdomains(const Eigen::SparseMatrix &sparse_A) @@ -823,7 +1173,7 @@ namespace polysolve::linear ++num_not_poorly_conditioned; } - SPDLOG_TRACE("[{}] [subdomain_filtering] [{}] [total_dofs_before={}] [total_dofs_after={}] [num_too_small={}] [num_too_large={}] [num_not_poorly_conditioned={}]", + SPDLOG_INFO("[{}] [subdomain_filtering] [{}] [total_dofs_before={}] [total_dofs_after={}] [num_too_small={}] [num_too_large={}] [num_not_poorly_conditioned={}]", name(), elapsed_seconds(phase_begin), original_num_bad_dofs, h_all_bad_dofs.size(), num_too_small, num_too_large, num_not_poorly_conditioned); } @@ -845,10 +1195,189 @@ namespace polysolve::linear h_all_bad_dofs = std::move(new_bad_dofs); - SPDLOG_TRACE("[{}] [subdomain_expansion] [{}] [num_dofs_before={}] [num_dofs_after={}]", + SPDLOG_INFO("[{}] [subdomain_expansion] [{}] [num_dofs_before={}] [num_dofs_after={}]", name(), elapsed_seconds(phase_begin), num_bad_dofs_before, h_all_bad_dofs.size()); } + // Under contact_patch_schwarz, bad_indices_arrays (the contact patches) + // may overlap, so a covered dof can belong to more than one patch. + // Compute a single global scale factor -- 1 / (average number of patches + // per covered dof) -- used to damp the summed correction in + // dense_contact_patch_precond_iter so overlapping patches don't + // over-correct on average. Also records the number of distinct covered + // dofs, needed to size the reduce_by_key output in that same function. + void GPUHybridSolver::compute_overlap_scale() + { + std::unordered_map dof_patch_count; + for (const auto &subdomain : bad_indices_arrays) + { + for (int dof : subdomain) + { + ++dof_patch_count[dof]; + } + } + + num_unique_overlap_dofs_ = static_cast(dof_patch_count.size()); + + if (dof_patch_count.empty()) + { + overlap_scale_ = 1.0; + return; + } + + long long total_incidences = 0; + for (const auto &kv : dof_patch_count) + { + total_incidences += kv.second; + } + + const double avg_overlap = static_cast(total_incidences) / static_cast(dof_patch_count.size()); + overlap_scale_ = 1.0 / avg_overlap; + + SPDLOG_INFO("[{}] [contact_patch_overlap] [avg_overlap={}] [overlap_scale={}] [num_unique_dofs={}]", + name(), avg_overlap, overlap_scale_, num_unique_overlap_dofs_); + } + + void GPUHybridSolver::factorize_contact_patches_dense(const Eigen::SparseMatrix &sparse_A) + { + auto phase_begin = clock::now(); + + free_contact_patch_dense_memory(); + + constexpr int N = kContactPatchDenseSize; + + num_contact_patches_ = bad_indices_arrays.size(); + + for (const auto &patch : bad_indices_arrays) + { + if ((int)patch.size() > N) + { + throw std::runtime_error( + "GPUHybridSolver: contact patch has " + std::to_string(patch.size()) + + " dofs, which exceeds the fixed dense batched solve size of " + std::to_string(N)); + } + total_contact_patch_real_dofs_ += patch.size(); + } + + if (num_contact_patches_ == 0) + { + SPDLOG_INFO("[{}] [factorize_contact_patches_dense] [{:.6f}] [n_patches=0]", name(), elapsed_seconds(phase_begin)); + return; + } + + // Assemble each patch's dense (padded) block on the host from the + // sparse system matrix, the same way the CPU hybrid solver's + // assemble_D does. This runs once per factorize() call, not on the + // per-iteration hot path, so host-side assembly is fine here even + // though the factorization and every solve happen on the GPU. + std::vector h_dense_A(static_cast(num_contact_patches_) * N * N, 0.0); + std::vector h_dof_map; + std::vector h_real_to_dense_slot; + h_dof_map.reserve(total_contact_patch_real_dofs_); + h_real_to_dense_slot.reserve(total_contact_patch_real_dofs_); + + for (int p = 0; p < num_contact_patches_; ++p) + { + const auto &patch = bad_indices_arrays[p]; + std::unordered_map local_index; + local_index.reserve(patch.size()); + for (int i = 0; i < (int)patch.size(); ++i) + { + local_index[patch[i]] = i; + } + + double *block = h_dense_A.data() + (size_t)p * N * N; + + for (int local_col = 0; local_col < (int)patch.size(); ++local_col) + { + int global_col = patch[local_col]; + for (Eigen::SparseMatrix::InnerIterator it(sparse_A, global_col); it; ++it) + { + auto found = local_index.find(it.row()); + if (found != local_index.end()) + { + // Column-major (lda = N), to match cublasgetrfBatched/getrsBatched. + block[found->second + local_col * N] = it.value(); + } + } + + h_dof_map.push_back(global_col); + h_real_to_dense_slot.push_back(p * N + local_col); + } + + // Pad the unused rows/cols with an identity block, so the padded + // system stays nonsingular and decouples from the real dofs -- + // the padded slots always solve to 0, leaving the real dofs' + // solve exactly as if the patch had been solved unpadded. + for (int i = (int)patch.size(); i < N; ++i) + { + block[i + i * N] = 1.0; + } + } + + d_dense_A = h_dense_A; + d_dense_pivot.resize((size_t)num_contact_patches_ * N); + d_dense_info.resize(num_contact_patches_); + d_dense_b.resize((size_t)num_contact_patches_ * N); + d_contact_dof_map = h_dof_map; + d_contact_real_to_dense_slot = h_real_to_dense_slot; + + std::vector h_dense_A_ptrs(num_contact_patches_); + std::vector h_dense_b_ptrs(num_contact_patches_); + double *raw_dense_A = thrust::raw_pointer_cast(d_dense_A.data()); + double *raw_dense_b = thrust::raw_pointer_cast(d_dense_b.data()); + for (int p = 0; p < num_contact_patches_; ++p) + { + h_dense_A_ptrs[p] = raw_dense_A + (size_t)p * N * N; + h_dense_b_ptrs[p] = raw_dense_b + (size_t)p * N; + } + d_dense_A_ptrs = h_dense_A_ptrs; + d_dense_b_ptrs = h_dense_b_ptrs; + + // Factorize every patch's padded block in one batched call. d_dense_info + // is a device array (one entry per patch), per cublasDgetrfBatched's + // documented convention -- unlike getrsBatched's info below, which is + // a single host scalar. + CHECK_CUBLAS(cublasDgetrfBatched( + cublas_handle, N, + thrust::raw_pointer_cast(d_dense_A_ptrs.data()), N, + thrust::raw_pointer_cast(d_dense_pivot.data()), + thrust::raw_pointer_cast(d_dense_info.data()), + num_contact_patches_)); + CHECK_CUDA(cudaDeviceSynchronize()); + + std::vector h_info(num_contact_patches_); + thrust::copy(d_dense_info.begin(), d_dense_info.end(), h_info.begin()); + for (int p = 0; p < num_contact_patches_; ++p) + { + if (h_info[p] != 0) + { + SPDLOG_WARN("[{}] [contact_patch_singular] [patch={}] [info={}]", name(), p, h_info[p]); + } + } + + compute_overlap_scale(); + + // Precompute the sort-by-global-dof plan used every iteration to sum + // overlapping patches' contributions to a shared dof with + // thrust::reduce_by_key, instead of racing to scatter into it. Since + // bad_indices_arrays (and so d_contact_dof_map) is fixed for the life + // of this factorization, both the permutation and the sorted keys + // only need to be computed once, here. + d_reduce_perm.resize(total_contact_patch_real_dofs_); + thrust::sequence(d_reduce_perm.begin(), d_reduce_perm.end()); + + d_reduce_keys = d_contact_dof_map; + thrust::sort_by_key(d_reduce_keys.begin(), d_reduce_keys.end(), d_reduce_perm.begin()); + + d_reduce_values.resize(total_contact_patch_real_dofs_); + d_reduce_unique_keys.resize(total_contact_patch_real_dofs_); + d_reduce_unique_values.resize(total_contact_patch_real_dofs_); + + SPDLOG_INFO("[{}] [factorize_contact_patches_dense] [{:.6f}] [n_patches={}] [n_real_dofs={}] [n_unique_dofs={}]", + name(), elapsed_seconds(phase_begin), num_contact_patches_, total_contact_patch_real_dofs_, num_unique_overlap_dofs_); + } + void GPUHybridSolver::factorize_submatrix() { auto phase_begin = clock::now(); @@ -1075,7 +1604,7 @@ namespace polysolve::linear } CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [factorize_submatrix] [{}] [n_sparse={}] [n_sparse_dofs={}]", + SPDLOG_INFO("[{}] [factorize_submatrix] [{}] [n_sparse={}] [n_sparse_dofs={}]", name(), elapsed_seconds(phase_begin), sparse_batch_count, total_sparse_dofs); } @@ -1118,7 +1647,7 @@ namespace polysolve::linear gamma = dot(r, z); old_gamma = gamma; - SPDLOG_TRACE("[{}] [pre_loop] [{:.6f}] [rhs_norm={}]", name(), elapsed_seconds(phase_begin), sqrt(bi_prod)); + SPDLOG_INFO("[{}] [pre_loop] [{:.6f}] [rhs_norm={}]", name(), elapsed_seconds(phase_begin), sqrt(bi_prod)); } for (int k = 0; k < max_iter_; ++k) @@ -1131,7 +1660,7 @@ namespace polysolve::linear if (sdotp == 0.0) { - SPDLOG_TRACE("[{}] [err_zero_sdotp] [0.000000]", name()); + SPDLOG_INFO("[{}] [err_zero_sdotp] [0.000000]", name()); break; } @@ -1139,12 +1668,12 @@ namespace polysolve::linear if (alpha <= 0.0) { - SPDLOG_TRACE("[{}] [err_negative_alpha] [0.000000]", name()); + SPDLOG_INFO("[{}] [err_negative_alpha] [0.000000]", name()); break; } else if (alpha < __DBL_MIN__) { - SPDLOG_TRACE("[{}] [err_subnormal_alpha] [0.000000]", name()); + SPDLOG_INFO("[{}] [err_subnormal_alpha] [0.000000]", name()); break; } @@ -1154,13 +1683,13 @@ namespace polysolve::linear if (rel_eps > 0 && (i_prod / bi_prod) < rel_eps) { - SPDLOG_TRACE("[{}] [converged_rel] [0.000000]", name()); + SPDLOG_INFO("[{}] [converged_rel] [0.000000]", name()); break; } if (abs_eps > 0 && i_prod < abs_eps) { - SPDLOG_TRACE("[{}] [converged_abs] [0.000000]", name()); + SPDLOG_INFO("[{}] [converged_abs] [0.000000]", name()); break; } @@ -1176,7 +1705,7 @@ namespace polysolve::linear vector_add(1.0, z, p); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[{}] [pcg_iter] [{:.6f}] [iter={}] [residual={}]", name(), elapsed_seconds(phase_begin), k, sqrt(i_prod)); + GPUHYBRID_LOG_ITER("[{}] [pcg_iter] [{:.6f}] [iter={}] [residual={}]", name(), elapsed_seconds(phase_begin), k, sqrt(i_prod)); } } @@ -1190,12 +1719,39 @@ namespace polysolve::linear } free_device_memory(); + free_contact_patch_dense_memory(); if (cudss_handle) { cudssDestroy(cudss_handle); cudss_handle = nullptr; } + if (cublas_handle) + { + cublasDestroy(cublas_handle); + cublas_handle = nullptr; + } + } + + void GPUHybridSolver::free_contact_patch_dense_memory() + { + num_contact_patches_ = 0; + total_contact_patch_real_dofs_ = 0; + num_unique_overlap_dofs_ = 0; + + d_dense_A.clear(); + d_dense_A_ptrs.clear(); + d_dense_pivot.clear(); + d_dense_info.clear(); + d_dense_b.clear(); + d_dense_b_ptrs.clear(); + d_contact_dof_map.clear(); + d_contact_real_to_dense_slot.clear(); + d_reduce_perm.clear(); + d_reduce_keys.clear(); + d_reduce_values.clear(); + d_reduce_unique_keys.clear(); + d_reduce_unique_values.clear(); } void GPUHybridSolver::free_device_memory() diff --git a/src/polysolve/linear/GPUHybridSolver.hpp b/src/polysolve/linear/GPUHybridSolver.hpp index b71d656a..a44b5681 100644 --- a/src/polysolve/linear/GPUHybridSolver.hpp +++ b/src/polysolve/linear/GPUHybridSolver.hpp @@ -13,6 +13,7 @@ #include #include +#include #include #include @@ -43,6 +44,12 @@ namespace polysolve::linear // Set solver parameters virtual void set_parameters(const json ¶ms) override; + // Set block size for multigrid solvers + virtual void set_block_size(int block_size) override; + + // Set the function (block) assigned to each row for multigrid solvers + virtual void set_block_mapping(const Eigen::VectorXi &block_mapping) override; + // Retrieve solve information virtual void get_info(json ¶ms) const override; @@ -57,14 +64,32 @@ namespace polysolve::linear // Name of the solver type (for debugging purposes) virtual std::string name() const override { - return "GPUHybrid"; + return select_bad_dofs_from_l1_row_norm ? "GPUHybrid" : "GPUAMGF"; } + virtual void set_problematic_dofs(const std::set &bad_dofs) override { h_all_bad_dofs = bad_dofs; } + + // Alternative strategies for picking which DOFs get the dense/direct + // subdomain treatment, in place of the default GMM row-norm clustering. + enum class SubdomainSelectionStrategy + { + KNEE, + GMM, + FD, + COST, + APOSTERIORI + }; + protected: // AMG settings double theta = 0.5; // Hybrid preconditioner settings + // When false, bad DOFs are taken only from whatever set_problematic_dofs + // supplies instead of being inferred from the row-norm heuristic below + // (AMGF mode). + bool select_bad_dofs_from_l1_row_norm = true; + SubdomainSelectionStrategy subdomain_selection_strategy = SubdomainSelectionStrategy::GMM; bool decompose_subdomains = true; int min_subdomain_size = 1; int max_subdomain_size = 1e9; @@ -73,14 +98,31 @@ namespace polysolve::linear int max_gmm_iterations = 20; bool expand_subdomains = true; bool additive_mode = false; + // Give each contact patch its own (possibly overlapping) subdomain, + // solved additively -- an additive Schwarz treatment of contact DOFs. + // Unlike the generic (row-norm/GMM-selected) subdomains, contact + // patches are factorized and solved with a fixed-size dense batched + // GPU solve rather than the sparse batched cuDSS path -- see + // factorize_contact_patches_dense(). + bool contact_patch_schwarz = false; // General solver settings int dimension_ = 1; // 1 = scalar (Laplace), 2 or 3 = vector (Elasticity) + + // Optional per-row function (block) assignment for HYPRE_BoomerAMGSetDofFunc. + // Empty means the default interleaved mapping (row i belongs to function i % dimension_). + Eigen::VectorXi block_mapping_; int max_iter_ = 10000; double rel_conv_tol_ = 1e-10; double abs_conv_tol_ = 0.0; double conditioning_threshold = 100.0; + // When false, skips the per-PCG-iteration timing logs (matmul, + // amg_v_cycle, subdomain_solve, contact_patch_dense_subdomain_solve, + // pcg_iter) so a long solve doesn't flood the log; every other log + // line is unaffected. + bool detailed_log = true; + // solve information HYPRE_Int num_iterations; HYPRE_Complex final_res_norm; @@ -108,11 +150,12 @@ namespace polysolve::linear cudssHandle_t cudss_handle = nullptr; cudssConfig_t cudss_config = nullptr; cudssData_t cudss_solver_data = nullptr; - cudssMatrix_t batch_A = nullptr; cudssMatrix_t batch_x = nullptr; cudssMatrix_t batch_b = nullptr; + cublasHandle_t cublas_handle = nullptr; + int sparse_batch_count = 0; thrust::device_vector d_sparse_dof_map; @@ -126,6 +169,45 @@ namespace polysolve::linear thrust::device_vector d_sparse_inner_void, d_sparse_outer_void, d_sparse_values_void; thrust::device_vector d_sparse_x_void, d_sparse_b_void; + // contact_patch_schwarz: fixed-size dense batched solve data. Every + // contact patch is padded (with an identity block on the unused + // rows/cols) up to this many dofs and factorized/solved as one dense + // block per patch via cublasgetrfBatched/getrsBatched -- see + // factorize_contact_patches_dense(). + static constexpr int kContactPatchDenseSize = 64; + int num_contact_patches_ = 0; + int total_contact_patch_real_dofs_ = 0; + + thrust::device_vector d_dense_A; // num_contact_patches_ * 64 * 64, column-major, LU factors in place + thrust::device_vector d_dense_A_ptrs; // num_contact_patches_ device pointers into d_dense_A + thrust::device_vector d_dense_pivot; // num_contact_patches_ * 64 + thrust::device_vector d_dense_info; // num_contact_patches_ (getrfBatched status, device) + + thrust::device_vector d_dense_b; // num_contact_patches_ * 64, padded rhs/solution + thrust::device_vector d_dense_b_ptrs; // num_contact_patches_ device pointers into d_dense_b + + // Per real (non-padding) dof, patch-major: its global dof index, and + // its flat slot (patch * 64 + local_index) in the padded layout above. + thrust::device_vector d_contact_dof_map; + thrust::device_vector d_contact_real_to_dense_slot; + + // Contact patches may overlap, so a dof's correction is the sum of + // every patch that covers it. This global (not per-dof) scale factor + // -- 1 / (average number of patches covering a covered dof) -- damps + // that sum so overlapping patches don't over-correct on average. + double overlap_scale_ = 1.0; + int num_unique_overlap_dofs_ = 0; + + // Scratch for summing overlapping patches' contributions to a shared + // dof without atomics: sort real entries by global dof once at + // factorize time (d_reduce_perm/d_reduce_keys), then each iteration + // gather solved values into that order and thrust::reduce_by_key them. + thrust::device_vector d_reduce_perm; + thrust::device_vector d_reduce_keys; + thrust::device_vector d_reduce_values; + thrust::device_vector d_reduce_unique_keys; + thrust::device_vector d_reduce_unique_values; + public: void free_device_memory(); @@ -141,6 +223,9 @@ namespace polysolve::linear void expand_subdomains_to_strongly_connected(const Eigen::SparseMatrix &sparse_A); void select_bad_dofs(); void factorize_submatrix(); + void factorize_contact_patches_dense(const Eigen::SparseMatrix &sparse_A); + void compute_overlap_scale(); + void free_contact_patch_dense_memory(); // linear algebra void set_hypre_vec(HYPRE_IJVector &ij_x, HYPRE_ParVector &par_x, const thrust::device_vector &x); @@ -154,6 +239,7 @@ namespace polysolve::linear void custom_mixed_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &r, thrust::device_vector &z, thrust::device_vector &buffer, thrust::device_vector &z2); void amg_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &b, thrust::device_vector &x); void dss_precond_iter(thrust::device_vector &z, thrust::device_vector &r, thrust::device_vector &next_z); + void dense_contact_patch_precond_iter(thrust::device_vector &z, thrust::device_vector &r, thrust::device_vector &next_z); // Krylov solve methods void pcg_solve(thrust::device_vector &rhs, thrust::device_vector &result, HYPRE_ParVector &par_b, HYPRE_ParVector &par_x, HYPRE_Solver &precond); diff --git a/src/polysolve/linear/HypreSolver.cpp b/src/polysolve/linear/HypreSolver.cpp index 974125d1..2645e849 100644 --- a/src/polysolve/linear/HypreSolver.cpp +++ b/src/polysolve/linear/HypreSolver.cpp @@ -5,6 +5,9 @@ #include #include +// Internal header, needed for hypre_CTAlloc/hypre_TFree so the dof_func array +// handed to HYPRE_BoomerAMGSetDofFunc uses the allocator HYPRE itself frees it with. +#include "_hypre_utilities.h" //////////////////////////////////////////////////////////////////////////////// namespace polysolve::linear @@ -61,13 +64,19 @@ namespace polysolve::linear { interp_rbms = params["Hypre"]["interp_rbms"]; } - if (params["Hypre"].contains("dimension")) - { - dimension_ = params["Hypre"]["dimension"]; - } } } + void HypreSolver::set_block_size(int block_size) + { + dimension_ = block_size; + } + + void HypreSolver::set_block_mapping(const Eigen::VectorXi &block_mapping) + { + block_mapping_ = block_mapping; + } + void HypreSolver::get_info(json ¶ms) const { params["num_iterations"] = num_iterations; @@ -308,6 +317,19 @@ namespace polysolve::linear Eigen::MatrixXd positions_; assert(!interp_rbms); HypreBoomerAMG_SetElasticityOptions(precond, dimension_, theta, nodal_coarsening, interp_rbms, positions_, rbms, par_rbms); + + if (block_mapping_.size() > 0) + { + assert(block_mapping_.size() == rhs.size()); + + // HYPRE takes ownership of this array and frees it when the AMG + // preconditioner is destroyed, so it must be allocated with + // HYPRE's own allocator rather than new[]/malloc. + HYPRE_Int *dof_func = hypre_CTAlloc(HYPRE_Int, rhs.size(), HYPRE_MEMORY_HOST); + for (HYPRE_Int i = 0; i < rhs.size(); ++i) + dof_func[i] = block_mapping_[i]; + HYPRE_BoomerAMGSetDofFunc(precond, dof_func); + } } /* Set the PCG preconditioner */ diff --git a/src/polysolve/linear/HypreSolver.hpp b/src/polysolve/linear/HypreSolver.hpp index e8b0c763..c7ab0b8c 100644 --- a/src/polysolve/linear/HypreSolver.hpp +++ b/src/polysolve/linear/HypreSolver.hpp @@ -38,6 +38,12 @@ namespace polysolve::linear // Set solver parameters virtual void set_parameters(const json ¶ms) override; + // Set block size for multigrid solvers + virtual void set_block_size(int block_size) override; + + // Set the function (block) assigned to each row for multigrid solvers + virtual void set_block_mapping(const Eigen::VectorXi &block_mapping) override; + // Retrieve memory information from Pardiso virtual void get_info(json ¶ms) const override; @@ -73,6 +79,10 @@ namespace polysolve::linear bool has_matrix_ = false; int precond_num_; + // Optional per-row function (block) assignment for HYPRE_BoomerAMGSetDofFunc. + // Empty means the default interleaved mapping (row i belongs to function i % dimension_). + Eigen::VectorXi block_mapping_; + HYPRE_IJMatrix A; HYPRE_ParCSRMatrix parcsr_A; }; diff --git a/src/polysolve/linear/MASSolver.cu b/src/polysolve/linear/MASSolver.cu index e73155f5..c738c661 100644 --- a/src/polysolve/linear/MASSolver.cu +++ b/src/polysolve/linear/MASSolver.cu @@ -197,8 +197,6 @@ namespace polysolve::linear void set_parameters(const json ¶ms) { - if (params.contains("block_dim")) - block_dim_ = params["block_dim"]; if (params.contains("max_iter")) max_iter_ = params["max_iter"]; if (params.contains("relative_tolerance")) @@ -211,6 +209,11 @@ namespace polysolve::linear use_preconditioned_residual_norm_ = params["use_preconditioned_residual_norm"]; } + void set_block_size(int block_size) + { + block_dim_ = block_size; + } + void get_info(json ¶ms) const { params["solver_iter"] = iterations_; @@ -308,16 +311,16 @@ namespace polysolve::linear // Build adjacency for graph partition. auto phase_begin = clock::now(); BSRAdjacency adj{A, block_dim_, rt}; - SPDLOG_TRACE("[MAS] [factorize_topology_adj] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_topology_adj] [{:.6f}]", elapsed_seconds(phase_begin)); // Sort nodes based on parition. phase_begin = clock::now(); build_partition_and_perm(adj); - SPDLOG_TRACE("[MAS] [factorize_graph_partition] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_graph_partition] [{:.6f}]", elapsed_seconds(phase_begin)); } else { - SPDLOG_TRACE("[MAS] [factorize_reuse_partition] [0.000000]"); + SPDLOG_INFO("[MAS] [factorize_reuse_partition] [0.000000]"); } // Build new sorted BSR matrix for MAS initialization. @@ -327,7 +330,7 @@ namespace polysolve::linear block_dim_, ctd::span(permutation_.data(), permutation_.size()), rt}; - SPDLOG_TRACE("[MAS] [factorize_permuted_bsr] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_permuted_bsr] [{:.6f}]", elapsed_seconds(phase_begin)); BSRView view = A_.view(); dim_ = A.rows(); @@ -340,7 +343,7 @@ namespace polysolve::linear ctd::span(part_offsets_.data(), part_offsets_.size()), rt); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_mas] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_mas] [{:.6f}]", elapsed_seconds(phase_begin)); // Copy permutation to device. phase_begin = clock::now(); @@ -363,14 +366,14 @@ namespace polysolve::linear scalar_rz_old_ = safe_alloc(1, rt, "device_buffers scalar_rz_old"); scalar_rr_ = safe_alloc(1, rt, "device_buffers scalar_rr"); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_device_buffers] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_device_buffers] [{:.6f}]", elapsed_seconds(phase_begin)); // Allocates buffers for CuSparse. phase_begin = clock::now(); setup_cusparse(rt); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_cusparse] [{:.6f}]", elapsed_seconds(phase_begin)); - SPDLOG_TRACE("[MAS] [factorize_total] [{:.6f}]", elapsed_seconds(total_begin)); + SPDLOG_INFO("[MAS] [factorize_cusparse] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_total] [{:.6f}]", elapsed_seconds(total_begin)); } void solve(const Eigen::Ref b, Eigen::Ref x) @@ -563,7 +566,7 @@ namespace polysolve::linear if (k % 100 == 0) { rt.stream.sync(); - SPDLOG_TRACE( + SPDLOG_INFO( "[MAS] [pcg_loop] [{:.6f}] [iter={}] [residual={:.6e}]", elapsed_seconds(iter_window_begin), k, @@ -587,7 +590,7 @@ namespace polysolve::linear status_ = MASSolverStatus::ReachMaxIterations; } - SPDLOG_TRACE( + SPDLOG_INFO( "[MAS] [pcg_loop] [{:.6f}] [iter={}] [residual={:.6e}]", elapsed_seconds(iter_window_begin), iterations_, @@ -613,6 +616,11 @@ namespace polysolve::linear impl_->set_parameters(params[solver_name]); } + void MASSolver::set_block_size(int block_size) + { + impl_->set_block_size(block_size); + } + void MASSolver::get_info(json ¶ms) const { impl_->get_info(params); diff --git a/src/polysolve/linear/MASSolver.hpp b/src/polysolve/linear/MASSolver.hpp index 8124f54b..b0ca7771 100644 --- a/src/polysolve/linear/MASSolver.hpp +++ b/src/polysolve/linear/MASSolver.hpp @@ -48,6 +48,9 @@ namespace polysolve::linear // Set solver parameters void set_parameters(const json ¶ms) override; + // Set block size for multigrid solvers + void set_block_size(int block_size) override; + // Retrieve information void get_info(json ¶ms) const override; diff --git a/src/polysolve/linear/RestartedAMGPCG.cu b/src/polysolve/linear/RestartedAMGPCG.cu new file mode 100644 index 00000000..b34b5ca3 --- /dev/null +++ b/src/polysolve/linear/RestartedAMGPCG.cu @@ -0,0 +1,637 @@ +#include "RestartedAMGPCG.hpp" + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +// Internal header, needed for hypre_CTAlloc/hypre_TMemcpy so the dof_func +// array handed to HYPRE_BoomerAMGSetDofFunc uses the allocator HYPRE itself +// frees it with (and lives in the device memory HYPRE expects here). +#include "_hypre_utilities.h" + +#ifdef HYPRE_ENABLE_MPI +#include +#endif + +#define CHECK_CUDA(call) \ + do \ + { \ + cudaError_t status = call; \ + if (status != cudaSuccess) \ + { \ + std::cerr << "CUDA Error at " << __FILE__ << ":" << __LINE__ \ + << " - " << cudaGetErrorName(status) \ + << " (" << cudaGetErrorString(status) << ")" << std::endl; \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + +// Per-PCG-iteration logs (matmul, amg_v_cycle, pcg_iter) go through this +// instead of SPDLOG_INFO directly, so they can be silenced with +// detailed_log=false without affecting any other log line. +#define RESTARTED_AMG_PCG_LOG_ITER(...) \ + do \ + { \ + if (detailed_log) \ + SPDLOG_INFO(__VA_ARGS__); \ + } while (0) + +namespace polysolve::linear +{ + + namespace + { + using clock = std::chrono::steady_clock; + + double elapsed_seconds(const std::chrono::time_point &begin) + { + return std::chrono::duration(clock::now() - begin).count(); + } + + // Error-free transformations: s + e == a + b and p + e == a * b + // exactly (Knuth's TwoSum; TwoProduct via FMA). The explicit + // __dadd_rn/__dsub_rn/__dmul_rn intrinsics are never contracted into + // an FMA by nvcc, unlike plain +, -, *; such a contraction (e.g. of + // a * b + c) would silently break the exactness these rely on. + __device__ inline void two_sum(const double a, const double b, double &s, double &e) + { + s = __dadd_rn(a, b); + const double bb = __dsub_rn(s, a); + e = __dadd_rn(__dsub_rn(a, __dsub_rn(s, bb)), __dsub_rn(b, bb)); + } + + __device__ inline void two_prod(const double a, const double b, double &p, double &e) + { + p = __dmul_rn(a, b); + e = __fma_rn(a, b, -p); + } + } // namespace + + RestartedAMGPCG::RestartedAMGPCG() + { +#ifdef HYPRE_ENABLE_MPI + int done_already; + + MPI_Initialized(&done_already); + if (!done_already) + { + MPI_Init(nullptr, nullptr); + } +#endif + if (!HYPRE_Initialized()) + { + HYPRE_Initialize(); + } + + HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE); + HYPRE_SetExecutionPolicy(HYPRE_EXEC_DEVICE); + HYPRE_SetSpGemmUseCusparse(false); + HYPRE_SetUseGpuRand(true); + } + + void RestartedAMGPCG::set_parameters(const json ¶ms) + { + if (params.contains("RestartedAMGPCG")) + { + if (params["RestartedAMGPCG"].contains("max_iter")) + { + max_iter_ = params["RestartedAMGPCG"]["max_iter"]; + } + if (params["RestartedAMGPCG"].contains("relative_tolerance")) + { + rel_conv_tol_ = params["RestartedAMGPCG"]["relative_tolerance"]; + } + if (params["RestartedAMGPCG"].contains("absolute_tolerance")) + { + abs_conv_tol_ = params["RestartedAMGPCG"]["absolute_tolerance"]; + } + if (params["RestartedAMGPCG"].contains("theta")) + { + theta = params["RestartedAMGPCG"]["theta"]; + } + if (params["RestartedAMGPCG"].contains("restart_interval")) + { + const int interval = params["RestartedAMGPCG"]["restart_interval"]; + if (interval < 1) + { + throw std::runtime_error("RestartedAMGPCG: restart_interval must be at least 1, got " + std::to_string(interval)); + } + restart_interval = interval; + } + if (params["RestartedAMGPCG"].contains("detailed_log")) + { + detailed_log = params["RestartedAMGPCG"]["detailed_log"]; + } + } + } + + void RestartedAMGPCG::set_block_size(int block_size) + { + dimension_ = block_size; + } + + void RestartedAMGPCG::set_block_mapping(const Eigen::VectorXi &block_mapping) + { + block_mapping_ = block_mapping; + } + + void RestartedAMGPCG::get_info(json ¶ms) const + { + params["num_iterations"] = num_iterations; + params["num_restarts"] = num_restarts; + params["final_res_norm"] = final_res_norm; + } + + void RestartedAMGPCG::factorize(const StiffnessMatrix &Ain) + { + SPDLOG_INFO("[{}] [start_solve] [0.000000] [problem_size={}]", name(), Ain.rows()); + + { + auto phase_begin = clock::now(); + + d_outer_indices.resize(Ain.rows() + 1); + d_inner_indices.resize(Ain.nonZeros()); + d_values.resize(Ain.nonZeros()); + + thrust::copy(Ain.outerIndexPtr(), Ain.outerIndexPtr() + d_outer_indices.size(), d_outer_indices.begin()); + thrust::copy(Ain.innerIndexPtr(), Ain.innerIndexPtr() + d_inner_indices.size(), d_inner_indices.begin()); + thrust::copy(Ain.valuePtr(), Ain.valuePtr() + d_values.size(), d_values.begin()); + + SPDLOG_INFO("[{}] [copy_matrix_to_gpu] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + if (has_matrix_) + { + HYPRE_IJMatrixDestroy(A); + has_matrix_ = false; + A = nullptr; + } + + copy_matrix_to_hypre(); + has_matrix_ = true; + } + + namespace + { + // Same BoomerAMG configuration as GPUHybridSolver. + void HypreBoomerAMG_SetDefaultOptions(HYPRE_Solver &amg_precond) + { + // AMG coarsening options: + int coarsen_type = 8; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP + int agg_levels = 1; // number of aggressive coarsening levels + double theta = 0.25; // strength threshold: 0.25, 0.5, 0.8 + + // AMG interpolation options: + int interp_type = 6; // 6 = extended+i, 0 = classical + int Pmax = 4; // max number of elements per row in P + + // AMG relaxation options: + int relax_type = 18; // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi + int relax_sweeps = 1; // relaxation sweeps on each level + + // Additional options: + int print_level = 0; // print AMG iterations? 1 = no, 2 = yes + int max_levels = 25; // max number of levels in AMG hierarchy + + int min_coarse_size = 5; + + HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); + HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels); + HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); + + HYPRE_BoomerAMGSetRelaxOrder(amg_precond, false); + HYPRE_BoomerAMGSetRAP2(amg_precond, true); + HYPRE_BoomerAMGSetKeepTranspose(amg_precond, true); + + HYPRE_BoomerAMGSetMinCoarseSize(amg_precond, min_coarse_size); + HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_type, 3); + HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps); + HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta); + HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); + HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax); + HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); + HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); + + // Use as a preconditioner (one V-cycle, zero tolerance) + HYPRE_BoomerAMGSetMaxIter(amg_precond, 1); + HYPRE_BoomerAMGSetTol(amg_precond, 0.0); + } + + void HypreBoomerAMG_SetElasticityOptions(HYPRE_Solver &amg_precond, int dim, double theta) + { + // Make sure the systems AMG options are set + HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim); + + // More robust options with respect to convergence + HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0); + HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta); + } + } // namespace + + void RestartedAMGPCG::solve(const Ref b, Ref x) + { + thrust::device_vector d_x(x.size()); + thrust::device_vector d_b(b.size()); + + thrust::copy(x.data(), x.data() + x.size(), d_x.begin()); + thrust::copy(b.data(), b.data() + b.size(), d_b.begin()); + + HYPRE_ParVector par_b; + HYPRE_ParVector par_x; + init_hypre_vectors(b.size()); + + set_hypre_vec(ij_b, par_b, d_b); + set_hypre_vec(ij_x, par_x, d_x); + + HYPRE_Solver precond; + + { + auto phase_begin = clock::now(); + + HYPRE_BoomerAMGCreate(&precond); + HypreBoomerAMG_SetDefaultOptions(precond); + if (dimension_ > 1) + { + HypreBoomerAMG_SetElasticityOptions( + precond, + dimension_, + theta); + + if (block_mapping_.size() > 0) + { + assert(block_mapping_.size() == b.size()); + + // HYPRE_MEMORY_DEVICE is the active memory location here (see the + // constructor), so the array must live on the GPU. Build it on the + // host first, then copy it over with HYPRE's own allocator/copy so + // HYPRE can free it (via hypre_TFree) when the AMG precond is destroyed. + std::vector dof_func_host(b.size()); + for (int i = 0; i < b.size(); ++i) + dof_func_host[i] = block_mapping_[i]; + + HYPRE_Int *dof_func = hypre_CTAlloc(HYPRE_Int, b.size(), HYPRE_MEMORY_DEVICE); + hypre_TMemcpy(dof_func, dof_func_host.data(), HYPRE_Int, b.size(), + HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + HYPRE_BoomerAMGSetDofFunc(precond, dof_func); + } + } + + HYPRE_BoomerAMGSetup(precond, parcsr_A, par_b, par_x); + CHECK_CUDA(cudaDeviceSynchronize()); + SPDLOG_INFO("[{}] [amg_setup] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + { + auto phase_begin = clock::now(); + + // Also sets final_res_norm, from the double-double residual it + // last checked for convergence. + pcg_solve(d_b, d_x, precond); + + thrust::copy(d_x.begin(), d_x.end(), x.data()); + + CHECK_CUDA(cudaDeviceSynchronize()); + SPDLOG_INFO("[{}] [pcg_solve] [{:.6f}] [pcg_iters={}] [restarts={}] [residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, num_restarts, final_res_norm); + } + + { + HYPRE_BoomerAMGDestroy(precond); + HYPRE_IJVectorDestroy(ij_x); + HYPRE_IJVectorDestroy(ij_b); + } + } + + void RestartedAMGPCG::copy_matrix_to_hypre() + { + auto phase_begin = clock::now(); + + const HYPRE_Int num_rows = d_outer_indices.size() - 1; + +#ifdef HYPRE_ENABLE_MPI + HYPRE_IJMatrixCreate(MPI_COMM_WORLD, 0, num_rows - 1, 0, num_rows - 1, &A); +#else + HYPRE_IJMatrixCreate(0, 0, num_rows - 1, 0, num_rows - 1, &A); +#endif + HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR); + HYPRE_IJMatrixInitialize(A); + + thrust::device_vector d_rows(num_rows); + thrust::sequence(d_rows.begin(), d_rows.end()); + + thrust::device_vector d_n_cols(num_rows); + const HYPRE_Int *raw_outer = thrust::raw_pointer_cast(d_outer_indices.data()); + HYPRE_Int *raw_n_cols = thrust::raw_pointer_cast(d_n_cols.data()); + + thrust::for_each(thrust::device, + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(num_rows), + [=] __device__(int i) { + raw_n_cols[i] = raw_outer[i + 1] - raw_outer[i]; + }); + + HYPRE_Int *gpu_n_cols = thrust::raw_pointer_cast(d_n_cols.data()); + HYPRE_Int *gpu_rows = thrust::raw_pointer_cast(d_rows.data()); + + HYPRE_Int *gpu_cols = thrust::raw_pointer_cast(d_inner_indices.data()); + double *gpu_vals = thrust::raw_pointer_cast(d_values.data()); + + HYPRE_IJMatrixSetValues(A, num_rows, gpu_n_cols, gpu_rows, gpu_cols, gpu_vals); + + HYPRE_IJMatrixAssemble(A); + + void *temp_A = nullptr; + HYPRE_IJMatrixGetObject(A, &temp_A); + parcsr_A = static_cast(temp_A); + + SPDLOG_INFO("[{}] [copy_matrix_to_hypre] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + void RestartedAMGPCG::init_hypre_vectors(const int size) + { +#ifdef HYPRE_ENABLE_MPI + HYPRE_IJVectorCreate(MPI_COMM_WORLD, 0, size - 1, &ij_x); +#else + HYPRE_IJVectorCreate(0, 0, size - 1, &ij_x); +#endif + HYPRE_IJVectorSetObjectType(ij_x, HYPRE_PARCSR); + HYPRE_IJVectorInitializeShell(ij_x); +#ifdef HYPRE_ENABLE_MPI + HYPRE_IJVectorCreate(MPI_COMM_WORLD, 0, size - 1, &ij_b); +#else + HYPRE_IJVectorCreate(0, 0, size - 1, &ij_b); +#endif + HYPRE_IJVectorSetObjectType(ij_b, HYPRE_PARCSR); + HYPRE_IJVectorInitializeShell(ij_b); + } + + void RestartedAMGPCG::matmul(const thrust::device_vector &x, thrust::device_vector &result) + { + auto phase_begin = clock::now(); + HYPRE_ParVector par_x; + HYPRE_ParVector par_result; + + set_hypre_vec(ij_x, par_x, x); + set_hypre_vec(ij_b, par_result, result); + + HYPRE_ParCSRMatrixMatvec(1.0, parcsr_A, par_x, 0.0, par_result); + CHECK_CUDA(cudaDeviceSynchronize()); + RESTARTED_AMG_PCG_LOG_ITER("[{}] [matmul] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + double RestartedAMGPCG::dot(const thrust::device_vector &a, const thrust::device_vector &b) + { + HYPRE_ParVector par_a; + HYPRE_ParVector par_b; + + set_hypre_vec(ij_x, par_a, a); + set_hypre_vec(ij_b, par_b, b); + + double result; + HYPRE_ParVectorInnerProd(par_a, par_b, &result); + return result; + } + + void RestartedAMGPCG::vector_copy(const thrust::device_vector &x, thrust::device_vector &y) + { + HYPRE_ParVector par_x; + HYPRE_ParVector par_y; + + set_hypre_vec(ij_x, par_x, x); + set_hypre_vec(ij_b, par_y, y); + + HYPRE_ParVectorCopy(par_x, par_y); + } + + void RestartedAMGPCG::vector_add(double alpha, const thrust::device_vector &x, thrust::device_vector &y) + { + HYPRE_ParVector par_x; + HYPRE_ParVector par_y; + + set_hypre_vec(ij_x, par_x, x); + set_hypre_vec(ij_b, par_y, y); + + hypre_ParVectorAxpy(alpha, par_x, par_y); + } + + void RestartedAMGPCG::vector_scale(double alpha, thrust::device_vector &x) + { + HYPRE_ParVector par_x; + + set_hypre_vec(ij_x, par_x, x); + + HYPRE_ParVectorScale(alpha, par_x); + } + + void RestartedAMGPCG::set_hypre_vec(HYPRE_IJVector &my_ij_x, HYPRE_ParVector &par_x, const thrust::device_vector &x) + { + double *raw_ptr = const_cast(thrust::raw_pointer_cast(x.data())); + + HYPRE_IJVectorSetData(my_ij_x, raw_ptr); + HYPRE_IJVectorAssemble(my_ij_x); + HYPRE_IJVectorGetObject(my_ij_x, (void **)&par_x); + } + + void RestartedAMGPCG::dd_residual(const thrust::device_vector &b, const thrust::device_vector &x, thrust::device_vector &r) + { + const int num_rows = d_outer_indices.size() - 1; + + const int *raw_outer = thrust::raw_pointer_cast(d_outer_indices.data()); + const int *raw_inner = thrust::raw_pointer_cast(d_inner_indices.data()); + const double *raw_values = thrust::raw_pointer_cast(d_values.data()); + const double *raw_b = thrust::raw_pointer_cast(b.data()); + const double *raw_x = thrust::raw_pointer_cast(x.data()); + double *raw_r = thrust::raw_pointer_cast(r.data()); + + // Ogita, Rump & Oishi's Dot2 ("Accurate sum and dot product", 2005), + // seeded with b_i: the rounding error of every product and partial + // sum is recovered exactly and accumulated into lo, so hi + lo is as + // accurate as if b_i - (Ax)_i were computed in twice the working + // precision and rounded once -- even when b_i and (Ax)_i nearly + // cancel, which is exactly when a double-precision residual fails. + thrust::for_each(thrust::device, + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(num_rows), + [=] __device__(int i) { + double hi = raw_b[i]; + double lo = 0.0; + for (int k = raw_outer[i]; k < raw_outer[i + 1]; ++k) + { + double prod, prod_err, sum_err; + two_prod(-raw_values[k], raw_x[raw_inner[k]], prod, prod_err); + two_sum(hi, prod, hi, sum_err); + lo += sum_err + prod_err; + } + raw_r[i] = hi + lo; + }); + } + + void RestartedAMGPCG::amg_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &b, thrust::device_vector &x) + { + auto phase_begin = clock::now(); + HYPRE_ParVector par_x; + HYPRE_ParVector par_b; + + set_hypre_vec(ij_x, par_x, x); + set_hypre_vec(ij_b, par_b, b); + + HYPRE_BoomerAMGSolve(precond, parcsr_A, par_b, par_x); + CHECK_CUDA(cudaDeviceSynchronize()); + RESTARTED_AMG_PCG_LOG_ITER("[{}] [amg_v_cycle] [{:.6f}]", name(), elapsed_seconds(phase_begin)); + } + + void RestartedAMGPCG::pcg_solve(thrust::device_vector &rhs, thrust::device_vector &result, HYPRE_Solver &precond) + { + thrust::device_vector r(rhs.size()); + thrust::device_vector p(rhs.size()); + thrust::device_vector z(rhs.size()); + thrust::device_vector buffer(rhs.size()); + + num_iterations = 0; + num_restarts = 0; + + const double bi_prod = dot(rhs, rhs); + if (bi_prod <= 0.0) + { + thrust::fill(result.begin(), result.end(), 0.0); + final_res_norm = 0; + return; + } + + const double rel_eps = rel_conv_tol_ * rel_conv_tol_; + const double abs_eps = abs_conv_tol_ * abs_conv_tol_; + + const auto converged_rel = [&](double r_prod) { return rel_eps > 0 && (r_prod / bi_prod) < rel_eps; }; + const auto converged_abs = [&](double r_prod) { return abs_eps > 0 && r_prod < abs_eps; }; + + // Squared norm of the recursively-updated residual as of the last PCG + // iteration, compared against the recomputed one at each restart. + double i_prod = 0.0; + bool breakdown = false; + + // Every pass through this loop's head (including the last, however + // PCG stopped) recomputes the true residual in double-double, so + // convergence is only declared on it, never on the recursive one. + for (int cycle = 0;; ++cycle) + { + auto phase_begin = clock::now(); + + dd_residual(rhs, result, r); + const double r_prod = dot(r, r); + final_res_norm = sqrt(r_prod); + + if (cycle == 0) + { + SPDLOG_INFO("[{}] [pre_loop] [{:.6f}] [rhs_norm={}] [residual={}]", name(), elapsed_seconds(phase_begin), sqrt(bi_prod), final_res_norm); + } + else + { + num_restarts = cycle; + SPDLOG_INFO("[{}] [restart] [{:.6f}] [iter={}] [residual={}] [recursive_residual={}]", name(), elapsed_seconds(phase_begin), num_iterations, final_res_norm, sqrt(i_prod)); + } + + if (converged_rel(r_prod)) + { + SPDLOG_INFO("[{}] [converged_rel] [0.000000]", name()); + break; + } + if (converged_abs(r_prod)) + { + SPDLOG_INFO("[{}] [converged_abs] [0.000000]", name()); + break; + } + if (breakdown || num_iterations >= max_iter_) + { + break; + } + + // Restart PCG from the recomputed residual: throw away the + // previous search direction and start over from p = M^-1 r. + thrust::fill(z.begin(), z.end(), 0.0); + amg_precond_iter(precond, r, z); + vector_copy(z, p); + double gamma = dot(r, z); + + for (int k = 0; k < restart_interval && num_iterations < max_iter_; ++k) + { + auto iter_begin = clock::now(); + ++num_iterations; + + matmul(p, buffer); + double sdotp = dot(p, buffer); + + if (sdotp == 0.0) + { + SPDLOG_INFO("[{}] [err_zero_sdotp] [0.000000]", name()); + breakdown = true; + break; + } + + double alpha = gamma / sdotp; + + if (alpha <= 0.0) + { + SPDLOG_INFO("[{}] [err_negative_alpha] [0.000000]", name()); + breakdown = true; + break; + } + else if (alpha < __DBL_MIN__) + { + SPDLOG_INFO("[{}] [err_subnormal_alpha] [0.000000]", name()); + breakdown = true; + break; + } + + vector_add(alpha, p, result); + vector_add(-1.0 * alpha, buffer, r); + i_prod = dot(r, r); + + // Only a hint: cut this cycle short and let the recomputed + // residual at the loop head confirm (or refute) convergence. + if (converged_rel(i_prod) || converged_abs(i_prod)) + { + break; + } + + thrust::fill(z.begin(), z.end(), 0.0); + amg_precond_iter(precond, r, z); + + double old_gamma = gamma; + gamma = dot(r, z); + + double beta = gamma / old_gamma; + + vector_scale(beta, p); + vector_add(1.0, z, p); + + CHECK_CUDA(cudaDeviceSynchronize()); + RESTARTED_AMG_PCG_LOG_ITER("[{}] [pcg_iter] [{:.6f}] [iter={}] [residual={}]", name(), elapsed_seconds(iter_begin), num_iterations - 1, sqrt(i_prod)); + } + } + } + + RestartedAMGPCG::~RestartedAMGPCG() + { + if (has_matrix_) + { + HYPRE_IJMatrixDestroy(A); + has_matrix_ = false; + A = nullptr; + } + } +} // namespace polysolve::linear diff --git a/src/polysolve/linear/RestartedAMGPCG.hpp b/src/polysolve/linear/RestartedAMGPCG.hpp new file mode 100644 index 00000000..4bc21269 --- /dev/null +++ b/src/polysolve/linear/RestartedAMGPCG.hpp @@ -0,0 +1,135 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "Solver.hpp" + +#include +#include + +#include +#include +#include + +#include + +extern "C" +{ + HYPRE_Int hypre_ParVectorAxpy(HYPRE_Complex alpha, HYPRE_ParVector x, HYPRE_ParVector y); +} + +namespace polysolve::linear +{ + + // BoomerAMG-preconditioned CG on the GPU. AMG is set up exactly as in + // GPUHybridSolver (minus its problematic-dof subspace correction). Every + // restart_interval iterations the true residual b - Ax is recomputed in + // double-double precision and CG is restarted from it, so the + // recursively-updated residual can't drift away from the true one; + // convergence is only ever declared on that recomputed residual. + class RestartedAMGPCG : public Solver + { + + public: + RestartedAMGPCG(); + ~RestartedAMGPCG(); + + private: + POLYSOLVE_DELETE_MOVE_COPY(RestartedAMGPCG) + + public: + ////////////////////// + // Public interface // + ////////////////////// + + // Set solver parameters + virtual void set_parameters(const json ¶ms) override; + + // Set block size for multigrid solvers + virtual void set_block_size(int block_size) override; + + // Set the function (block) assigned to each row for multigrid solvers + virtual void set_block_mapping(const Eigen::VectorXi &block_mapping) override; + + // Retrieve solve information + virtual void get_info(json ¶ms) const override; + + // 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 "RestartedAMGPCG"; } + + protected: + // AMG settings + double theta = 0.5; + + // Number of PCG iterations between double-double residual + // recomputations (each of which restarts PCG). + int restart_interval = 50; + + // General solver settings + int dimension_ = 1; // 1 = scalar (Laplace), 2 or 3 = vector (Elasticity) + + // Optional per-row function (block) assignment for HYPRE_BoomerAMGSetDofFunc. + // Empty means the default interleaved mapping (row i belongs to function i % dimension_). + Eigen::VectorXi block_mapping_; + int max_iter_ = 10000; + double rel_conv_tol_ = 1e-10; + double abs_conv_tol_ = 0.0; + + // When false, skips the per-PCG-iteration logs (matmul, amg_v_cycle, + // pcg_iter) so a long solve doesn't flood the log; every other log + // line, including the per-restart one, is unaffected. + bool detailed_log = true; + + // solve information + int num_iterations = 0; + int num_restarts = 0; + double final_res_norm = 0.0; + + private: + bool has_matrix_ = false; + + // Hypre variables + HYPRE_IJMatrix A; + HYPRE_ParCSRMatrix parcsr_A; + HYPRE_IJVector ij_x; + HYPRE_IJVector ij_b; + + // Device copy of the system matrix (CSR, same layout handed to + // HYPRE). Unlike GPUHybridSolver, this is kept after factorize() for + // the double-double residual, at the cost of a second copy of the + // matrix alongside HYPRE's own. + thrust::device_vector d_inner_indices; + thrust::device_vector d_outer_indices; + thrust::device_vector d_values; + + public: + // factorization helpers + void copy_matrix_to_hypre(); + + // solve helpers + void init_hypre_vectors(const int size); + + // linear algebra + void set_hypre_vec(HYPRE_IJVector &ij_x, HYPRE_ParVector &par_x, const thrust::device_vector &x); + void matmul(const thrust::device_vector &x, thrust::device_vector &result); + double dot(const thrust::device_vector &x, const thrust::device_vector &y); + void vector_copy(const thrust::device_vector &x, thrust::device_vector &y); + void vector_add(double alpha, const thrust::device_vector &x, thrust::device_vector &y); + void vector_scale(double alpha, thrust::device_vector &x); + + // r = b - Ax, accumulated in double-double precision and rounded to double + void dd_residual(const thrust::device_vector &b, const thrust::device_vector &x, thrust::device_vector &r); + + // preconditioning functions + void amg_precond_iter(const HYPRE_Solver &precond, thrust::device_vector &b, thrust::device_vector &x); + + // Krylov solve methods + void pcg_solve(thrust::device_vector &rhs, thrust::device_vector &result, HYPRE_Solver &precond); + }; + +} // namespace polysolve::linear diff --git a/src/polysolve/linear/Solver.cpp b/src/polysolve/linear/Solver.cpp index f112b702..544d71c3 100644 --- a/src/polysolve/linear/Solver.cpp +++ b/src/polysolve/linear/Solver.cpp @@ -71,6 +71,7 @@ namespace polysolve::linear #endif #ifdef POLYSOLVE_WITH_GPU_HYBRID #include "GPUHybridSolver.hpp" +#include "RestartedAMGPCG.hpp" #endif #ifdef POLYSOLVE_WITH_CUDSS #include "cuDSS.hpp" @@ -146,7 +147,7 @@ namespace polysolve::linear } } - std::unique_ptr Solver::create(const json ¶ms_in, spdlog::logger &logger, const bool strict_validation) + std::unique_ptr Solver::create(const json ¶ms_in, spdlog::logger &logger, const bool strict_validation, const int dimension) { json params = params_in; // mutable copy @@ -165,6 +166,7 @@ namespace polysolve::linear params = jse.inject_defaults(params, rules); auto res = create(params["solver"], params["precond"]); + res->set_block_size(dimension); res->set_parameters(params); return res; @@ -423,6 +425,10 @@ namespace polysolve::linear else if (solver == "GPUHybrid") { return std::make_unique(); + } + else if (solver == "RestartedAMGPCG") + { + return std::make_unique(); #endif #ifdef POLYSOLVE_WITH_CPU_HYBRID } @@ -577,6 +583,7 @@ namespace polysolve::linear #endif #ifdef POLYSOLVE_WITH_GPU_HYBRID "GPUHybrid", + "RestartedAMGPCG", #endif #ifdef POLYSOLVE_WITH_CPU_HYBRID "CPUHybrid", diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index c6837056..aad2db96 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -3,6 +3,8 @@ #include #include +#include +#include #define POLYSOLVE_DELETE_MOVE_COPY(Base) \ Base(Base &&) = delete; \ @@ -56,11 +58,13 @@ namespace polysolve::linear /// @param[in] params Parameter of the solver, including name and preconditioner /// @param[in] logger Logger used for error /// @param[in] strict_validation strict validation of the input paraams + /// @param[in] dimension block size for multigrid solvers, set before set_parameters /// @return a pointer to a linear solver // static std::unique_ptr create(const json ¶ms, spdlog::logger &logger, - const bool strict_validation = true); + const bool strict_validation = true, + const int dimension = 1); /// @brief Static constructor /// @@ -110,12 +114,24 @@ namespace polysolve::linear /// Set block size for multigrid solvers virtual void set_block_size(int block_size) {} + /// Set block mapping for multigrid solvers + virtual void set_block_mapping(const Eigen::VectorXi &block_mapping) {} + /// If the problem is nullspace for multigrid solvers virtual void set_is_nullspace(const VectorXd &x) {} /// Set solver tolerance virtual void set_tolerance(const double tol) {} + /// Set the set of DOFs flagged as problematic (e.g. poorly conditioned or + /// in contact), letting a caller drive subdomain selection instead of the + /// solver's own heuristics. + virtual void set_problematic_dofs(const std::set &bad_dofs) {} + + /// Contact patches (e.g. one per contact pair), used by solvers that + /// support Schwarz-style overlapping subdomain treatment of contact DOFs. + std::vector> contact_patches; + /// /// @brief { Solve the linear system Ax = b } /// diff --git a/src/polysolve/linear/cuDSS.cpp b/src/polysolve/linear/cuDSS.cpp index 348ff771..1d9b87fa 100644 --- a/src/polysolve/linear/cuDSS.cpp +++ b/src/polysolve/linear/cuDSS.cpp @@ -123,7 +123,7 @@ namespace polysolve::linear auto phase_begin = clock::now(); CHECK_CUDSS(cudssConfigCreate(&config)); CHECK_CUDSS(cudssDataCreate(cudss_handle, &solverData)); - SPDLOG_TRACE("[cuDSS] [create_solver_objects] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [create_solver_objects] [{:.6f}]", elapsed_seconds(phase_begin)); } { @@ -147,7 +147,7 @@ namespace polysolve::linear CUDSS_MVIEW_FULL, CUDSS_BASE_ZERO)); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[cuDSS] [copy_sparse_matrix] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [copy_sparse_matrix] [{:.6f}]", elapsed_seconds(phase_begin)); } { @@ -159,7 +159,7 @@ namespace polysolve::linear MatrixA, nullptr, nullptr)); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[cuDSS] [pattern_analysis] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [pattern_analysis] [{:.6f}]", elapsed_seconds(phase_begin)); } } @@ -169,7 +169,7 @@ namespace polysolve::linear CHECK_CUDSS(cudssExecute(cudss_handle, CUDSS_PHASE_FACTORIZATION, config, solverData, MatrixA, nullptr, nullptr)); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[cuDSS] [numerical_factorization] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [numerical_factorization] [{:.6f}]", elapsed_seconds(phase_begin)); } void cuDSSSolver::solve(const Ref b, Ref x) @@ -194,7 +194,7 @@ namespace polysolve::linear CHECK_CUDA(cudaMemcpy(d_x, x.data(), m_nrows * sizeof(double), cudaMemcpyHostToDevice)); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[cuDSS] [copy_vectors] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [copy_vectors] [{:.6f}]", elapsed_seconds(phase_begin)); } { @@ -202,7 +202,7 @@ namespace polysolve::linear CHECK_CUDSS(cudssExecute(cudss_handle, CUDSS_PHASE_SOLVE, config, solverData, MatrixA, MatrixX, MatrixB)); CHECK_CUDA(cudaDeviceSynchronize()); - SPDLOG_TRACE("[cuDSS] [solve] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[cuDSS] [solve] [{:.6f}]", elapsed_seconds(phase_begin)); } CHECK_CUDA(cudaMemcpy(x.data(), d_x, m_nrows * sizeof(double), cudaMemcpyDeviceToHost)); diff --git a/src/polysolve/linear/mas_utils/MASPreconditioner.cu b/src/polysolve/linear/mas_utils/MASPreconditioner.cu index 5b099c7f..c11f490d 100644 --- a/src/polysolve/linear/mas_utils/MASPreconditioner.cu +++ b/src/polysolve/linear/mas_utils/MASPreconditioner.cu @@ -1049,7 +1049,7 @@ namespace polysolve::linear::mas *(padded_topology_.cols)); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_build_padded_topology_device] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_build_padded_topology_device] [{:.6f}]", elapsed_seconds(phase_begin)); // Build coarse space hierarchy. phase_begin = clock::now(); @@ -1060,7 +1060,7 @@ namespace polysolve::linear::mas *(padded_topology_.padded_to_real), rt); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_build_coarse_space] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_build_coarse_space] [{:.6f}]", elapsed_seconds(phase_begin)); // Build coarse space matrices by 1. gather coarse space hessian from fine nodes 2. invert phase_begin = clock::now(); @@ -1071,7 +1071,7 @@ namespace polysolve::linear::mas padded_topology_.padded_node_num, rt); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_build_coarse_matrices] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_build_coarse_matrices] [{:.6f}]", elapsed_seconds(phase_begin)); // Allocate memory for coarse space residual (r) and preconditioned residual (z). phase_begin = clock::now(); @@ -1090,8 +1090,8 @@ namespace polysolve::linear::mas coarse_vectors_.multi_level_z = safe_alloc(total_level_scalars, 0.0, rt, "MAS coarse_vectors z"); rt.stream.sync(); - SPDLOG_TRACE("[MAS] [factorize_allocate_coarse_vectors] [{:.6f}]", elapsed_seconds(phase_begin)); - SPDLOG_TRACE("[MAS] [factorize_total] [{:.6f}]", elapsed_seconds(total_begin)); + SPDLOG_INFO("[MAS] [factorize_allocate_coarse_vectors] [{:.6f}]", elapsed_seconds(phase_begin)); + SPDLOG_INFO("[MAS] [factorize_total] [{:.6f}]", elapsed_seconds(total_begin)); initialized_ = true; } diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 92711ef7..f421015a 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -6,6 +6,7 @@ #include "PostStepData.hpp" #include +#include #include namespace polysolve::nonlinear @@ -18,6 +19,21 @@ namespace polysolve::nonlinear Linf = 2 }; + // Defined here (rather than in Solver.cpp, where it previously lived) so + // that every translation unit that converts a JSON "norm_type" string to + // a NormType sees this mapping via ADL. NLOHMANN_JSON_SERIALIZE_ENUM + // expands to template<> to_json/from_json overloads, so defining it in a + // header is safe (no ODR/duplicate-symbol issues); previously, code + // outside Solver.cpp's own translation unit silently fell back to + // nlohmann's generic enum-as-integer conversion, which rejects the + // string values ("L2", etc.) that json-specs actually declares for this + // field. + NLOHMANN_JSON_SERIALIZE_ENUM( + NormType, + {{NormType::EUCLIDEAN, "Euclidean"}, + {NormType::L2, "L2"}, + {NormType::Linf, "Linf"}}) + /// @brief Class defining optimization problem to be solved. To be defined by user code class Problem { @@ -34,6 +50,28 @@ namespace polysolve::nonlinear virtual bool is_residual() const { return false; } + /// @brief Number of scalar components per node (e.g., 1 for scalar problems, + /// 2 or 3 for vector-valued problems). Used as the block size for + /// multigrid-style linear solvers. + virtual int dimension() const { return 1; } + + /// @brief Per-row function (block) assignment for multigrid-style linear + /// solvers, passed through to Solver::set_block_mapping. An empty vector + /// (the default) means the solver should fall back to its own default + /// mapping instead of a problem-specific one. + virtual Eigen::VectorXi block_mapping() const { return Eigen::VectorXi(); } + + /// @brief Report DOFs the linear solver should treat as problematic + /// (e.g. poorly conditioned or in contact), letting solvers that support + /// externally-driven subdomain selection (e.g. AMGF) use them instead of + /// their own row-norm heuristics. Default reports none. + virtual void get_problematic_dofs(std::set &bad_dofs) {} + + /// @brief Contact patches (one set of DOFs per patch), passed through to + /// Solver::contact_patches for solvers supporting Schwarz-style handling + /// of contact DOFs. + std::vector> contact_patches; + /// @brief Initialize the problem. /// @param x0 Initial guess. virtual void init(const TVector &x0) {} diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index b7c3e8b3..d3b6d114 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -10,6 +10,7 @@ #include "descent_strategies/ADAM.hpp" #include "descent_strategies/GradientDescent.hpp" #include "descent_strategies/LBFGS.hpp" +#include "descent_strategies/NonlinearCG.hpp" #include @@ -37,47 +38,52 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) + const NormType norm_type, + const int dimension) { if (solver_name == "BFGS") { return std::make_shared( solver_params, linear_solver_params, - characteristic_length, logger); + characteristic_length, logger, dimension); } else if (solver_name == "DenseNewton" || solver_name == "dense_newton") { - return std::make_shared(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "DenseProjectedNewton") { - return std::make_shared(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "DenseRegularizedNewton") { - return std::make_shared(false, false, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(false, false, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "DenseRegularizedProjectedNewton") { - return std::make_shared(false, true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(false, true, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "Newton" || solver_name == "SparseNewton" || solver_name == "sparse_newton") { - return std::make_shared(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "ProjectedNewton") { - return std::make_shared(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "RegularizedNewton") { - return std::make_shared(true, false, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(true, false, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); } else if (solver_name == "RegularizedProjectedNewton") { - return std::make_shared(true, true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + return std::make_shared(true, true, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); + } + else if (solver_name == "NewtonCG") + { + return std::make_shared(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); } else if (solver_name == "LBFGS" || solver_name == "L-BFGS") @@ -94,6 +100,11 @@ namespace polysolve::nonlinear return std::make_shared(solver_params, false, characteristic_length, logger); } + else if (solver_name == "NonlinearCG" || solver_name == "nonlinear_cg") + { + return std::make_shared(solver_params, characteristic_length, logger); + } + else if (solver_name == "ADAM" || solver_name == "adam") { return std::make_shared(solver_params, false, characteristic_length, logger); @@ -114,11 +125,9 @@ namespace polysolve::nonlinear {FiniteDiffStrategy::DIRECTIONAL_DERIVATIVE, "DirectionalDerivative"}, {FiniteDiffStrategy::FULL_FINITE_DIFF, "FullFiniteDiff"}}) - NLOHMANN_JSON_SERIALIZE_ENUM( - NormType, - {{NormType::EUCLIDEAN, "Euclidean"}, - {NormType::L2, "L2"}, - {NormType::Linf, "Linf"}}) + // NormType's serializer now lives in Problem.hpp (next to the enum + // itself), so it's visible to every translation unit that converts a + // JSON "norm_type" field, not just this one. // Static constructor std::unique_ptr Solver::create( @@ -127,7 +136,8 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger, const bool strict_validation, - const NormType norm_type) + const NormType norm_type, + const int dimension) { json solver_params = solver_params_in; // mutable copy @@ -149,7 +159,7 @@ namespace polysolve::nonlinear for (const auto &j : solver_params["solver"]) { const std::string solver_name = j["type"]; - solver->add_strategy(create_solver(solver_name, j, linear_solver_params, characteristic_length, logger, norm_type)); + solver->add_strategy(create_solver(solver_name, j, linear_solver_params, characteristic_length, logger, norm_type, dimension)); } } else @@ -158,19 +168,19 @@ namespace polysolve::nonlinear if (solver_name == "DenseNewton" || solver_name == "dense_newton") { - auto n = Newton::create_solver(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + auto n = Newton::create_solver(false, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); for (auto &s : n) solver->add_strategy(s); } else if (solver_name == "Newton" || solver_name == "SparseNewton" || solver_name == "sparse_newton") { - auto n = Newton::create_solver(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type); + auto n = Newton::create_solver(true, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension); for (auto &s : n) solver->add_strategy(s); } else { - solver->add_strategy(create_solver(solver_name, solver_params, linear_solver_params, characteristic_length, logger, norm_type)); + solver->add_strategy(create_solver(solver_name, solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension)); } if (solver_name != "GradientDescent" && solver_name != "gradient_descent") @@ -194,6 +204,7 @@ namespace polysolve::nonlinear "StochasticADAM", "GradientDescent", "StochasticGradientDescent", + "NonlinearCG", "L-BFGS"}; } diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index 1ff49771..b0d88c88 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -40,7 +40,8 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger, const bool strict_validation = true, - const NormType norm_type = NormType::EUCLIDEAN); + const NormType norm_type = NormType::EUCLIDEAN, + const int dimension = 1); /// @brief List available solvers static std::vector available_solvers(); diff --git a/src/polysolve/nonlinear/descent_strategies/BFGS.cpp b/src/polysolve/nonlinear/descent_strategies/BFGS.cpp index 4dccd3ff..65a062b8 100644 --- a/src/polysolve/nonlinear/descent_strategies/BFGS.cpp +++ b/src/polysolve/nonlinear/descent_strategies/BFGS.cpp @@ -8,10 +8,11 @@ namespace polysolve::nonlinear BFGS::BFGS(const json &solver_params, const json &linear_solver_params, const double characteristic_length, - spdlog::logger &logger) + spdlog::logger &logger, + const int dimension) : Superclass(solver_params, characteristic_length, logger) { - linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger); + linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger, true, dimension); if (!linear_solver->is_dense()) log_and_throw_error(logger, "BFGS linear solver must be dense, instead got {}", linear_solver->name()); } @@ -42,6 +43,10 @@ namespace polysolve::nonlinear } else { + const Eigen::VectorXi block_mapping = objFunc.block_mapping(); + if (block_mapping.size() > 0) + linear_solver->set_block_mapping(block_mapping); + try { linear_solver->analyze_pattern_dense(hess, hess.rows()); diff --git a/src/polysolve/nonlinear/descent_strategies/BFGS.hpp b/src/polysolve/nonlinear/descent_strategies/BFGS.hpp index e9083efc..c2e76f02 100644 --- a/src/polysolve/nonlinear/descent_strategies/BFGS.hpp +++ b/src/polysolve/nonlinear/descent_strategies/BFGS.hpp @@ -19,7 +19,8 @@ namespace polysolve::nonlinear BFGS(const json &solver_params, const json &linear_solver_params, const double characteristic_length, - spdlog::logger &logger); + spdlog::logger &logger, + const int dimension = 1); std::string name() const override { return "BFGS"; } diff --git a/src/polysolve/nonlinear/descent_strategies/CMakeLists.txt b/src/polysolve/nonlinear/descent_strategies/CMakeLists.txt index fc7b1e3f..4179c765 100644 --- a/src/polysolve/nonlinear/descent_strategies/CMakeLists.txt +++ b/src/polysolve/nonlinear/descent_strategies/CMakeLists.txt @@ -10,6 +10,8 @@ set(SOURCES ADAM.hpp Newton.hpp Newton.cpp + NonlinearCG.hpp + NonlinearCG.cpp ) source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES}) diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index ff8d015a..d1d25a74 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.cpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.cpp @@ -8,6 +8,11 @@ #include #endif +#include +#include +#include +#include + namespace polysolve::nonlinear { @@ -17,7 +22,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) + const NormType norm_type, + const int dimension) { // Copies stuff from main newton json proj_solver_params = R"({"ProjectedNewton": {}})"_json; @@ -35,21 +41,21 @@ namespace polysolve::nonlinear res.push_back(std::make_unique( sparse, solver_params, linear_solver_params, - characteristic_length, logger, norm_type)); + characteristic_length, logger, norm_type, dimension)); const bool use_psd_projection = solver_params["Newton"]["use_psd_projection"]; if (use_psd_projection) res.push_back(std::make_unique( sparse, proj_solver_params, linear_solver_params, - characteristic_length, logger, norm_type)); + characteristic_length, logger, norm_type, dimension)); const double reg_weight_min = solver_params["Newton"]["reg_weight_min"]; if (reg_weight_min > 0) res.push_back(std::make_unique( sparse, solver_params["Newton"]["use_psd_projection_in_regularized"], reg_solver_params, linear_solver_params, - characteristic_length, logger, norm_type)); + characteristic_length, logger, norm_type, dimension)); if (res.empty()) log_and_throw_error(logger, "Newton needs to have at least one of force_psd_projection=false, reg_weight_min>0, or use_psd_projection=true"); @@ -63,11 +69,12 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) + const NormType norm_type, + const int dimension) : Superclass(solver_params, characteristic_length, logger), is_sparse(sparse), characteristic_length(characteristic_length), residual_tolerance(residual_tolerance), norm_type(norm_type) { - linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger); + linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger, true, dimension); if (linear_solver->is_dense() == sparse) log_and_throw_error(logger, "Newton linear solver must be {}, instead got {}", sparse ? "sparse" : "dense", linear_solver->name()); @@ -82,8 +89,9 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) - : Newton(sparse, extract_param("Newton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type) + const NormType norm_type, + const int dimension) + : Newton(sparse, extract_param("Newton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension) { } @@ -93,8 +101,9 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) - : Superclass(sparse, extract_param("ProjectedNewton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type) + const NormType norm_type, + const int dimension) + : Superclass(sparse, extract_param("ProjectedNewton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension) { } @@ -105,8 +114,9 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type) - : Superclass(sparse, extract_param("RegularizedNewton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type), + const NormType norm_type, + const int dimension) + : Superclass(sparse, extract_param("RegularizedNewton", "residual_tolerance", solver_params), solver_params, linear_solver_params, characteristic_length, logger, norm_type, dimension), project_to_psd(project_to_psd) { reg_weight_min = extract_param("RegularizedNewton", "reg_weight_min", solver_params); @@ -125,6 +135,28 @@ namespace polysolve::nonlinear log_and_throw_error(logger, "Newton reg_weight_max must be > {}, instead got {}", reg_weight_min, reg_weight_max); } + NewtonCG::NewtonCG( + const bool sparse, + const json &solver_params, + const json &linear_solver_params, + const double characteristic_length, + spdlog::logger &logger, + const NormType norm_type) + : Superclass(sparse, std::numeric_limits::infinity(), solver_params, linear_solver_params, characteristic_length, logger, norm_type), + force_psd_projection(solver_params.contains("NewtonCG") + ? solver_params["NewtonCG"].value("force_psd_projection", false) + : solver_params.value("force_psd_projection", false)) + { + const std::string solver_name = linear_solver_name(); + // CPUHybridSolver/GPUHybridSolver report themselves as "CPUAMGF"/"GPUAMGF" + // (rather than "CPUHybrid"/"GPUHybrid") when select_bad_dofs_from_l1_norm + // is false, i.e. the informed/AMGF-Schwarz subdomain-selection mode -- + // same underlying class, just a different display name. + if (solver_name != "CPUHybrid" && solver_name != "GPUHybrid" + && solver_name != "CPUAMGF" && solver_name != "GPUAMGF") + log_and_throw_error(logger, "NewtonCG requires the linear solver to be CPUHybridSolver or GPUHybridSolver, instead got {}", solver_name); + } + // ======================================================================= void Newton::reset(const int ndof) @@ -168,6 +200,45 @@ namespace polysolve::nonlinear return true; } + void Newton::set_linear_solver_relative_tolerance(const double tol) + { + json params; + params[linear_solver->name()]["relative_tolerance"] = tol; + linear_solver->set_parameters(params); + } + + // ======================================================================= + + bool NewtonCG::compute_update_direction( + Problem &objFunc, + const TVector &x, + const TVector &grad, + TVector &direction) + { + const double grad_norm = objFunc.grad_norm(grad, get_norm_type()); + set_linear_solver_relative_tolerance(std::min(0.5, std::sqrt(grad_norm))); + + return Superclass::compute_update_direction(objFunc, x, grad, direction); + } + + void NewtonCG::compute_hessian(Problem &objFunc, + const TVector &x, + polysolve::StiffnessMatrix &hessian) + + { + objFunc.set_project_to_psd(force_psd_projection); + objFunc.hessian(x, hessian); + } + + void NewtonCG::compute_hessian(Problem &objFunc, + const TVector &x, + Eigen::MatrixXd &hessian) + + { + objFunc.set_project_to_psd(force_psd_projection); + objFunc.hessian(x, hessian); + } + // ======================================================================= double Newton::solve_sparse_linear_system(Problem &objFunc, @@ -185,6 +256,15 @@ namespace polysolve::nonlinear { POLYSOLVE_SCOPED_STOPWATCH("linear solve", this->inverting_time, m_logger); + const Eigen::VectorXi block_mapping = objFunc.block_mapping(); + if (block_mapping.size() > 0) + linear_solver->set_block_mapping(block_mapping); + + std::set bad_dofs; + objFunc.get_problematic_dofs(bad_dofs); + linear_solver->set_problematic_dofs(bad_dofs); + linear_solver->contact_patches = objFunc.contact_patches; + // TODO: get the correct size linear_solver->analyze_pattern(hessian, hessian.rows()); @@ -228,6 +308,10 @@ namespace polysolve::nonlinear { POLYSOLVE_SCOPED_STOPWATCH("linear solve", this->inverting_time, m_logger); + const Eigen::VectorXi block_mapping = objFunc.block_mapping(); + if (block_mapping.size() > 0) + linear_solver->set_block_mapping(block_mapping); + try { linear_solver->analyze_pattern_dense(hessian, hessian.rows()); diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.hpp b/src/polysolve/nonlinear/descent_strategies/Newton.hpp index 507f280f..fb2c0bee 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.hpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.hpp @@ -18,7 +18,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type); + const NormType norm_type, + const int dimension = 1); protected: Newton(const bool sparse, @@ -27,7 +28,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type); + const NormType norm_type, + const int dimension = 1); public: Newton(const bool sparse, @@ -35,7 +37,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type); + const NormType norm_type, + const int dimension = 1); std::string name() const override { return internal_name() + "Newton"; } @@ -61,6 +64,10 @@ namespace polysolve::nonlinear protected: std::string internal_name() const { return is_sparse ? "Sparse" : "Dense"; } + NormType get_norm_type() const { return norm_type; } + + std::string linear_solver_name() const { return linear_solver->name(); } + void set_linear_solver_relative_tolerance(const double tol); virtual void compute_hessian(Problem &objFunc, const TVector &x, @@ -89,7 +96,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type); + const NormType norm_type, + const int dimension = 1); std::string name() const override { return internal_name() + "ProjectedNewton"; } @@ -113,7 +121,8 @@ namespace polysolve::nonlinear const json &linear_solver_params, const double characteristic_length, spdlog::logger &logger, - const NormType norm_type); + const NormType norm_type, + const int dimension = 1); std::string name() const override { @@ -143,4 +152,37 @@ namespace polysolve::nonlinear Eigen::MatrixXd &hessian) override; }; + /// @brief Inexact Newton solver that requires a Hybrid CG-based linear + /// solver (CPUHybridSolver or GPUHybridSolver). The relative tolerance of + /// the linear solve is adapted to the current gradient norm, and the + /// resulting (possibly inexact) direction is always accepted. + class NewtonCG : public Newton + { + public: + using Superclass = Newton; + + NewtonCG(const bool sparse, + const json &solver_params, + const json &linear_solver_params, + const double characteristic_length, + spdlog::logger &logger, + const NormType norm_type); + + std::string name() const override { return internal_name() + "NewtonCG"; } + + bool compute_update_direction(Problem &objFunc, const TVector &x, const TVector &grad, TVector &direction) override; + + private: + const bool force_psd_projection; + + protected: + void compute_hessian(Problem &objFunc, + const TVector &x, + polysolve::StiffnessMatrix &hessian) override; + + void compute_hessian(Problem &objFunc, + const TVector &x, + Eigen::MatrixXd &hessian) override; + }; + } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/descent_strategies/NonlinearCG.cpp b/src/polysolve/nonlinear/descent_strategies/NonlinearCG.cpp new file mode 100644 index 00000000..fd25e3d5 --- /dev/null +++ b/src/polysolve/nonlinear/descent_strategies/NonlinearCG.cpp @@ -0,0 +1,98 @@ +#include "NonlinearCG.hpp" + +#include + +#include + +namespace polysolve::nonlinear +{ + NonlinearCG::NonlinearCG(const json &solver_params, + const double characteristic_length, + spdlog::logger &logger) + : Superclass(solver_params, characteristic_length, logger) + { + const json ¶ms = solver_params.find("NonlinearCG") != solver_params.end() + ? solver_params["NonlinearCG"] + : solver_params; + + const std::string formula = params["formula"]; + if (formula == "FletcherReeves") + beta_formula_ = BetaFormula::FletcherReeves; + else if (formula == "PolakRibiere") + beta_formula_ = BetaFormula::PolakRibiere; + else if (formula == "HestenesStiefel") + beta_formula_ = BetaFormula::HestenesStiefel; + else + log_and_throw_error(logger, "Unknown NonlinearCG formula: {}", formula); + + restart_frequency_ = params["restart_frequency"]; + if (restart_frequency_ < 0) + log_and_throw_error(logger, "NonlinearCG restart_frequency must be >= 0, instead got {}", restart_frequency_); + } + + void NonlinearCG::reset(const int ndof) + { + Superclass::reset(ndof); + ndof_ = ndof; + prev_direction_.resize(0); + iters_since_restart_ = 0; + } + + bool NonlinearCG::compute_update_direction( + Problem &objFunc, + const TVector &x, + const TVector &grad, + TVector &direction) + { + const int effective_restart_frequency = restart_frequency_ > 0 ? restart_frequency_ : ndof_; + const bool restart = prev_direction_.size() == 0 + || (effective_restart_frequency > 0 && iters_since_restart_ >= effective_restart_frequency); + + if (restart) + { + direction = -grad; + } + else + { + const double prev_grad_sq_norm = prev_grad_.squaredNorm(); + double beta = 0; + + if (prev_grad_sq_norm > 0) + { + switch (beta_formula_) + { + case BetaFormula::FletcherReeves: + beta = grad.squaredNorm() / prev_grad_sq_norm; + break; + case BetaFormula::PolakRibiere: + beta = grad.dot(grad - prev_grad_) / prev_grad_sq_norm; + break; + case BetaFormula::HestenesStiefel: + { + const TVector y = grad - prev_grad_; + const double denom = prev_direction_.dot(y); + beta = denom != 0 ? grad.dot(y) / denom : 0; + break; + } + default: + log_and_throw_error(m_logger, "Unknown NonlinearCG beta formula"); + } + } + + // Negative beta triggers an automatic restart to steepest descent. + beta = std::max(beta, 0.0); + + direction = -grad + beta * prev_direction_; + + // Guard against loss of descent from numerical error. + if (grad.dot(direction) >= 0) + direction = -grad; + } + + prev_grad_ = grad; + prev_direction_ = direction; + iters_since_restart_ = restart ? 1 : iters_since_restart_ + 1; + + return true; + } +} // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/descent_strategies/NonlinearCG.hpp b/src/polysolve/nonlinear/descent_strategies/NonlinearCG.hpp new file mode 100644 index 00000000..4251617e --- /dev/null +++ b/src/polysolve/nonlinear/descent_strategies/NonlinearCG.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "DescentStrategy.hpp" +#include + +namespace polysolve::nonlinear +{ + /// @brief Unpreconditioned nonlinear conjugate gradient descent strategy. + /// + /// Computes the search direction as d_k = -g_k + beta_k * d_{k-1}, where + /// beta_k is computed using the Fletcher-Reeves, Polak-Ribiere, or + /// Hestenes-Stiefel formula (clipped to be non-negative, so a negative + /// beta triggers an automatic restart to steepest descent). The direction + /// is also periodically restarted to steepest descent every + /// restart_frequency iterations. + class NonlinearCG : public DescentStrategy + { + public: + using Superclass = DescentStrategy; + + enum class BetaFormula + { + FletcherReeves, + PolakRibiere, + HestenesStiefel + }; + + NonlinearCG(const json &solver_params, + const double characteristic_length, + spdlog::logger &logger); + + std::string name() const override { return "NonlinearCG"; } + + void reset(const int ndof) override; + + bool compute_update_direction( + Problem &objFunc, + const TVector &x, + const TVector &grad, + TVector &direction) override; + + private: + BetaFormula beta_formula_; + /// Number of iterations between forced restarts to steepest descent. + /// A value of 0 defaults to the number of degrees of freedom. + int restart_frequency_; + + int ndof_ = 0; + int iters_since_restart_ = 0; + + TVector prev_grad_; + TVector prev_direction_; + }; +} // namespace polysolve::nonlinear diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b4c33102..6fd03020 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -10,6 +10,14 @@ set(test_sources ) add_executable(unit_tests ${test_sources}) +# Standalone CLI driver: load a Matrix Market matrix/rhs, run any solver from +# a JSON config, and report timing/residual -- for reproducing/benchmarking a +# solve outside of the full polyfem integration. +# Uses POSIX-only APIs (unistd.h, sys/resource.h), so it can't build on Windows. +if(NOT WIN32) + add_executable(linear_solve linear_solve.cpp) +endif() + ################################################################################ # Required Libraries ################################################################################ @@ -19,11 +27,28 @@ target_link_libraries(unit_tests PUBLIC Catch2::Catch2) target_link_libraries(unit_tests PUBLIC polysolve::polysolve) +include(argparse) +if(TARGET linear_solve) + target_link_libraries(linear_solve PUBLIC polysolve::polysolve) + if(TARGET argparse::argparse) + target_link_libraries(linear_solve PRIVATE argparse::argparse) + endif() + if(DEFINED argparse_SOURCE_DIR) + target_include_directories(linear_solve SYSTEM PRIVATE "${argparse_SOURCE_DIR}/include") + endif() +endif() + include(polysolve_warnings) target_link_libraries(unit_tests PRIVATE polysolve::warnings) +if(TARGET linear_solve) + target_link_libraries(linear_solve PRIVATE polysolve::warnings) +endif() if(POLYSOLVE_WITH_SANITIZERS) add_sanitizers(unit_tests) + if(TARGET linear_solve) + add_sanitizers(linear_solve) + endif() endif() if(POLYSOLVE_WITH_AMGCL) @@ -40,6 +65,7 @@ target_link_libraries(unit_tests PRIVATE polyfem::data) foreach(source IN ITEMS ${test_sources}) source_group("tests" FILES "${source}") endforeach() +source_group("tests" FILES "linear_solve.cpp") # Register tests set(PARSE_CATCH_TESTS_ADD_TO_CONFIGURE_DEPENDS ON) diff --git a/tests/linear_solve.cpp b/tests/linear_solve.cpp new file mode 100644 index 00000000..4e7460ec --- /dev/null +++ b/tests/linear_solve.cpp @@ -0,0 +1,453 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace polysolve; +using namespace polysolve::linear; + +namespace +{ + constexpr std::uint32_t DEFAULT_RAND_SEED = 0; + + // Silences a solver's own logging during warmup iterations, so only the + // timed (post-warmup) iterations' [start_analyze_pattern]/[start_factorize]/ + // [start_solve]/[residual=...] lines show up. + class ScopedOutputSilencer + { + public: + ScopedOutputSilencer() + { + flush_output(); + + null_fd_ = open("/dev/null", O_WRONLY); + if (null_fd_ == -1) + { + throw std::runtime_error("failed to open /dev/null"); + } + + stdout_fd_ = dup(STDOUT_FILENO); + if (stdout_fd_ == -1) + { + close(null_fd_); + throw std::runtime_error("failed to duplicate stdout"); + } + + stderr_fd_ = dup(STDERR_FILENO); + if (stderr_fd_ == -1) + { + close(stdout_fd_); + close(null_fd_); + throw std::runtime_error("failed to duplicate stderr"); + } + + if (dup2(null_fd_, STDOUT_FILENO) == -1 || dup2(null_fd_, STDERR_FILENO) == -1) + { + restore(); + throw std::runtime_error("failed to redirect output to /dev/null"); + } + } + + ~ScopedOutputSilencer() + { + restore(); + } + + ScopedOutputSilencer(const ScopedOutputSilencer &) = delete; + ScopedOutputSilencer &operator=(const ScopedOutputSilencer &) = delete; + + private: + void restore() noexcept + { + flush_output(); + + if (stdout_fd_ != -1) + { + (void)dup2(stdout_fd_, STDOUT_FILENO); + close(stdout_fd_); + stdout_fd_ = -1; + } + if (stderr_fd_ != -1) + { + (void)dup2(stderr_fd_, STDERR_FILENO); + close(stderr_fd_); + stderr_fd_ = -1; + } + if (null_fd_ != -1) + { + close(null_fd_); + null_fd_ = -1; + } + } + + static void flush_output() noexcept + { + std::fflush(stdout); + std::fflush(stderr); + } + + int null_fd_ = -1; + int stdout_fd_ = -1; + int stderr_fd_ = -1; + }; + + Eigen::VectorXd make_random_rhs(const Eigen::Index size, const std::optional seed) + { + std::mt19937 generator; + if (seed.has_value()) + { + generator.seed(*seed); + } + else + { + std::random_device random_device; + generator.seed(random_device()); + } + + std::uniform_real_distribution distribution(-1.0, 1.0); + Eigen::VectorXd rhs(size); + for (Eigen::Index i = 0; i < size; ++i) + { + rhs[i] = distribution(generator); + } + return rhs; + } + + // A * x with every product and partial sum's rounding error recovered + // exactly (TwoProduct via std::fma, Knuth's TwoSum; Ogita, Rump & Oishi's + // Dot2) and accumulated separately, so each entry is as accurate as if + // computed in twice double precision and rounded once. Used for b = A * + // x_true, to keep x_true as close as possible to the exact solution of + // Ax = b. (Assumes no FP contraction of the plain + and * below into an + // FMA, which holds unless built with -mfma/-march=native.) + Eigen::VectorXd multiply_accurately(const Eigen::SparseMatrix &A, const Eigen::VectorXd &x) + { + Eigen::VectorXd hi = Eigen::VectorXd::Zero(A.rows()); + Eigen::VectorXd lo = Eigen::VectorXd::Zero(A.rows()); + for (int k = 0; k < A.outerSize(); ++k) + { + for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) + { + const double prod = it.value() * x[it.col()]; + const double prod_err = std::fma(it.value(), x[it.col()], -prod); + + const double old_hi = hi[it.row()]; + const double sum = old_hi + prod; + const double bb = sum - old_hi; + const double sum_err = (old_hi - (sum - bb)) + (prod - bb); + + hi[it.row()] = sum; + lo[it.row()] += sum_err + prod_err; + } + } + return hi + lo; + } + + size_t getPeakRSS() + { + struct rusage rusage; + getrusage(RUSAGE_SELF, &rusage); + return (size_t)(rusage.ru_maxrss * 1024L); + } +} // namespace + +int main(int argc, char *argv[]) +{ + argparse::ArgumentParser program("linear_solve"); + program.add_description("Run a PolySolve linear solve from Matrix Market inputs."); + program.add_argument("-A") + .required() + .metavar("matrix.mtx") + .help("Input Matrix Market sparse matrix."); + auto &rhs_group = program.add_mutually_exclusive_group(); + rhs_group.add_argument("-b") + .metavar("rhs.mtx") + .help("Optional Matrix Market RHS vector. Defaults to a deterministic random vector with seed 0 unless --rand is used."); + rhs_group.add_argument("--rand") + .metavar("seed") + .scan<'i', long long>() + .nargs(argparse::nargs_pattern::optional) + .help("Generate a random RHS instead of loading -b. Optionally provide a seed."); + rhs_group.add_argument("--x_true") + .metavar("x_true.mtx") + .help("Matrix Market vector used as the known solution: the RHS is set to A * x_true " + "(accumulated in twice double precision, rounded once) instead of loading -b, " + "and each solve's error against x_true is reported."); + program.add_argument("-o", "--output") + .metavar("x.mtx") + .help("Optional Matrix Market vector file to write the solution of the last (timed) solve to."); + program.add_argument("-i") + .metavar("problem_info.mtx") + .help("Optional Matrix Market vector used to derive the problematic DOFs passed to " + "Solver::set_problematic_dofs (e.g. Hybrid's AMGF mode)."); + program.add_argument("-t", "--threshold") + .help("Threshold value to apply against -i.") + .scan<'g', double>() + .default_value(0.0); + program.add_argument("--less_than") + .default_value(false) + .implicit_value(true) + .help("Condition to select problematic DoFs relative to threshold. (Defaults to problematic_info(i) > thresh.)"); + program.add_argument("-j") + .metavar("spec.json") + .help("Optional solver JSON config. Defaults to {\"solver\":\"Eigen::SimplicialLDLT\"}."); + program.add_argument("--block_size") + .scan<'i', int>() + .metavar("block_size") + .help("Optional block size passed to Solver::set_block_size (e.g. 3 for 3D elasticity). " + "Defaults to the solver's own."); + program.add_argument("-w") + .default_value(1) + .scan<'i', int>() + .metavar("warmup") + .help("Number of warmup iterations. Default: 1."); + program.add_argument("-r") + .default_value(1) + .scan<'i', int>() + .metavar("repeat") + .help("Number of iterations after warmup. Default: 1."); + program.add_argument("--force_symmetry") + .default_value(false) + .implicit_value(true) + .help("Force matrix to be symmetric. Default: false."); + + program.parse_args(argc, argv); + + const std::string matrix_path = program.get("-A"); + const int warmup = program.get("-w"); + const int repeat = program.get("-r"); + const bool use_random_rhs = program.is_used("--rand"); + const auto rand_seed_argument = program.present("--rand"); + const bool force_symmetry = program.get("--force_symmetry"); + + if (warmup < 0) + { + throw std::runtime_error("warmup count must be non-negative"); + } + if (repeat < 1) + { + throw std::runtime_error("repeat count must be at least 1"); + } + + std::optional random_seed; + if (rand_seed_argument.has_value()) + { + if (*rand_seed_argument < 0) + { + throw std::runtime_error("random seed must be non-negative"); + } + if (static_cast(*rand_seed_argument) > std::numeric_limits::max()) + { + throw std::runtime_error("random seed exceeds uint32_t range"); + } + random_seed = static_cast(*rand_seed_argument); + } + + json solver_config = json::object({{"solver", "Eigen::SimplicialLDLT"}}); + if (program.present("-j")) + { + const std::string json_path = program.get("-j"); + std::ifstream in(json_path); + if (!in.is_open()) + { + throw std::runtime_error("failed to open json config: " + json_path); + } + in >> solver_config; + } + + auto logger = spdlog::stderr_color_mt("linear_solve"); + logger->set_level(spdlog::level::off); + + // Validate solver configuration before loading potentially large inputs. + { + auto solver = Solver::create(solver_config, *logger); + (void)solver; + } + + Eigen::SparseMatrix A; + if (!loadMarket(A, matrix_path)) + { + throw std::runtime_error("failed to load matrix market matrix: " + matrix_path); + } + if (A.rows() != A.cols()) + { + throw std::runtime_error("matrix must be square"); + } + if (force_symmetry) + { + std::vector> triplets; + triplets.reserve(2 * (A.nonZeros() - A.rows()) + A.rows()); + + for (int k = 0; k < A.outerSize(); ++k) + { + for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) + { + triplets.push_back(Eigen::Triplet(it.row(), it.col(), it.value())); + if (it.col() != it.row()) + { + triplets.push_back(Eigen::Triplet(it.col(), it.row(), it.value())); + } + } + } + + A.setFromTriplets(triplets.begin(), triplets.end()); + } + A.makeCompressed(); + + Eigen::VectorXd b(A.rows()); + std::optional x_true; + if (program.present("--x_true")) + { + const std::string x_true_path = program.get("--x_true"); + x_true.emplace(A.rows()); + if (!Eigen::loadMarketVector(*x_true, x_true_path)) + { + throw std::runtime_error("failed to load matrix market x_true: " + x_true_path); + } + if (x_true->size() != A.rows()) + { + throw std::runtime_error("x_true dimension mismatch"); + } + b = multiply_accurately(A, *x_true); + } + else if (program.present("-b")) + { + const std::string rhs_path = program.get("-b"); + if (!Eigen::loadMarketVector(b, rhs_path)) + { + throw std::runtime_error("failed to load matrix market rhs: " + rhs_path); + } + if (b.size() != A.rows()) + { + throw std::runtime_error("rhs dimension mismatch"); + } + } + else + { + b = make_random_rhs(A.rows(), use_random_rhs ? random_seed : std::optional(DEFAULT_RAND_SEED)); + } + + std::set problem_specific_bad_dofs; + if (program.present("-i")) + { + const bool less_than = program.get("--less_than"); + const double problematic_threshold = program.get("-t"); + const std::string p_path = program.get("-i"); + Eigen::VectorXd problematic_info(A.rows()); + + if (!Eigen::loadMarketVector(problematic_info, p_path)) + { + throw std::runtime_error("failed to load matrix market problem info: " + p_path); + } + if (problematic_info.size() != A.rows()) + { + throw std::runtime_error("problematic info dimension mismatch"); + } + + for (int i = 0; i < problematic_info.size(); ++i) + { + const bool is_bad = less_than + ? problematic_info(i) < problematic_threshold + : problematic_info(i) > problematic_threshold; + if (is_bad) + { + problem_specific_bad_dofs.insert(i); + } + } + } + + std::optional dense_A; + Eigen::VectorXd x; + + // Created once and reused for every iteration (rather than once per + // iteration) so that solvers with a process-wide setup/teardown cost -- + // e.g. CPUHybrid, whose construction/destruction respectively spins up + // and fully tears down a multi-thread nanompi rank team -- don't pay + // that cost (and, for CPUHybrid, don't race tearing the team down and + // immediately rebuilding it) on every -w/-r iteration. + auto solver = Solver::create(solver_config, *logger); + if (const auto block_size = program.present("--block_size")) + { + solver->set_block_size(*block_size); + } + solver->set_problematic_dofs(problem_specific_bad_dofs); + + const int iterations = warmup + repeat; + for (int iteration = 0; iteration < iterations; ++iteration) + { + const std::unique_ptr silencer = + iteration < warmup ? std::make_unique() : nullptr; + + SPDLOG_INFO("[{}] [matrix_info] [size={}] [nnzs={}]", solver->name(), A.rows(), A.nonZeros()); + x.resize(A.cols()); + x.setZero(); + + if (solver->is_dense()) + { + if (!dense_A.has_value()) + { + dense_A.emplace(A); + } + solver->analyze_pattern_dense(*dense_A, A.rows()); + } + else + { + SPDLOG_INFO("[{}] [start_analyze_pattern]", solver->name()); + solver->analyze_pattern(A, A.rows()); + } + + if (solver->is_dense()) + { + solver->factorize_dense(*dense_A); + } + else + { + SPDLOG_INFO("[{}] [start_factorize]", solver->name()); + solver->factorize(A); + } + + SPDLOG_INFO("[{}] [start_solve]", solver->name()); + solver->solve(b, x); + double residual = (b - A * x).norm(); + SPDLOG_INFO("[{}] [residual={}] [peak_memory={}]", solver->name(), residual, getPeakRSS()); + + if (x_true.has_value()) + { + // The energy norm is the one CG minimizes the error in. + const Eigen::VectorXd error = x - *x_true; + const double rel_error = error.norm() / x_true->norm(); + const double rel_error_inf = error.lpNorm() / x_true->lpNorm(); + const double rel_error_energy = std::sqrt(error.dot(A * error) / x_true->dot(A * *x_true)); + SPDLOG_INFO("[{}] [x_error] [rel_error={}] [rel_error_inf={}] [rel_error_energy={}]", solver->name(), rel_error, rel_error_inf, rel_error_energy); + } + } + + if (program.present("-o")) + { + const std::string output_path = program.get("-o"); + if (!Eigen::saveMarketVector(x, output_path)) + { + throw std::runtime_error("failed to write matrix market solution: " + output_path); + } + } + + return 0; +} diff --git a/tests/test_linear_solver.cpp b/tests/test_linear_solver.cpp index 318b6e89..2b37f207 100644 --- a/tests/test_linear_solver.cpp +++ b/tests/test_linear_solver.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -134,7 +135,7 @@ TEST_CASE("all", "[solver]") params[s]["use_preconditioned_residual_norm"] = false; solver->set_parameters(params); } - else if (s == "GPUHybrid" || s == "CPUHybrid") + else if (s == "GPUHybrid" || s == "CPUHybrid" || s == "RestartedAMGPCG") { params[s]["relative_tolerance"] = 0.0; params[s]["absolute_tolerance"] = 1e-8; @@ -223,10 +224,10 @@ TEST_CASE("mas_block_dim", "[.][solver]") { auto solver = Solver::create("MAS", ""); json params; - params["MAS"]["block_dim"] = block_dim; params["MAS"]["relative_tolerance"] = 0.0; params["MAS"]["absolute_tolerance"] = 1e-6; params["MAS"]["use_preconditioned_residual_norm"] = false; + solver->set_block_size(block_dim); solver->set_parameters(params); // Eigen::VectorXd b(A.rows()); @@ -264,10 +265,10 @@ TEST_CASE("hybrid_convergence", "[.][solver]") { auto solver = Solver::create(s, ""); json params; - params[s]["block_size"] = 3; params[s]["relative_tolerance"] = 0; params[s]["absolute_tolerance"] = 1e-12; + solver->set_block_size(3); solver->set_parameters(params); Eigen::VectorXd b(A.rows()); b.setRandom(); @@ -289,6 +290,313 @@ TEST_CASE("hybrid_convergence", "[.][solver]") } } +// Contact patches may overlap (e.g. a dof shared by two contact pairs). With +// contact_patch_schwarz, each patch gets its own subdomain solve and their +// corrections are summed for shared dofs, damped by the average overlap -- +// see CPUHybridSolver::dss_precond_iter / GPUHybridSolver::dense_contact_patch_precond_iter. +// This checks that the solve still converges to the right answer with +// genuinely overlapping patches, and (GPU only) that a patch bigger than the +// fixed-size dense batched solve is rejected rather than silently corrupting +// memory. +TEST_CASE("contact_patch_schwarz", "[.][solver]") +{ + const std::string path = POLYFEM_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_contact.mtx"); + REQUIRE(ok); + + std::vector> patches; + // Two overlapping patches (dofs 3,4,5 are shared by both). + patches.push_back({0, 1, 2, 3, 4, 5}); + patches.push_back({3, 4, 5, 6, 7, 8}); + // A third, disjoint patch elsewhere in the system. + patches.push_back({30, 31, 32, 33, 34, 35}); + + std::vector solvers; +#ifdef POLYSOLVE_WITH_CPU_HYBRID + solvers.push_back("CPUHybrid"); +#endif +#ifdef POLYSOLVE_WITH_GPU_HYBRID + solvers.push_back("GPUHybrid"); +#endif + + for (const auto &s : solvers) + { + auto solver = Solver::create(s, ""); + json params; + params[s]["relative_tolerance"] = 0; + params[s]["absolute_tolerance"] = 1e-12; + params[s]["select_bad_dofs_from_l1_norm"] = false; + params[s]["decompose_subdomains"] = false; + params[s]["contact_patch_schwarz"] = true; + + solver->set_block_size(3); + solver->set_parameters(params); + solver->contact_patches = patches; + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyze_pattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + // A_contact.mtx is quite ill-conditioned, and here only a handful of + // dofs (the contact patches) get any subspace correction at all + // (select_bad_dofs_from_l1_norm is off), so convergence is slow; + // loosen the tolerance accordingly relative to the other hybrid tests. + const double err = (A * x - b).norm(); + INFO("solver: " + s); + REQUIRE(err < 1e-3); + } + +#ifdef POLYSOLVE_WITH_GPU_HYBRID + { + auto solver = Solver::create("GPUHybrid", ""); + json params; + params["GPUHybrid"]["select_bad_dofs_from_l1_norm"] = false; + params["GPUHybrid"]["decompose_subdomains"] = false; + params["GPUHybrid"]["contact_patch_schwarz"] = true; + + solver->set_block_size(3); + solver->set_parameters(params); + + std::set oversized_patch; + for (int i = 0; i < 65; ++i) + oversized_patch.insert(i); + solver->contact_patches = {oversized_patch}; + + solver->analyze_pattern(A, A.rows()); + REQUIRE_THROWS(solver->factorize(A)); + } +#endif +} + +// Checks set_block_mapping (HYPRE_BoomerAMGSetDofFunc) for Hypre and the CPU/GPU +// hybrid solvers: an explicit mapping that reproduces HYPRE's own default +// interleaved dof_func (row i -> function i % block_size) should converge in the +// same number of iterations as never calling set_block_mapping at all. For +// CPUHybrid in particular this also exercises slicing the (globally-indexed) +// mapping down to each MPI rank's local row range. +TEST_CASE("block_mapping", "[.][solver]") +{ + const std::string path = POLYFEM_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_contact.mtx"); + REQUIRE(ok); + + const int block_size = 3; + REQUIRE(A.rows() % block_size == 0); + + Eigen::VectorXi block_mapping(A.rows()); + for (int i = 0; i < A.rows(); ++i) + block_mapping[i] = i % block_size; + + std::vector solvers; +#ifdef POLYSOLVE_WITH_HYPRE + solvers.push_back("Hypre"); +#endif +#ifdef POLYSOLVE_WITH_CPU_HYBRID + solvers.push_back("CPUHybrid"); +#endif +#ifdef POLYSOLVE_WITH_GPU_HYBRID + solvers.push_back("GPUHybrid"); +#endif + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + + for (const auto &s : solvers) + { + json params; + if (s == "Hypre") + params[s]["tolerance"] = 1e-12; + else + { + params[s]["relative_tolerance"] = 0; + params[s]["absolute_tolerance"] = 1e-12; + } + + // Baseline: default interleaved dof_func, set_block_mapping never called. + auto baseline_solver = Solver::create(s, ""); + baseline_solver->set_block_size(block_size); + baseline_solver->set_parameters(params); + + Eigen::VectorXd x_baseline(b.size()); + x_baseline.setZero(); + baseline_solver->analyze_pattern(A, A.rows()); + baseline_solver->factorize(A); + baseline_solver->solve(b, x_baseline); + + json baseline_info; + baseline_solver->get_info(baseline_info); + + // Explicit mapping reproducing that same interleaved assignment. + auto mapped_solver = Solver::create(s, ""); + mapped_solver->set_block_size(block_size); + mapped_solver->set_block_mapping(block_mapping); + mapped_solver->set_parameters(params); + + Eigen::VectorXd x_mapped(b.size()); + x_mapped.setZero(); + mapped_solver->analyze_pattern(A, A.rows()); + mapped_solver->factorize(A); + mapped_solver->solve(b, x_mapped); + + json mapped_info; + mapped_solver->get_info(mapped_info); + + INFO("solver: " + s); + const double err = (A * x_mapped - b).norm(); + REQUIRE(err < 1e-3); + // An explicit dof_func matching HYPRE's own default interleaved mapping + // should make AMG behave identically to never calling set_block_mapping. + REQUIRE(mapped_info["num_iterations"] == baseline_info["num_iterations"]); + } +} + +#ifdef POLYSOLVE_WITH_GPU_HYBRID +// RestartedAMGPCG recomputes b - Ax in double-double precision every +// restart_interval iterations, restarts CG from it, and only declares +// convergence on that recomputed residual. Checks that it converges both with +// frequent restarts and with (effectively) none, that it restarts at least as +// often as restart_interval says, and that the residual it reports matches +// one computed in quad precision on the host. +TEST_CASE("restarted_amg_pcg", "[solver]") +{ + const std::string path = POLYFEM_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_2.mat"); + REQUIRE(ok); + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + + for (const int restart_interval : {5, 50, 1000000}) + { + auto solver = Solver::create("RestartedAMGPCG", ""); + json params; + params["RestartedAMGPCG"]["relative_tolerance"] = 0; + params["RestartedAMGPCG"]["absolute_tolerance"] = 1e-10; + params["RestartedAMGPCG"]["restart_interval"] = restart_interval; + solver->set_parameters(params); + + Eigen::VectorXd x(b.size()); + x.setZero(); + + solver->analyze_pattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + json info; + solver->get_info(info); + const int num_iterations = info["num_iterations"]; + const int num_restarts = info["num_restarts"]; + const double final_res_norm = info["final_res_norm"]; + + INFO("restart_interval: " << restart_interval << ", num_iterations: " << num_iterations << ", num_restarts: " << num_restarts); + REQUIRE(final_res_norm < 1e-10); + REQUIRE(num_restarts >= num_iterations / restart_interval); + +#ifdef __SIZEOF_FLOAT128__ + // A product of two doubles is exact in binary128, so this reference + // residual is accurate to ~1e-34 relative to |A||x|, far beyond both + // the double-double residual and the tolerance here. + const Eigen::SparseMatrix A_rows = A; + __float128 ref_sq_norm = 0; + for (int i = 0; i < A_rows.outerSize(); ++i) + { + __float128 r_i = b[i]; + for (Eigen::SparseMatrix::InnerIterator it(A_rows, i); it; ++it) + r_i -= static_cast<__float128>(it.value()) * static_cast<__float128>(x[it.col()]); + ref_sq_norm += r_i * r_i; + } + const double ref_res_norm = std::sqrt(static_cast(ref_sq_norm)); + + INFO("reported residual: " << final_res_norm << ", quad-precision residual: " << ref_res_norm << ", double-precision residual: " << (A * x - b).norm()); + REQUIRE(std::abs(final_res_norm - ref_res_norm) <= 1e-10 * ref_res_norm); +#endif + } + + { + auto solver = Solver::create("RestartedAMGPCG", ""); + json params; + params["RestartedAMGPCG"]["restart_interval"] = 0; + REQUIRE_THROWS(solver->set_parameters(params)); + } +} +#endif + +#if defined(POLYSOLVE_WITH_CPU_HYBRID) && defined(POLYSOLVE_WITH_STRUMPACK) +// use_strumpack hands every subdomain over 1000 dofs to STRUMPACK, but the +// default selection only picks ~700 dofs of A_contact.mtx. So pick the +// subspace here instead -- the 1200 dofs with the largest L1 row norm, the +// same heuristic with a wider cut -- and hand it over as a single subdomain +// (AMGF mode, no decomposition or expansion). Uncompressed, STRUMPACK's LU and +// the default LDLT factorize the same subdomain matrix, so PCG should need the +// same number of iterations up to rounding; compressed, it only has to converge. +TEST_CASE("cpu_hybrid_strumpack", "[.][solver]") +{ + const std::string path = POLYFEM_DATA_DIR; + Eigen::SparseMatrix A; + const bool ok = loadMarket(A, path + "/A_contact.mtx"); + REQUIRE(ok); + + // A is symmetric, so its column norms are its row norms. + std::vector> row_norms; + for (int k = 0; k < A.outerSize(); ++k) + row_norms.emplace_back(A.col(k).cwiseAbs().sum(), k); + std::sort(row_norms.rbegin(), row_norms.rend()); + + std::set bad_dofs; + for (int i = 0; i < 1200; ++i) + bad_dofs.insert(row_norms[i].second); + + Eigen::VectorXd b(A.rows()); + b.setRandom(); + + const auto solve = [&](const bool use_strumpack, const std::string &compression, json &info) { + auto solver = Solver::create("CPUHybrid", ""); + json params; + params["CPUHybrid"]["relative_tolerance"] = 0; + params["CPUHybrid"]["absolute_tolerance"] = 1e-12; + params["CPUHybrid"]["select_bad_dofs_from_l1_norm"] = false; + params["CPUHybrid"]["decompose_subdomains"] = false; + params["CPUHybrid"]["expand_subdomains"] = false; + params["CPUHybrid"]["use_strumpack"] = use_strumpack; + params["CPUHybrid"]["strumpack_compression"] = compression; + + solver->set_block_size(3); + solver->set_parameters(params); + solver->set_problematic_dofs(bad_dofs); + + Eigen::VectorXd x(b.size()); + x.setZero(); + solver->analyze_pattern(A, A.rows()); + solver->factorize(A); + solver->solve(b, x); + + solver->get_info(info); + return (A * x - b).norm(); + }; + + json baseline_info, strumpack_info, blr_info; + const double baseline_err = solve(false, "none", baseline_info); + const double strumpack_err = solve(true, "none", strumpack_info); + const double blr_err = solve(true, "blr", blr_info); + + REQUIRE(baseline_err < 1e-4); + REQUIRE(strumpack_err < 1e-4); + REQUIRE(blr_err < 1e-4); + + const int baseline_iters = baseline_info["num_iterations"]; + const int strumpack_iters = strumpack_info["num_iterations"]; + REQUIRE(std::abs(strumpack_iters - baseline_iters) <= 2); +} +#endif + TEST_CASE("pre_factor", "[solver]") { const std::string path = POLYFEM_DATA_DIR; @@ -576,7 +884,7 @@ TEST_CASE("amgcl_blocksolver_small_scale", "[solver]") auto solver = Solver::create("AMGCL", ""); json params; params["AMGCL"]["tolerance"] = 1e-8; - params["AMGCL"]["block_size"] = 3; + solver->set_block_size(3); solver->set_parameters(params); solver->analyze_pattern(A, A.rows()); solver->factorize(A); @@ -634,7 +942,7 @@ TEST_CASE("amgcl_blocksolver_b2", "[solver]") json params; params["AMGCL"]["tolerance"] = 1e-8; params["AMGCL"]["max_iter"] = 1000; - params["AMGCL"]["block_size"] = 2; + solver->set_block_size(2); solver->set_parameters(params); solver->analyze_pattern(A, A.rows()); solver->factorize(A); @@ -676,7 +984,7 @@ TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]") json params; params["AMGCL"]["tolerance"] = 1e-8; params["AMGCL"]["max_iter"] = 1000; - params["AMGCL"]["block_size"] = 3; + solver->set_block_size(3); solver->set_parameters(params); solver->analyze_pattern(A, A.rows()); solver->factorize(A); @@ -740,8 +1048,8 @@ TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]") json params; params["AMGCL"]["tolerance"] = 1e-8; params["AMGCL"]["max_iter"] = 10000; - params["AMGCL"]["block_size"] = 3; params["AMGCL"]["solver_type"] = "bicgstab"; + solver->set_block_size(3); solver->set_parameters(params); solver->analyze_pattern(A, A.rows()); solver->factorize(A); diff --git a/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index 805ba7d3..f8913dc7 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -482,6 +482,33 @@ TEST_CASE("nonlinear-fallbacks", "[solver]") } } +#ifdef POLYSOLVE_WITH_CPU_HYBRID +// NewtonCG requires a Hybrid CG-based linear solver (see Newton.hpp). This +// just checks that the new force_psd_projection option is accepted by the +// JSON spec and threaded through NewtonCG's constructor without throwing -- +// not a full solve, since the underlying CPUHybridSolver's MPI-rank +// partitioning isn't meant for a problem this tiny. +TEST_CASE("newton_cg_force_psd_projection", "[.][solver]") +{ + json solver_params = R"({"solver": [ + { + "type": "NewtonCG", + "force_psd_projection": true + } + ] + })"_json; + json linear_solver_params; + linear_solver_params["solver"] = "CPUHybrid"; + + const double characteristic_length = 1; + + static std::shared_ptr logger = spdlog::stdout_color_mt("newton-cg-psd-test-logger"); + logger->set_level(spdlog::level::info); + + REQUIRE_NOTHROW(Solver::create(solver_params, linear_solver_params, characteristic_length, *logger)); +} +#endif + TEST_CASE("nonlinear-gradient-fd", "[solver]") { test_solvers_gradient_fd(false);