diff --git a/linear-solver-spec.json b/linear-solver-spec.json index fbf698c4..700ceca0 100644 --- a/linear-solver-spec.json +++ b/linear-solver-spec.json @@ -148,8 +148,7 @@ "tolerance", "theta", "nodal_coarsening", - "interp_rbms", - "dimension" + "interp_rbms" ], "doc": "Settings for the Hypre solver." }, @@ -288,12 +287,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 +454,6 @@ "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -470,17 +462,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 +497,6 @@ "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -533,12 +513,6 @@ ], "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, @@ -622,7 +596,6 @@ "default": null, "type": "object", "optional": [ - "block_dim", "max_iter", "relative_tolerance", "absolute_tolerance", @@ -639,12 +612,6 @@ ], "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, 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/CPUHybridSolver.cpp b/src/polysolve/linear/CPUHybridSolver.cpp index 39cf9d6e..361e8daf 100644 --- a/src/polysolve/linear/CPUHybridSolver.cpp +++ b/src/polysolve/linear/CPUHybridSolver.cpp @@ -131,10 +131,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"]; @@ -174,6 +170,47 @@ namespace polysolve::linear } } + // 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 { if (myid != 0) @@ -483,6 +520,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); @@ -1450,6 +1502,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..c346a291 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; @@ -137,6 +145,14 @@ namespace polysolve::linear // 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; diff --git a/src/polysolve/linear/GPUHybridSolver.cu b/src/polysolve/linear/GPUHybridSolver.cu index e34ded99..f7b92ebb 100644 --- a/src/polysolve/linear/GPUHybridSolver.cu +++ b/src/polysolve/linear/GPUHybridSolver.cu @@ -30,6 +30,11 @@ #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 @@ -116,10 +121,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"]; @@ -159,6 +160,16 @@ namespace polysolve::linear } } + 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; @@ -296,8 +307,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,6 +341,24 @@ 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); diff --git a/src/polysolve/linear/GPUHybridSolver.hpp b/src/polysolve/linear/GPUHybridSolver.hpp index b71d656a..97504e4f 100644 --- a/src/polysolve/linear/GPUHybridSolver.hpp +++ b/src/polysolve/linear/GPUHybridSolver.hpp @@ -43,6 +43,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; @@ -76,6 +82,10 @@ namespace polysolve::linear // 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; 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..ddb08bc9 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_; @@ -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/Solver.cpp b/src/polysolve/linear/Solver.cpp index 9bb5070c..8cc4cce3 100644 --- a/src/polysolve/linear/Solver.cpp +++ b/src/polysolve/linear/Solver.cpp @@ -146,7 +146,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 +165,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; diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index c6837056..a5da674f 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -56,11 +56,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,6 +112,9 @@ 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) {} diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 92711ef7..cad6a072 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -34,6 +34,17 @@ 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 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..1c05c704 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -37,47 +37,48 @@ 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 == "LBFGS" || solver_name == "L-BFGS") @@ -127,7 +128,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 +151,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 +160,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") 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/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index ff8d015a..f5d45685 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.cpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.cpp @@ -17,7 +17,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 +36,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 +64,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 +84,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 +96,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 +109,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); @@ -185,6 +190,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); + // TODO: get the correct size linear_solver->analyze_pattern(hessian, hessian.rows()); @@ -228,6 +237,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..252cfe77 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"; } @@ -89,7 +92,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 +117,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 { diff --git a/tests/test_linear_solver.cpp b/tests/test_linear_solver.cpp index 318b6e89..a36eeac4 100644 --- a/tests/test_linear_solver.cpp +++ b/tests/test_linear_solver.cpp @@ -223,10 +223,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 +264,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 +289,89 @@ TEST_CASE("hybrid_convergence", "[.][solver]") } } +// 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"]); + } +} + TEST_CASE("pre_factor", "[solver]") { const std::string path = POLYFEM_DATA_DIR; @@ -576,7 +659,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 +717,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 +759,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 +823,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);