Skip to content

Block structure #119

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
maxpaik16 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Block structure #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2f8eff4
Refactor iterative solvers to set block size via set_block_size
maxpaik16 Sep 9, 2026
32e0baf
Implement set_block_mapping for Hypre and CPU/GPU hybrid solvers
maxpaik16 Sep 9, 2026
ecb86ee
Let Problem supply a block mapping for the linear solver
maxpaik16 Sep 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 1 addition & 34 deletions linear-solver-spec.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,8 +148,7 @@
"tolerance",
"theta",
"nodal_coarsening",
"interp_rbms",
"dimension"
"interp_rbms"
],
"doc": "Settings for the Hypre solver."
},
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -461,7 +454,6 @@
"default": null,
"type": "object",
"optional": [
"block_dim",
"max_iter",
"relative_tolerance",
"absolute_tolerance",
Expand All @@ -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,
Expand Down Expand Up @@ -516,7 +497,6 @@
"default": null,
"type": "object",
"optional": [
"block_dim",
"max_iter",
"relative_tolerance",
"absolute_tolerance",
Expand All @@ -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,
Expand Down Expand Up @@ -622,7 +596,6 @@
"default": null,
"type": "object",
"optional": [
"block_dim",
"max_iter",
"relative_tolerance",
"absolute_tolerance",
Expand All @@ -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,
Expand Down
10 changes: 5 additions & 5 deletions src/polysolve/linear/AMGCL.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -127,6 +122,11 @@ namespace polysolve::linear
}
}

void AMGCL::set_block_size(int block_size)
{
block_size_ = block_size;
}

void AMGCL::get_info(json &params) const
{
if (block_size_ == 2)
Expand Down
3 changes: 3 additions & 0 deletions src/polysolve/linear/AMGCL.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,9 @@ namespace polysolve::linear
// Set solver parameters
virtual void set_parameters(const json &params) override;

// Set block size for multigrid solvers
virtual void set_block_size(int block_size) override;

// Retrieve information
virtual void get_info(json &params) const override;

Expand Down
72 changes: 68 additions & 4 deletions src/polysolve/linear/CPUHybridSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"];
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down
16 changes: 16 additions & 0 deletions src/polysolve/linear/CPUHybridSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -104,6 +106,12 @@ namespace polysolve::linear
// Set solver parameters
virtual void set_parameters(const json &params) 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 &params) const override;

Expand Down Expand Up @@ -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;
Expand Down
39 changes: 33 additions & 6 deletions src/polysolve/linear/GPUHybridSolver.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,11 @@

#include <spdlog/spdlog.h>

// 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 <mpi.h>
#endif
Expand Down Expand Up @@ -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"];
Expand Down Expand Up @@ -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 &params) const
{
params["num_iterations"] = num_iterations;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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<HYPRE_Int> 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);
Expand Down
10 changes: 10 additions & 0 deletions src/polysolve/linear/GPUHybridSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,12 @@ namespace polysolve::linear
// Set solver parameters
virtual void set_parameters(const json &params) 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 &params) const override;

Expand Down Expand Up @@ -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;
Expand Down
Loading
Loading