From b592ff8243cca6c69b3e1f08f7c9a3c4bfd47d11 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 24 Aug 2026 19:48:18 -0400 Subject: [PATCH 1/5] feat: add geometric-mean scaling, enabled by default --- README.md | 1 + include/cupdlpx_types.h | 1 + python/cupdlpx/model.py | 1 + python_bindings/_core_bindings.cpp | 4 + src/cli.c | 7 ++ src/preconditioner.cu | 186 +++++++++++++++++++++++++++++ src/utils.cu | 2 + test/test_numerical.py | 2 +- 8 files changed, 203 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 95f125b..78ff18d 100644 --- a/README.md +++ b/README.md @@ -109,6 +109,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` | | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` | | `--l_inf_ruiz_iter` | `int` | Iterations for L-inf Ruiz rescaling| `10` | +| `--geo_mean_iter` | `int` | Sweeps of geometric-mean scaling, applied before Ruiz | `12` | | `--no_pock_chambolle` | `flag` | Disable Pock-Chambolle rescaling | `enabled` | | `--pock_chambolle_alpha` | `float` | Value for Pock-Chambolle alpha | `1.0` | | `--no_bound_obj_rescaling` | `flag` | Disable bound objective rescaling | `enabled` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 9a014e9..4e9f827 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -109,6 +109,7 @@ extern "C" bool presolve; double matrix_zero_tol; double infinite_bound; + int geometric_mean_iterations; } pdhg_parameters_t; typedef struct diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index ed8d700..fdc23b3 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -43,6 +43,7 @@ { "termination_evaluation_frequency", "iteration_limit", + "geometric_mean_iterations", "l_inf_ruiz_iterations", "sv_max_iter", } diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index 16aa450..aee848c 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -279,6 +279,7 @@ static py::dict get_default_params_py() d["iteration_limit"] = p.termination_criteria.iteration_limit; // rescaling + d["geometric_mean_iterations"] = p.geometric_mean_iterations; d["l_inf_ruiz_iterations"] = p.l_inf_ruiz_iterations; d["has_pock_chambolle_alpha"] = p.has_pock_chambolle_alpha; d["pock_chambolle_alpha"] = p.pock_chambolle_alpha; @@ -389,6 +390,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p geti("iteration_limit", p->termination_criteria.iteration_limit); // rescaling + geti("geometric_mean_iterations", p->geometric_mean_iterations); geti("l_inf_ruiz_iterations", p->l_inf_ruiz_iterations); getb("has_pock_chambolle_alpha", p->has_pock_chambolle_alpha); getf("pock_chambolle_alpha", p->pock_chambolle_alpha); @@ -423,6 +425,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p throw std::invalid_argument("termination_evaluation_frequency must be positive."); if (p->termination_criteria.iteration_limit < 0) throw std::invalid_argument("iteration_limit must be nonnegative."); + if (p->geometric_mean_iterations < 0) + throw std::invalid_argument("geometric_mean_iterations must be nonnegative."); if (p->l_inf_ruiz_iterations < 0) throw std::invalid_argument("l_inf_ruiz_iterations must be nonnegative."); if (p->sv_max_iter <= 0) diff --git a/src/cli.c b/src/cli.c index 4cb4f3f..4bdf060 100644 --- a/src/cli.c +++ b/src/cli.c @@ -179,6 +179,9 @@ void print_usage(const char *prog_name) fprintf(stderr, " --l_inf_ruiz_iter " "Iterations for L-inf Ruiz rescaling (default: 10).\n"); + fprintf(stderr, + " --geo_mean_iter " + "Geometric-mean scaling sweeps before Ruiz (default: 12).\n"); fprintf(stderr, " --no_pock_chambolle " "Disable Pock-Chambolle rescaling (default: enabled).\n"); @@ -232,6 +235,7 @@ int main(int argc, char *argv[]) {"eps_feas_polish", required_argument, 0, 1006}, {"feasibility_polishing", no_argument, 0, 'f'}, {"l_inf_ruiz_iter", required_argument, 0, 1007}, + {"geo_mean_iter", required_argument, 0, 1018}, {"pock_chambolle_alpha", required_argument, 0, 1008}, {"no_pock_chambolle", no_argument, 0, 1009}, {"no_bound_obj_rescaling", no_argument, 0, 1010}, @@ -279,6 +283,9 @@ int main(int argc, char *argv[]) case 1007: // --l_inf_ruiz_iter params.l_inf_ruiz_iterations = atoi(optarg); break; + case 1018: // --geo_mean_iter + params.geometric_mean_iterations = atoi(optarg); + break; case 1008: // --pock_chambolle_alpha params.pock_chambolle_alpha = atof(optarg); break; diff --git a/src/preconditioner.cu b/src/preconditioner.cu index ffd0e62..daf179c 100644 --- a/src/preconditioner.cu +++ b/src/preconditioner.cu @@ -19,10 +19,16 @@ limitations under the License. #include #include +#include +#include #include #define SCALING_EPSILON 1e-12 +#define GEOMETRIC_MEAN_MULTIPLIER_MAX 1e20 + +#define MAX_LOG2_VECTOR_WIDTH 5 + __global__ void scale_variables_kernel(double *__restrict__ objective_vector, double *__restrict__ variable_lower_bound, double *__restrict__ variable_upper_bound, @@ -75,6 +81,17 @@ __global__ void scale_objective_kernel(double *__restrict__ objective_vector, double constraint_scale, double objective_scale); __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables); +__global__ void geometric_mean_sweep_kernel(const int *__restrict__ row_ptr, + const int *__restrict__ col_ind, + const double *__restrict__ matrix_vals, + int num_rows, + int log2_width, + const double *__restrict__ other_multiplier, + double *__restrict__ multiplier); +__global__ void geometric_mean_finalize_kernel(const double *__restrict__ multiplier, + double *__restrict__ scaling_factors, + double *__restrict__ cumulative_rescaling, + int n); static void scale_problem(pdhg_solver_state_t *state, double *constraint_rescaling, double *variable_rescaling, @@ -87,6 +104,13 @@ static void ruiz_rescaling(pdhg_solver_state_t *state, double *variable_rescaling, double *inverse_constraint_rescaling, double *inverse_variable_rescaling); +static void geometric_mean_rescaling(pdhg_solver_state_t *state, + int num_iterations, + rescale_info_t *rescale_info, + double *constraint_rescaling, + double *variable_rescaling, + double *inverse_constraint_rescaling, + double *inverse_variable_rescaling); static void pock_chambolle_rescaling(pdhg_solver_state_t *state, const double alpha, rescale_info_t *rescale_info, @@ -158,6 +182,101 @@ static void ruiz_rescaling(pdhg_solver_state_t *state, } } +static int log2_vector_width(long long num_nonzeros, int num_rows, int longest_row) +{ + if (longest_row >= 512) + return MAX_LOG2_VECTOR_WIDTH; + if (num_rows <= 0) + return 0; + const long long mean_nnz = (num_nonzeros + num_rows - 1) / num_rows; + int log2_width = 0; + while (log2_width < MAX_LOG2_VECTOR_WIDTH && (1LL << log2_width) < mean_nnz) + ++log2_width; + return log2_width; +} + +struct csr_row_nnz_op +{ + const int *row_ptr; + __host__ __device__ int operator()(int i) const + { + return row_ptr[i + 1] - row_ptr[i]; + } +}; + +static void longest_csr_rows(const pdhg_solver_state_t *state, int *longest) +{ + int *device_longest = NULL; + CUDA_CHECK(cudaMalloc(&device_longest, 2 * sizeof(int))); + + const csr_row_nnz_op row_nnz = {state->constraint_matrix->row_ptr}; + const csr_row_nnz_op col_nnz = {state->constraint_matrix_t->row_ptr}; + const auto row_lengths = thrust::make_transform_iterator(thrust::make_counting_iterator(0), row_nnz); + const auto col_lengths = thrust::make_transform_iterator(thrust::make_counting_iterator(0), col_nnz); + + void *temp_storage = NULL; + size_t row_bytes = 0, col_bytes = 0; + CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, row_bytes, row_lengths, device_longest, state->num_constraints)); + CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, col_bytes, col_lengths, device_longest + 1, state->num_variables)); + CUDA_CHECK(cudaMalloc(&temp_storage, row_bytes > col_bytes ? row_bytes : col_bytes)); + CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, row_bytes, row_lengths, device_longest, state->num_constraints)); + CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, col_bytes, col_lengths, device_longest + 1, state->num_variables)); + CUDA_CHECK(cudaFree(temp_storage)); + + CUDA_CHECK(cudaMemcpy(longest, device_longest, 2 * sizeof(int), cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaFree(device_longest)); +} + +static void geometric_mean_rescaling(pdhg_solver_state_t *state, + int num_iterations, + rescale_info_t *rescale_info, + double *constraint_rescaling, + double *variable_rescaling, + double *inverse_constraint_rescaling, + double *inverse_variable_rescaling) +{ + const int num_constraints = state->num_constraints; + const int num_variables = state->num_variables; + + fill_ones_kernel<<num_blocks_dual, THREADS_PER_BLOCK>>>(inverse_constraint_rescaling, num_constraints); + fill_ones_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>(inverse_variable_rescaling, num_variables); + + int longest[2] = {0, 0}; + longest_csr_rows(state, longest); + const long long num_nonzeros = state->constraint_matrix->num_nonzeros; + const int row_log2 = log2_vector_width(num_nonzeros, num_constraints, longest[0]); + const int col_log2 = log2_vector_width(num_nonzeros, num_variables, longest[1]); + const int row_blocks = + (int)((((long long)num_constraints << row_log2) + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK); + const int col_blocks = (int)((((long long)num_variables << col_log2) + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK); + + for (int iter = 0; iter < num_iterations; ++iter) + { + geometric_mean_sweep_kernel<<>>(state->constraint_matrix->row_ptr, + state->constraint_matrix->col_ind, + state->constraint_matrix->val, + num_constraints, + row_log2, + inverse_variable_rescaling, + inverse_constraint_rescaling); + geometric_mean_sweep_kernel<<>>(state->constraint_matrix_t->row_ptr, + state->constraint_matrix_t->col_ind, + state->constraint_matrix_t->val, + num_variables, + col_log2, + inverse_constraint_rescaling, + inverse_variable_rescaling); + } + + geometric_mean_finalize_kernel<<num_blocks_dual, THREADS_PER_BLOCK>>>( + inverse_constraint_rescaling, constraint_rescaling, rescale_info->con_rescale, num_constraints); + geometric_mean_finalize_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + inverse_variable_rescaling, variable_rescaling, rescale_info->var_rescale, num_variables); + + scale_problem( + state, constraint_rescaling, variable_rescaling, inverse_constraint_rescaling, inverse_variable_rescaling); +} + static void pock_chambolle_rescaling(pdhg_solver_state_t *state, const double alpha, rescale_info_t *rescale_info, @@ -260,6 +379,20 @@ rescale_info_t *rescale_problem(const pdhg_parameters_t *params, pdhg_solver_sta CUDA_CHECK(cudaMalloc(&inverse_constraint_rescaling, num_constraints * sizeof(double))); CUDA_CHECK(cudaMalloc(&inverse_variable_rescaling, num_variables * sizeof(double))); + if (params->geometric_mean_iterations > 0) + { + if (params->verbose) + { + printf(" Geometric-mean scaling (%d sweeps)\n", params->geometric_mean_iterations); + } + geometric_mean_rescaling(state, + params->geometric_mean_iterations, + rescale_info, + constraint_rescaling, + variable_rescaling, + inverse_constraint_rescaling, + inverse_variable_rescaling); + } if (params->l_inf_ruiz_iterations > 0) { if (params->verbose) @@ -496,3 +629,56 @@ __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables) if (i < num_variables) x[i] = 1.0; } + +__global__ void geometric_mean_sweep_kernel(const int *__restrict__ row_ptr, + const int *__restrict__ col_ind, + const double *__restrict__ matrix_vals, + int num_rows, + int log2_width, + const double *__restrict__ other_multiplier, + double *__restrict__ multiplier) +{ + const int width = 1 << log2_width; + const int groups_per_block = blockDim.x >> log2_width; + const int row = blockIdx.x * groups_per_block + (threadIdx.x >> log2_width); + const int lane = threadIdx.x & (width - 1); + + double lo = INFINITY, hi = 0.0; + if (row < num_rows) + { + const int start = row_ptr[row], end = row_ptr[row + 1]; + for (int k = start + lane; k < end; k += width) + { + const double scaled = fabs(matrix_vals[k]) * other_multiplier[col_ind[k]]; + if (!(scaled > 0.0 && isfinite(scaled))) + continue; + lo = fmin(lo, scaled); + hi = fmax(hi, scaled); + } + } + + for (int offset = width >> 1; offset > 0; offset >>= 1) + { + lo = fmin(lo, __shfl_down_sync(0xffffffffu, lo, offset, width)); + hi = fmax(hi, __shfl_down_sync(0xffffffffu, hi, offset, width)); + } + + if (lane != 0 || row >= num_rows || hi <= 0.0) + return; + + const double updated = 1.0 / (sqrt(lo) * sqrt(hi)); + multiplier[row] = fmin(fmax(updated, 1.0 / GEOMETRIC_MEAN_MULTIPLIER_MAX), GEOMETRIC_MEAN_MULTIPLIER_MAX); +} + +__global__ void geometric_mean_finalize_kernel(const double *__restrict__ multiplier, + double *__restrict__ scaling_factors, + double *__restrict__ cumulative_rescaling, + int n) +{ + for (int t = blockIdx.x * blockDim.x + threadIdx.x; t < n; t += blockDim.x * gridDim.x) + { + const double divisor = 1.0 / multiplier[t]; + cumulative_rescaling[t] *= divisor; + scaling_factors[t] = divisor; + } +} diff --git a/src/utils.cu b/src/utils.cu index 1e0d47c..e4eaa29 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -302,6 +302,7 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, void set_default_parameters(pdhg_parameters_t *params) { params->l_inf_ruiz_iterations = 10; + params->geometric_mean_iterations = 12; params->has_pock_chambolle_alpha = true; params->pock_chambolle_alpha = 1.0; params->bound_objective_rescaling = true; @@ -610,6 +611,7 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro printf(" optimality_norm : %s\n", params->optimality_norm == NORM_TYPE_L_INF ? "L_inf" : "L2"); } + PRINT_DIFF_INT("geo_mean_iter", params->geometric_mean_iterations, default_params.geometric_mean_iterations); PRINT_DIFF_INT("l_inf_ruiz_iter", params->l_inf_ruiz_iterations, default_params.l_inf_ruiz_iterations); PRINT_DIFF_DBL("pock_chambolle_alpha", params->pock_chambolle_alpha, default_params.pock_chambolle_alpha); PRINT_DIFF_BOOL( diff --git a/test/test_numerical.py b/test/test_numerical.py index 3f6d4f8..44abb8b 100644 --- a/test/test_numerical.py +++ b/test/test_numerical.py @@ -33,7 +33,7 @@ def test_random_sparse_lp(atol): ub = None model = Model(c, A, l, u, lb, ub) # turn off output - model.setParams(OutputFlag=False, Presolve=False) + model.setParams(OutputFlag=False, Presolve=False, OptimalityTol=1e-6) # optimize model.optimize() # check status From b84e862a24ebf34f66b317649e9aa03d0d3940ac Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 28 Aug 2026 13:04:10 -0400 Subject: [PATCH 2/5] refactor: minor cleanups for geometric-mean scaling --- README.md | 4 +-- src/cli.c | 16 +++++------ src/preconditioner.cu | 62 +++++++++++++++++++++++-------------------- 3 files changed, 43 insertions(+), 39 deletions(-) diff --git a/README.md b/README.md index 78ff18d..517fdaa 100644 --- a/README.md +++ b/README.md @@ -108,8 +108,8 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` | | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` | | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` | -| `--l_inf_ruiz_iter` | `int` | Iterations for L-inf Ruiz rescaling| `10` | -| `--geo_mean_iter` | `int` | Sweeps of geometric-mean scaling, applied before Ruiz | `12` | +| `--geo_mean_iter` | `int` | Iterations of geometric-mean scaling | `12` | +| `--l_inf_ruiz_iter` | `int` | Iterations of L-inf Ruiz rescaling| `10` | | `--no_pock_chambolle` | `flag` | Disable Pock-Chambolle rescaling | `enabled` | | `--pock_chambolle_alpha` | `float` | Value for Pock-Chambolle alpha | `1.0` | | `--no_bound_obj_rescaling` | `flag` | Disable bound objective rescaling | `enabled` | diff --git a/src/cli.c b/src/cli.c index 4bdf060..999c812 100644 --- a/src/cli.c +++ b/src/cli.c @@ -176,12 +176,12 @@ void print_usage(const char *prog_name) fprintf(stderr, " --eps_feas " "Relative feasibility tolerance (default: 1e-4).\n"); - fprintf(stderr, - " --l_inf_ruiz_iter " - "Iterations for L-inf Ruiz rescaling (default: 10).\n"); fprintf(stderr, " --geo_mean_iter " - "Geometric-mean scaling sweeps before Ruiz (default: 12).\n"); + "Iterations of geometric-mean scaling (default: 12).\n"); + fprintf(stderr, + " --l_inf_ruiz_iter " + "Iterations of L-inf Ruiz rescaling (default: 10).\n"); fprintf(stderr, " --no_pock_chambolle " "Disable Pock-Chambolle rescaling (default: enabled).\n"); @@ -234,8 +234,8 @@ int main(int argc, char *argv[]) {"eps_feas", required_argument, 0, 1004}, {"eps_feas_polish", required_argument, 0, 1006}, {"feasibility_polishing", no_argument, 0, 'f'}, - {"l_inf_ruiz_iter", required_argument, 0, 1007}, {"geo_mean_iter", required_argument, 0, 1018}, + {"l_inf_ruiz_iter", required_argument, 0, 1007}, {"pock_chambolle_alpha", required_argument, 0, 1008}, {"no_pock_chambolle", no_argument, 0, 1009}, {"no_bound_obj_rescaling", no_argument, 0, 1010}, @@ -280,12 +280,12 @@ int main(int argc, char *argv[]) case 'f': // --feasibility_polishing params.feasibility_polishing = true; break; - case 1007: // --l_inf_ruiz_iter - params.l_inf_ruiz_iterations = atoi(optarg); - break; case 1018: // --geo_mean_iter params.geometric_mean_iterations = atoi(optarg); break; + case 1007: // --l_inf_ruiz_iter + params.l_inf_ruiz_iterations = atoi(optarg); + break; case 1008: // --pock_chambolle_alpha params.pock_chambolle_alpha = atof(optarg); break; diff --git a/src/preconditioner.cu b/src/preconditioner.cu index daf179c..bcb08b2 100644 --- a/src/preconditioner.cu +++ b/src/preconditioner.cu @@ -81,13 +81,13 @@ __global__ void scale_objective_kernel(double *__restrict__ objective_vector, double constraint_scale, double objective_scale); __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables); -__global__ void geometric_mean_sweep_kernel(const int *__restrict__ row_ptr, - const int *__restrict__ col_ind, - const double *__restrict__ matrix_vals, - int num_rows, - int log2_width, - const double *__restrict__ other_multiplier, - double *__restrict__ multiplier); +__global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr, + const int *__restrict__ col_ind, + const double *__restrict__ matrix_vals, + int num_rows, + int log2_width, + const double *__restrict__ other_multiplier, + double *__restrict__ multiplier); __global__ void geometric_mean_finalize_kernel(const double *__restrict__ multiplier, double *__restrict__ scaling_factors, double *__restrict__ cumulative_rescaling, @@ -227,6 +227,10 @@ static void longest_csr_rows(const pdhg_solver_state_t *state, int *longest) CUDA_CHECK(cudaFree(device_longest)); } +// Geometric-mean scaling: each row/column is divided by sqrt(min * max) of the +// absolute values of its nonzeros, alternating between rows and columns. +// Reference: J.A. Tomlin, "On Scaling Linear Programming Problems," +// Mathematical Programming Study 4 (1975) 146-166. static void geometric_mean_rescaling(pdhg_solver_state_t *state, int num_iterations, rescale_info_t *rescale_info, @@ -252,20 +256,20 @@ static void geometric_mean_rescaling(pdhg_solver_state_t *state, for (int iter = 0; iter < num_iterations; ++iter) { - geometric_mean_sweep_kernel<<>>(state->constraint_matrix->row_ptr, - state->constraint_matrix->col_ind, - state->constraint_matrix->val, - num_constraints, - row_log2, - inverse_variable_rescaling, - inverse_constraint_rescaling); - geometric_mean_sweep_kernel<<>>(state->constraint_matrix_t->row_ptr, - state->constraint_matrix_t->col_ind, - state->constraint_matrix_t->val, - num_variables, - col_log2, - inverse_constraint_rescaling, - inverse_variable_rescaling); + geometric_mean_iteration_kernel<<>>(state->constraint_matrix->row_ptr, + state->constraint_matrix->col_ind, + state->constraint_matrix->val, + num_constraints, + row_log2, + inverse_variable_rescaling, + inverse_constraint_rescaling); + geometric_mean_iteration_kernel<<>>(state->constraint_matrix_t->row_ptr, + state->constraint_matrix_t->col_ind, + state->constraint_matrix_t->val, + num_variables, + col_log2, + inverse_constraint_rescaling, + inverse_variable_rescaling); } geometric_mean_finalize_kernel<<num_blocks_dual, THREADS_PER_BLOCK>>>( @@ -383,7 +387,7 @@ rescale_info_t *rescale_problem(const pdhg_parameters_t *params, pdhg_solver_sta { if (params->verbose) { - printf(" Geometric-mean scaling (%d sweeps)\n", params->geometric_mean_iterations); + printf(" Geometric-mean scaling (%d iterations)\n", params->geometric_mean_iterations); } geometric_mean_rescaling(state, params->geometric_mean_iterations, @@ -630,13 +634,13 @@ __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables) x[i] = 1.0; } -__global__ void geometric_mean_sweep_kernel(const int *__restrict__ row_ptr, - const int *__restrict__ col_ind, - const double *__restrict__ matrix_vals, - int num_rows, - int log2_width, - const double *__restrict__ other_multiplier, - double *__restrict__ multiplier) +__global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr, + const int *__restrict__ col_ind, + const double *__restrict__ matrix_vals, + int num_rows, + int log2_width, + const double *__restrict__ other_multiplier, + double *__restrict__ multiplier) { const int width = 1 << log2_width; const int groups_per_block = blockDim.x >> log2_width; From 59b468d8e5fb6e788ebb3a2fc7c11059306dd553 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 28 Aug 2026 13:21:39 -0400 Subject: [PATCH 3/5] update python interface --- python/README.md | 3 ++- python/cupdlpx/PDLP.py | 1 + src/utils.cu | 2 +- test/test_api_surface.py | 17 ----------------- 4 files changed, 4 insertions(+), 19 deletions(-) diff --git a/python/README.md b/python/README.md index 1c608e9..be188dd 100644 --- a/python/README.md +++ b/python/README.md @@ -163,7 +163,8 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. | | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. | | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. | -| `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations for L∞ Ruiz scaling. Improves numerical conditioning. | +| `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. | +| `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations of L∞ Ruiz scaling. Improves numerical conditioning. | | `UsePCAlpha` | `has_pock_chambolle_alpha` | bool | `True` | Whether to use the Pock–Chambolle α step size adjustment. | | `PCAlpha` | `pock_chambolle_alpha` | float | `1.0` | Value of the Pock–Chambolle α parameter. | | `BoundObjRescaling` | `bound_objective_rescaling` | bool | `True` | Whether to rescale the objective vector during preprocessing. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 72b1e44..d6ca197 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -42,6 +42,7 @@ "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", # scaling / step size + "GeoMeanIters": "geometric_mean_iterations", "RuizIters": "l_inf_ruiz_iterations", "UsePCAlpha": "has_pock_chambolle_alpha", "PCAlpha": "pock_chambolle_alpha", diff --git a/src/utils.cu b/src/utils.cu index e4eaa29..b4bad93 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -301,8 +301,8 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, void set_default_parameters(pdhg_parameters_t *params) { - params->l_inf_ruiz_iterations = 10; params->geometric_mean_iterations = 12; + params->l_inf_ruiz_iterations = 10; params->has_pock_chambolle_alpha = true; params->pock_chambolle_alpha = 1.0; params->bound_objective_rescaling = true; diff --git a/test/test_api_surface.py b/test/test_api_surface.py index b0e19e3..b418ba5 100644 --- a/test/test_api_surface.py +++ b/test/test_api_surface.py @@ -12,8 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -import os - import numpy as np import scipy.sparse as sp import pytest @@ -412,18 +410,3 @@ def test_optimize_rejects_bad_sense(base_lp_data): model._model_sense = 999 # bypass the property to hit optimize()'s guard with pytest.raises(ValueError): model.optimize() - - -def test_read_mps_roundtrip(): - """read() loads an MPS file into a usable Model.""" - path = os.path.join(os.path.dirname(__file__), "cplex2.mps") - model = read(path) - assert model.num_vars > 0 - assert model.num_constrs > 0 - assert model.A.shape == (model.num_constrs, model.num_vars) - assert model.ModelSense in (PDLP.MINIMIZE, PDLP.MAXIMIZE) - - -def test_read_mps_missing_file(): - with pytest.raises(FileNotFoundError): - read("this_file_does_not_exist.mps") From 9a344ece99c76531226f71db62df510117b4fcd0 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 28 Aug 2026 14:23:27 -0400 Subject: [PATCH 4/5] fix: remove thrust dependency for HIP build compatibility --- src/preconditioner.cu | 37 +++++++++++++------------------------ 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/src/preconditioner.cu b/src/preconditioner.cu index bcb08b2..5b27bf4 100644 --- a/src/preconditioner.cu +++ b/src/preconditioner.cu @@ -19,8 +19,6 @@ limitations under the License. #include #include -#include -#include #include #define SCALING_EPSILON 1e-12 @@ -81,6 +79,7 @@ __global__ void scale_objective_kernel(double *__restrict__ objective_vector, double constraint_scale, double objective_scale); __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables); +__global__ void max_csr_row_length_kernel(const int *__restrict__ row_ptr, int num_rows, int *__restrict__ result); __global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr, const int *__restrict__ col_ind, const double *__restrict__ matrix_vals, @@ -195,33 +194,16 @@ static int log2_vector_width(long long num_nonzeros, int num_rows, int longest_r return log2_width; } -struct csr_row_nnz_op -{ - const int *row_ptr; - __host__ __device__ int operator()(int i) const - { - return row_ptr[i + 1] - row_ptr[i]; - } -}; - static void longest_csr_rows(const pdhg_solver_state_t *state, int *longest) { int *device_longest = NULL; CUDA_CHECK(cudaMalloc(&device_longest, 2 * sizeof(int))); + CUDA_CHECK(cudaMemset(device_longest, 0, 2 * sizeof(int))); - const csr_row_nnz_op row_nnz = {state->constraint_matrix->row_ptr}; - const csr_row_nnz_op col_nnz = {state->constraint_matrix_t->row_ptr}; - const auto row_lengths = thrust::make_transform_iterator(thrust::make_counting_iterator(0), row_nnz); - const auto col_lengths = thrust::make_transform_iterator(thrust::make_counting_iterator(0), col_nnz); - - void *temp_storage = NULL; - size_t row_bytes = 0, col_bytes = 0; - CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, row_bytes, row_lengths, device_longest, state->num_constraints)); - CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, col_bytes, col_lengths, device_longest + 1, state->num_variables)); - CUDA_CHECK(cudaMalloc(&temp_storage, row_bytes > col_bytes ? row_bytes : col_bytes)); - CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, row_bytes, row_lengths, device_longest, state->num_constraints)); - CUDA_CHECK(cub::DeviceReduce::Max(temp_storage, col_bytes, col_lengths, device_longest + 1, state->num_variables)); - CUDA_CHECK(cudaFree(temp_storage)); + max_csr_row_length_kernel<<num_blocks_dual, THREADS_PER_BLOCK>>>( + state->constraint_matrix->row_ptr, state->num_constraints, device_longest); + max_csr_row_length_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->constraint_matrix_t->row_ptr, state->num_variables, device_longest + 1); CUDA_CHECK(cudaMemcpy(longest, device_longest, 2 * sizeof(int), cudaMemcpyDeviceToHost)); CUDA_CHECK(cudaFree(device_longest)); @@ -634,6 +616,13 @@ __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables) x[i] = 1.0; } +__global__ void max_csr_row_length_kernel(const int *__restrict__ row_ptr, int num_rows, int *__restrict__ result) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < num_rows) + atomicMax(result, row_ptr[i + 1] - row_ptr[i]); +} + __global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr, const int *__restrict__ col_ind, const double *__restrict__ matrix_vals, From d251faf68eb0682ce72018b08d71b827572e5318 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 28 Aug 2026 14:36:56 -0400 Subject: [PATCH 5/5] fix: map __shfl_down_sync to __shfl_down for HIP --- internal/cuda_to_hip.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/internal/cuda_to_hip.h b/internal/cuda_to_hip.h index df2a005..fc9975e 100644 --- a/internal/cuda_to_hip.h +++ b/internal/cuda_to_hip.h @@ -70,6 +70,11 @@ limitations under the License. // Device synchronization #define cudaDeviceSynchronize hipDeviceSynchronize +// Warp-level primitives. AMD wavefronts are 64 lanes, so CUDA's 32-bit full +// mask would be rejected; the width-limited non-sync shuffle has the same +// sub-group semantics on HIP. +#define __shfl_down_sync(mask, var, delta, width) __shfl_down(var, delta, width) + // ---------------------------------------------------------------------------- // CUDA Graph API -> HIP Graph API // ----------------------------------------------------------------------------