From 545df91779ba24c9b735cd8fe9acc76748c79d3f Mon Sep 17 00:00:00 2001 From: Lhongpei <1453244320@qq.com> Date: Thu, 30 Apr 2026 11:47:05 +0000 Subject: [PATCH 1/4] add space --- distributed/distributed_utils.cu | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/distributed/distributed_utils.cu b/distributed/distributed_utils.cu index 74efc55..adc088b 100644 --- a/distributed/distributed_utils.cu +++ b/distributed/distributed_utils.cu @@ -979,11 +979,12 @@ void print_distributed_params(const pdhg_parameters_t *params) if (params->grid_size.decided) { - printf(" Grid Size : %d x %d (Rows x Cols)\n", params->grid_size.row_dims, params->grid_size.col_dims); + printf( + " Grid Size : %d x %d (Rows x Cols)\n", params->grid_size.row_dims, params->grid_size.col_dims); } else { - printf(" Grid Size : Auto-detect (implementation dependent)\n"); + printf(" Grid Size : Auto-detect (implementation dependent)\n"); } printf(" Partition Method : "); From a92a1b5e1b38b02f08a45df48b49256751bd1f50 Mon Sep 17 00:00:00 2001 From: Lhongpei <1453244320@qq.com> Date: Fri, 1 May 2026 01:22:49 +0000 Subject: [PATCH 2/4] update doc --- docs/index.md | 2 +- docs/installation.md | 2 +- python/README.md | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/index.md b/docs/index.md index a102e43..48de299 100644 --- a/docs/index.md +++ b/docs/index.md @@ -28,7 +28,7 @@ Where: - **GPU Acceleration**: Fully leverages NVIDIA CUDA for extreme-scale QP problems - **Flexible Problem Structure**: Supports sparse quadratic terms, low-rank quadratic terms, or both - **High Performance**: Competitive with commercial solvers on large-scale problems -- **SpMVOp Auto-Detection**: Automatically uses cuSPARSE SpMVOp (fused SpMV+elementwise operations) on CUDA 13+ while falling back to standard SpMV on CUDA 12.x +- **SpMVOp Auto-Detection**: Automatically uses cuSPARSE SpMVOp on CUDA 13+ while falling back to standard SpMV on CUDA 12.x - **Multi-GPU Distributed Solving**: Supports parallel solving across multiple GPUs via MPI and NCCL (optional, enabled at compile time) ## Quick Links diff --git a/docs/installation.md b/docs/installation.md index 60e4a76..70be8e6 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -10,7 +10,7 @@ !!! note "CUDA Version and SpMVOp" PDHCG automatically detects your CUDA version at compile time: - - **CUDA 13+**: Uses cuSPARSE **SpMVOp** (fused SpMV + elementwise operations) for improved performance. + - **CUDA 13+**: Uses cuSPARSE **SpMVOp** for improved performance. - **CUDA 12.x**: Falls back to the standard **SpMV** API. No manual intervention is required. ## C++ Executable diff --git a/python/README.md b/python/README.md index fd52d31..84ced3b 100644 --- a/python/README.md +++ b/python/README.md @@ -20,7 +20,7 @@ It provides a high-level, Pythonic API for constructing, modifying, and solving !!! note "CUDA Version and SpMVOp" PDHCG automatically detects your CUDA version at compile time: - - **CUDA 13+**: Uses cuSPARSE **SpMVOp** (fused SpMV + elementwise operations) for improved performance. + - **CUDA 13+**: Uses cuSPARSE **SpMVOp** for improved performance. - **CUDA 12.x**: Falls back to the standard **SpMV** API. No manual intervention is required. ### Install From d0b037c1285de5b2b7fc3a6ba2d42e75b76653a5 Mon Sep 17 00:00:00 2001 From: Lhongpei <1453244320@qq.com> Date: Fri, 22 May 2026 09:14:52 +0000 Subject: [PATCH 3/4] Support Jacobi Precondition in inner solving --- distributed/distributed_solver.cu | 12 ++- include/pdhcg_types.h | 1 + internal/internal_types.h | 7 ++ internal/pdhcg_kernels.cuh | 34 +++++++ src/cli.c | 11 +++ src/pdhcg_kernels.cu | 106 +++++++++++++++++++++ src/pdhg_core_op.cu | 147 +++++++++++++++++++++--------- src/presolve_wrapper.c | 42 +++++---- src/solver_state.cu | 68 +++++++++++++- src/utils.cu | 4 +- 10 files changed, 365 insertions(+), 67 deletions(-) diff --git a/distributed/distributed_solver.cu b/distributed/distributed_solver.cu index 98b5fa4..959dfaf 100644 --- a/distributed/distributed_solver.cu +++ b/distributed/distributed_solver.cu @@ -433,8 +433,16 @@ pdhcg_result_t *distributed_optimize(const pdhg_parameters_t *params, const qp_p row_perm = (int *)malloc(working_problem->num_constraints * sizeof(int)); col_perm = (int *)malloc(working_problem->num_variables * sizeof(int)); - generate_random_permutation(working_problem->num_constraints, row_perm); - generate_random_permutation(working_problem->num_variables, col_perm); + if (params->permute_method == FULL_RANDOM_PERMUTATION) + { + generate_random_permutation(working_problem->num_variables, col_perm); + generate_random_permutation(working_problem->num_constraints, row_perm); + } + else if (params->permute_method == BLOCK_RANDOM_PERMUTATION) + { + generate_block_permutation(working_problem->num_variables, params->permute_block_size, col_perm); + generate_block_permutation(working_problem->num_constraints, params->permute_block_size, row_perm); + } permuted_problem = permute_problem_return_new(working_problem, row_perm, col_perm); working_problem = permuted_problem; diff --git a/include/pdhcg_types.h b/include/pdhcg_types.h index 901ba03..38559de 100644 --- a/include/pdhcg_types.h +++ b/include/pdhcg_types.h @@ -151,6 +151,7 @@ extern "C" norm_type_t optimality_norm; inner_solver_parameters_t inner_solver_parameters; bool presolve; + bool diag_jacobi_precond; partition_method_t partition_method; permute_method_t permute_method; grid_size_t grid_size; diff --git a/internal/internal_types.h b/internal/internal_types.h index ba1bf57..28c1811 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -70,6 +70,13 @@ typedef struct double *gradient; double *direction; double *scalar_buffer; + bool precond_enabled; + double *diag_h_static; + double *m_diag; + double *m_inv; + double *Ms_buffer; + double cached_inv_tau; + double tol_scale; } bb_step_size_t; typedef struct diff --git a/internal/pdhcg_kernels.cuh b/internal/pdhcg_kernels.cuh index 69d27d5..787bea6 100644 --- a/internal/pdhcg_kernels.cuh +++ b/internal/pdhcg_kernels.cuh @@ -170,6 +170,18 @@ extern "C" __global__ void compute_bb_alpha_safeguard_kernel(const double *d_norm_gtg, const double *d_tmp, double *d_alpha); + __global__ void compute_bb_alpha_M_kernel(const double *d_stMs, const double *d_tmp, double *d_alpha); + + __global__ void scalar_sqrt_copy_kernel(const double *src, double *dst); + + __global__ void + compute_csr_diag_kernel(const int *row_ptr, const int *col_ind, const double *val, double *diag, int num_rows); + + __global__ void compute_csr_row_sq_norm_kernel(const int *row_ptr, const double *val, double *out, int num_rows); + + __global__ void refresh_inner_precond_kernel( + const double *diag_h_static, double inv_tau, double *m_diag, double *m_inv, int n_vars); + __global__ void primal_gradient_descent_kernel_bb_init(const double *dual_product, double *gradient, double *direction, @@ -200,6 +212,28 @@ extern "C" const double *d_alpha, const int n_vars); + __global__ void primal_gradient_descent_kernel_bb_init_precond(const double *dual_product, + double *gradient, + double *direction, + const double *current_primal_solution, + double *pdhg_primal_solution, + const double *objective_vector, + const double *objective_product, + const double *var_lb, + const double *var_ub, + const double *m_inv, + const double stepsize, + const int n_vars); + + __global__ void primal_bb_update_direction_kernel_precond(double *pdhg_primal_solution, + const double *gradient, + double *direction, + const double *var_lb, + const double *var_ub, + const double *m_inv, + const double *d_alpha, + const int n_vars); + __global__ void primal_bb_final_kernel(const double *current_primal_solution, const double *pdhg_primal_solution, double *reflected_primal_solution, diff --git a/src/cli.c b/src/cli.c index 63912a9..e7ef240 100644 --- a/src/cli.c +++ b/src/cli.c @@ -169,6 +169,9 @@ void print_usage(const char *prog_name) fprintf(stderr, " --inner_init_tol Initial tolerance for the inner solver (default: 1e-3).\n"); fprintf(stderr, " --inner_min_tol Minimum tolerance for the inner solver (default: 1e-9).\n"); fprintf(stderr, " --presolve Enable (1) or disable (0) presolve (default: 1).\n"); + fprintf( + stderr, + " --no_diag_precond Disable Jacobi diagonal preconditioner for inner subproblem (default: enabled).\n"); #ifdef PDHCG_COMPILE_DISTRIBUTED fprintf(stderr, "\nDistributed Options (MPI & NCCL):\n"); @@ -206,6 +209,7 @@ int run_pdhcg(int argc, char *argv[]) {"inner_init_tol", required_argument, 0, 1016}, {"inner_min_tol", required_argument, 0, 1017}, {"presolve", required_argument, 0, 1018}, + {"no_diag_precond", no_argument, 0, 1019}, {0, 0, 0, 0}}; int opt; @@ -285,6 +289,9 @@ int run_pdhcg(int argc, char *argv[]) case 1018: params.presolve = (atoi(optarg) != 0); break; + case 1019: + params.diag_jacobi_precond = false; + break; case '?': return 1; } @@ -376,6 +383,7 @@ int run_d_pdhcg(int argc, char *argv[]) {"inner_init_tol", required_argument, 0, 1016}, {"inner_min_tol", required_argument, 0, 1017}, {"presolve", required_argument, 0, 1018}, + {"no_diag_precond", no_argument, 0, 1019}, {"grid_size", required_argument, 0, 2001}, {"partition_method", required_argument, 0, 2002}, {"permute_method", required_argument, 0, 2003}, @@ -462,6 +470,9 @@ int run_d_pdhcg(int argc, char *argv[]) case 1018: params.presolve = (atoi(optarg) != 0); break; + case 1019: + params.diag_jacobi_precond = false; + break; case 2001: // --grid_size r,c { int r, c; diff --git a/src/pdhcg_kernels.cu b/src/pdhcg_kernels.cu index dfb9e42..ad35f14 100644 --- a/src/pdhcg_kernels.cu +++ b/src/pdhcg_kernels.cu @@ -277,6 +277,65 @@ __global__ void compute_bb_alpha_safeguard_kernel(const double *d_norm_gtg, cons *d_alpha = (*d_norm_gtg * *d_norm_gtg) / *d_tmp; } +__global__ void compute_bb_alpha_M_kernel(const double *d_stMs, const double *d_tmp, double *d_alpha) +{ + *d_alpha = *d_stMs / *d_tmp; +} + +__global__ void scalar_sqrt_copy_kernel(const double *src, double *dst) +{ + *dst = sqrt(*src); +} + +__global__ void +compute_csr_diag_kernel(const int *row_ptr, const int *col_ind, const double *val, double *diag, int num_rows) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < num_rows) + { + double sum = 0.0; + int start = row_ptr[i]; + int end = row_ptr[i + 1]; + for (int k = start; k < end; ++k) + { + if (col_ind[k] == i) + sum += val[k]; + } + diag[i] = sum; + } +} + +__global__ void compute_csr_row_sq_norm_kernel(const int *row_ptr, const double *val, double *out, int num_rows) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < num_rows) + { + double sum = 0.0; + int start = row_ptr[i]; + int end = row_ptr[i + 1]; + for (int k = start; k < end; ++k) + { + double v = val[k]; + sum += v * v; + } + out[i] = sum; + } +} + +__global__ void +refresh_inner_precond_kernel(const double *diag_h_static, double inv_tau, double *m_diag, double *m_inv, int n_vars) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n_vars) + { + double m = diag_h_static[i] + inv_tau; + if (m <= 0.0) + m = 1.0; + m_diag[i] = m; + m_inv[i] = 1.0 / m; + } +} + __global__ void primal_gradient_descent_kernel_bb_init(const double *dual_product, double *gradient, double *direction, @@ -357,6 +416,53 @@ __global__ void primal_bb_final_kernel(const double *current_primal_solution, } } +__global__ void primal_gradient_descent_kernel_bb_init_precond(const double *dual_product, + double *gradient, + double *direction, + const double *current_primal_solution, + double *pdhg_primal_solution, + const double *objective_vector, + const double *objective_product, + const double *var_lb, + const double *var_ub, + const double *m_inv, + const double stepsize, + const int n_vars) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n_vars) + { + double current_grad = objective_product[i] + objective_vector[i] - dual_product[i]; + double current_primal_sol = current_primal_solution[i]; + double next_primal_sol = current_primal_sol - stepsize * m_inv[i] * current_grad; + next_primal_sol = fmax(var_lb[i], fmin(next_primal_sol, var_ub[i])); + pdhg_primal_solution[i] = next_primal_sol; + gradient[i] = current_grad; + direction[i] = next_primal_sol - current_primal_sol; + } +} + +__global__ void primal_bb_update_direction_kernel_precond(double *pdhg_primal_solution, + const double *gradient, + double *direction, + const double *var_lb, + const double *var_ub, + const double *m_inv, + const double *d_alpha, + const int n_vars) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n_vars) + { + double alpha = *d_alpha; + double cur_sol = pdhg_primal_solution[i]; + double next_sol = cur_sol - alpha * m_inv[i] * gradient[i]; + next_sol = fmax(var_lb[i], fmin(next_sol, var_ub[i])); + direction[i] = next_sol - cur_sol; + pdhg_primal_solution[i] = next_sol; + } +} + __global__ void compute_lp_residual_kernel(double *primal_residual, const double *primal_product, const double *constraint_lower_bound, diff --git a/src/pdhg_core_op.cu b/src/pdhg_core_op.cu index 5f914b3..3145b11 100644 --- a/src/pdhg_core_op.cu +++ b/src/pdhg_core_op.cu @@ -216,23 +216,65 @@ void primal_BB_step_size_update(pdhg_solver_state_t *state, double step_size) int inner_solver_iter = 1; double initial_alpha = 1.0 / inv_step_size; - double *d_norm_gtg = state->inner_solver->bb_step_size->scalar_buffer; - double *d_tmp = state->inner_solver->bb_step_size->scalar_buffer + 1; - double *d_alpha = state->inner_solver->bb_step_size->scalar_buffer + 2; + bb_step_size_t *bb = state->inner_solver->bb_step_size; + bool precond = bb->precond_enabled; + + double *d_norm_gtg = bb->scalar_buffer; + double *d_tmp = bb->scalar_buffer + 1; + double *d_alpha = bb->scalar_buffer + 2; + double *d_stMs = bb->scalar_buffer + 3; + + if (precond && bb->cached_inv_tau != inv_step_size) + { + refresh_inner_precond_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + bb->diag_h_static, inv_step_size, bb->m_diag, bb->m_inv, state->num_variables); + bb->cached_inv_tau = inv_step_size; + + double sum_m = 0.0; + CUBLAS_CHECK(cublasDasum(state->blas_handle, state->num_variables, bb->m_diag, 1, &sum_m)); + pdhcg_all_reduce_scalar(state->grid_context, &sum_m, PDHCG_OP_SUM, PDHCG_SCOPE_ROW, false); + int n_global = get_global_n(state); + if (n_global > 0) + bb->tol_scale = sqrt(sum_m / (double)n_global); + else + bb->tol_scale = 1.0; + } + + if (precond) + initial_alpha *= bb->tol_scale * bb->tol_scale; update_obj_product(state, state->current_primal_solution); - primal_gradient_descent_kernel_bb_init<<num_blocks_primal, THREADS_PER_BLOCK>>>( - state->dual_product, - state->inner_solver->bb_step_size->gradient, - state->inner_solver->bb_step_size->direction, - state->current_primal_solution, - state->pdhg_primal_solution, - state->objective_vector, - state->quadratic_objective_term->primal_obj_product, - state->variable_lower_bound, - state->variable_upper_bound, - initial_alpha, - state->num_variables); + if (precond) + { + primal_gradient_descent_kernel_bb_init_precond<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->dual_product, + bb->gradient, + bb->direction, + state->current_primal_solution, + state->pdhg_primal_solution, + state->objective_vector, + state->quadratic_objective_term->primal_obj_product, + state->variable_lower_bound, + state->variable_upper_bound, + bb->m_inv, + initial_alpha, + state->num_variables); + } + else + { + primal_gradient_descent_kernel_bb_init<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->dual_product, + bb->gradient, + bb->direction, + state->current_primal_solution, + state->pdhg_primal_solution, + state->objective_vector, + state->quadratic_objective_term->primal_obj_product, + state->variable_lower_bound, + state->variable_upper_bound, + initial_alpha, + state->num_variables); + } cublasSetPointerMode(state->blas_handle, CUBLAS_POINTER_MODE_DEVICE); @@ -241,17 +283,22 @@ void primal_BB_step_size_update(pdhg_solver_state_t *state, double step_size) while (inner_solver_iter < state->inner_solver->iteration_limit) { - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->inner_solver->bb_step_size->direction, - 1, - state->inner_solver->bb_step_size->direction, - 1, - d_norm_gtg)); - - pdhcg_all_reduce_scalar(state->grid_context, d_norm_gtg, PDHCG_OP_SUM, PDHCG_SCOPE_ROW, true); - - sqrt_scalar_kernel<<<1, 1>>>(d_norm_gtg); + if (precond) + { + element_wise_mul_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + bb->m_diag, bb->direction, bb->Ms_buffer, state->num_variables); + CUBLAS_CHECK( + cublasDdot(state->blas_handle, state->num_variables, bb->direction, 1, bb->Ms_buffer, 1, d_stMs)); + pdhcg_all_reduce_scalar(state->grid_context, d_stMs, PDHCG_OP_SUM, PDHCG_SCOPE_ROW, true); + scalar_sqrt_copy_kernel<<<1, 1>>>(d_stMs, d_norm_gtg); + } + else + { + CUBLAS_CHECK( + cublasDdot(state->blas_handle, state->num_variables, bb->direction, 1, bb->direction, 1, d_norm_gtg)); + pdhcg_all_reduce_scalar(state->grid_context, d_norm_gtg, PDHCG_OP_SUM, PDHCG_SCOPE_ROW, true); + sqrt_scalar_kernel<<<1, 1>>>(d_norm_gtg); + } if (inner_solver_iter == 1 || inner_solver_iter % check_frequency == 0) { @@ -267,31 +314,43 @@ void primal_BB_step_size_update(pdhg_solver_state_t *state, double step_size) state->objective_vector, state->dual_product, state->quadratic_objective_term->primal_obj_product, - state->inner_solver->bb_step_size->gradient, + bb->gradient, state->inner_solver->primal_buffer, inv_step_size, state->num_variables); - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->inner_solver->bb_step_size->direction, - 1, - state->inner_solver->primal_buffer, - 1, - d_tmp)); + CUBLAS_CHECK(cublasDdot( + state->blas_handle, state->num_variables, bb->direction, 1, state->inner_solver->primal_buffer, 1, d_tmp)); pdhcg_all_reduce_scalar(state->grid_context, d_tmp, PDHCG_OP_SUM, PDHCG_SCOPE_ROW, true); - compute_bb_alpha_safeguard_kernel<<<1, 1>>>(d_norm_gtg, d_tmp, d_alpha); - - primal_bb_update_direction_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( - state->pdhg_primal_solution, - state->inner_solver->bb_step_size->gradient, - state->inner_solver->bb_step_size->direction, - state->variable_lower_bound, - state->variable_upper_bound, - d_alpha, - state->num_variables); + if (precond) + { + compute_bb_alpha_M_kernel<<<1, 1>>>(d_stMs, d_tmp, d_alpha); + + primal_bb_update_direction_kernel_precond<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->pdhg_primal_solution, + bb->gradient, + bb->direction, + state->variable_lower_bound, + state->variable_upper_bound, + bb->m_inv, + d_alpha, + state->num_variables); + } + else + { + compute_bb_alpha_safeguard_kernel<<<1, 1>>>(d_norm_gtg, d_tmp, d_alpha); + + primal_bb_update_direction_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->pdhg_primal_solution, + bb->gradient, + bb->direction, + state->variable_lower_bound, + state->variable_upper_bound, + d_alpha, + state->num_variables); + } inner_solver_iter++; } diff --git a/src/presolve_wrapper.c b/src/presolve_wrapper.c index c141de3..4098085 100644 --- a/src/presolve_wrapper.c +++ b/src/presolve_wrapper.c @@ -122,7 +122,8 @@ pdhcg_presolve_info_t *pdhcg_presolve(const qp_problem_t *original_prob, const p return NULL; info->settings = default_settings(); - ((Settings *)info->settings)->verbose = false; + ((Settings *)info->settings)->verbose = true; + ((Settings *)info->settings)->dual_fix = false; bool has_q = (original_prob->objective_sparse_matrix != NULL); bool has_r = (original_prob->objective_lowrank_matrix != NULL); @@ -165,24 +166,29 @@ pdhcg_presolve_info_t *pdhcg_presolve(const qp_problem_t *original_prob, const p } else if (has_q) { - size_t Pnnz = (size_t)original_prob->objective_sparse_matrix_num_nonzeros; + size_t Qnnz = (size_t)original_prob->objective_sparse_matrix_num_nonzeros; presolver = - new_qp_presolver(original_prob->constraint_matrix ? original_prob->constraint_matrix->val : NULL, - original_prob->constraint_matrix ? original_prob->constraint_matrix->col_ind : NULL, - original_prob->constraint_matrix ? original_prob->constraint_matrix->row_ptr : NULL, - m, - n, - nnz, - original_prob->constraint_lower_bound, - original_prob->constraint_upper_bound, - original_prob->variable_lower_bound, - original_prob->variable_upper_bound, - original_prob->objective_vector, - original_prob->objective_sparse_matrix->val, - original_prob->objective_sparse_matrix->col_ind, - original_prob->objective_sparse_matrix->row_ptr, - Pnnz, - info->settings); + new_qp_presolver_qr(original_prob->constraint_matrix ? original_prob->constraint_matrix->val : NULL, + original_prob->constraint_matrix ? original_prob->constraint_matrix->col_ind : NULL, + original_prob->constraint_matrix ? original_prob->constraint_matrix->row_ptr : NULL, + m, + n, + nnz, + original_prob->constraint_lower_bound, + original_prob->constraint_upper_bound, + original_prob->variable_lower_bound, + original_prob->variable_upper_bound, + original_prob->objective_vector, + original_prob->objective_sparse_matrix->val, + original_prob->objective_sparse_matrix->col_ind, + original_prob->objective_sparse_matrix->row_ptr, + Qnnz, + NULL, // no low-rank component + NULL, + NULL, + 0, + 0, + info->settings); } else { diff --git a/src/solver_state.cu b/src/solver_state.cu index fc4052a..bbd1309 100644 --- a/src/solver_state.cu +++ b/src/solver_state.cu @@ -17,6 +17,7 @@ limitations under the License. #include "internal_types.h" #include "pdhcg.h" +#include "pdhcg_kernels.cuh" #include "pdhg_core_op.h" #include "preconditioner.h" #include "solver.h" @@ -317,14 +318,49 @@ static void initialize_inner_solver(pdhg_solver_state_t *state, const pdhg_param case PDHCG_LOW_RANK_Q: case PDHCG_LOW_RANK_PLUS_SPARSE_Q: state->inner_solver->has_inner_loop = true; - state->inner_solver->bb_step_size = (bb_step_size_t *)safe_malloc(sizeof(bb_step_size_t)); + state->inner_solver->bb_step_size = (bb_step_size_t *)safe_calloc(1, sizeof(bb_step_size_t)); ALLOC_ZERO(state->inner_solver->bb_step_size->gradient, state->num_variables * sizeof(double)); ALLOC_ZERO(state->inner_solver->bb_step_size->direction, state->num_variables * sizeof(double)); - ALLOC_ZERO(state->inner_solver->bb_step_size->scalar_buffer, 3 * sizeof(double)) + ALLOC_ZERO(state->inner_solver->bb_step_size->scalar_buffer, 4 * sizeof(double)); state->inner_solver->iteration_limit = iteration_limit; state->inner_solver->tol = initial_tol; state->inner_solver->min_tol = min_tol; + + state->inner_solver->bb_step_size->precond_enabled = params->diag_jacobi_precond; + if (params->diag_jacobi_precond) + { + int n = state->num_variables; + ALLOC_ZERO(state->inner_solver->bb_step_size->diag_h_static, n * sizeof(double)); + ALLOC_ZERO(state->inner_solver->bb_step_size->m_diag, n * sizeof(double)); + ALLOC_ZERO(state->inner_solver->bb_step_size->m_inv, n * sizeof(double)); + ALLOC_ZERO(state->inner_solver->bb_step_size->Ms_buffer, n * sizeof(double)); + state->inner_solver->bb_step_size->cached_inv_tau = -1.0; + state->inner_solver->bb_step_size->tol_scale = 1.0; + + if (state->quadratic_objective_term->quad_obj_type == PDHCG_SPARSE_Q || + state->quadratic_objective_term->quad_obj_type == PDHCG_LOW_RANK_PLUS_SPARSE_Q) + { + cu_sparse_matrix_csr_t *Q = state->quadratic_objective_term->objective_sparse_matrix; + compute_csr_diag_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + Q->row_ptr, Q->col_ind, Q->val, state->inner_solver->bb_step_size->diag_h_static, n); + CUDA_CHECK(cudaGetLastError()); + } + + if (state->quadratic_objective_term->quad_obj_type == PDHCG_LOW_RANK_Q || + state->quadratic_objective_term->quad_obj_type == PDHCG_LOW_RANK_PLUS_SPARSE_Q) + { + cu_sparse_matrix_csr_t *Rt = state->quadratic_objective_term->objective_lowrank_matrix_t; + double *out = state->inner_solver->bb_step_size->Ms_buffer; + compute_csr_row_sq_norm_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + Rt->row_ptr, Rt->val, out, n); + CUDA_CHECK(cudaGetLastError()); + const double one = 1.0; + CUBLAS_CHECK(cublasDaxpy( + state->blas_handle, n, &one, out, 1, state->inner_solver->bb_step_size->diag_h_static, 1)); + CUDA_CHECK(cudaMemset(out, 0, n * sizeof(double))); + } + } break; default: fprintf(stderr, "Error: Unknown Quadratic Objective Type detected.\n"); @@ -807,6 +843,34 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) free(state->quadratic_objective_term); } + if (state->inner_solver) + { + if (state->inner_solver->bb_step_size) + { + bb_step_size_t *bb = state->inner_solver->bb_step_size; + if (bb->gradient) + CUDA_CHECK(cudaFree(bb->gradient)); + if (bb->direction) + CUDA_CHECK(cudaFree(bb->direction)); + if (bb->scalar_buffer) + CUDA_CHECK(cudaFree(bb->scalar_buffer)); + if (bb->diag_h_static) + CUDA_CHECK(cudaFree(bb->diag_h_static)); + if (bb->m_diag) + CUDA_CHECK(cudaFree(bb->m_diag)); + if (bb->m_inv) + CUDA_CHECK(cudaFree(bb->m_inv)); + if (bb->Ms_buffer) + CUDA_CHECK(cudaFree(bb->Ms_buffer)); + free(bb); + } + if (state->inner_solver->primal_buffer) + CUDA_CHECK(cudaFree(state->inner_solver->primal_buffer)); + if (state->inner_solver->dual_buffer) + CUDA_CHECK(cudaFree(state->inner_solver->dual_buffer)); + free(state->inner_solver); + } + free(state); } diff --git a/src/utils.cu b/src/utils.cu index 94e3478..e617234 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -338,8 +338,10 @@ void set_default_parameters(pdhg_parameters_t *params) params->inner_solver_parameters.min_tolerance = 1e-9; params->grid_size.decided = false; - params->partition_method = UNIFORM_PARTITION; + params->partition_method = NNZ_BALANCE_PARTITION; params->permute_method = BLOCK_RANDOM_PERMUTATION; + + params->diag_jacobi_precond = true; } #define PRINT_DIFF_INT(name, current, default_val) \ From 2a832e264fd9fa2924a8f9714d79f5840bfb05ed Mon Sep 17 00:00:00 2001 From: Lhongpei <1453244320@qq.com> Date: Fri, 22 May 2026 09:27:07 +0000 Subject: [PATCH 4/4] add jacobi doc --- README.md | 1 + docs/c/types.md | 1 + docs/installation.md | 3 +++ docs/python/parameters.md | 1 + python/README.md | 4 ++++ python/pdhcg/PDHCG.py | 1 + python_bindings/_core_bindings.cpp | 2 ++ 7 files changed, 13 insertions(+) diff --git a/README.md b/README.md index 554b7e6..cad0353 100644 --- a/README.md +++ b/README.md @@ -99,6 +99,7 @@ Solver Parameters: | --inner_init_tol | double | Initial tolerance for the inner solver. | 1e-3 | | --inner_min_tol | double | Minimum tolerance for the inner solver. | 1e-9 | | --presolve | int | Enable (1) or disable (0) presolve. | 1 | +| --no_diag_precond | flag | Disable the Jacobi diagonal preconditioner used in the inner subproblem (enabled by default). | false | **Distributed Options** (only available when built with `-DPDHCG_COMPILE_DISTRIBUTED=ON`): | Option | Type | Description | Default | diff --git a/docs/c/types.md b/docs/c/types.md index 831a4f2..d83b9c0 100644 --- a/docs/c/types.md +++ b/docs/c/types.md @@ -235,6 +235,7 @@ typedef struct { norm_type_t optimality_norm; inner_solver_parameters_t inner_solver_parameters; bool presolve; + bool diag_jacobi_precond; partition_method_t partition_method; permute_method_t permute_method; grid_size_t grid_size; diff --git a/docs/installation.md b/docs/installation.md index 70be8e6..2e1b2e7 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -51,6 +51,9 @@ When enabled, the solver binary automatically detects whether it is launched wit ## Python Package +!!! note "Multi-GPU support" + The Python interface currently supports single-GPU solving only. For multi-GPU distributed solving, build the C++ executable with `-DPDHCG_COMPILE_DISTRIBUTED=ON` and launch it via `mpirun`. + ### From PyPI (Recommended) ```bash diff --git a/docs/python/parameters.md b/docs/python/parameters.md index af00fd9..6e6238f 100644 --- a/docs/python/parameters.md +++ b/docs/python/parameters.md @@ -40,6 +40,7 @@ m.setParams(TimeLimit=3600, LogLevel=1) | `InnerIterLimit` | int | 1000 | Max iterations for inner CG solver | | `InnerInitTol` | float | 1e-3 | Initial tolerance for inner solver | | `InnerMinTol` | float | 1e-9 | Minimum tolerance for inner solver | +| `DiagJacobiPrecond` | bool | True | Use the Jacobi diagonal preconditioner in the inner subproblem. Set to `False` to disable. | ### Singular Value Estimation diff --git a/python/README.md b/python/README.md index 84ced3b..cf94550 100644 --- a/python/README.md +++ b/python/README.md @@ -23,6 +23,9 @@ It provides a high-level, Pythonic API for constructing, modifying, and solving - **CUDA 13+**: Uses cuSPARSE **SpMVOp** for improved performance. - **CUDA 12.x**: Falls back to the standard **SpMV** API. No manual intervention is required. +!!! note "Multi-GPU support" + The Python interface currently supports single-GPU solving only. For multi-GPU distributed solving, build the C++ executable with `-DPDHCG_COMPILE_DISTRIBUTED=ON` and launch it via `mpirun` (see the [main README](../README.md)). + ### Install Install from PyPI: @@ -189,6 +192,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `InnerIterLimit` | `inner_iter_limit` | int | `1000` | Maximum number of iterations for the inner solver. | | `InnerInitTol` | `inner_init_tol` | float | `1e-3` | Initial tolerance for the inner solver. | | `InnerMinTol` | `inner_min_tol` | float | `1e-9` | Minimum tolerance for the inner solver. | +| `DiagJacobiPrecond` | `diag_jacobi_precond` | bool | `True` | Whether to use the Jacobi diagonal preconditioner in the inner subproblem. Set to `False` to disable. | They can be set in multiple ways: diff --git a/python/pdhcg/PDHCG.py b/python/pdhcg/PDHCG.py index ca2a9dc..3057211 100644 --- a/python/pdhcg/PDHCG.py +++ b/python/pdhcg/PDHCG.py @@ -61,6 +61,7 @@ "InnerIterLimit": "inner_iter_limit", "InnerInitTol": "inner_init_tol", "InnerMinTol": "inner_min_tol", + "DiagJacobiPrecond": "diag_jacobi_precond", # presolve "Presolve": "presolve", } diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index 6d4cdcd..1919a42 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -296,6 +296,7 @@ static py::dict get_default_params_py() d["inner_iter_limit"] = p.inner_solver_parameters.iteration_limit; d["inner_init_tol"] = p.inner_solver_parameters.initial_tolerance; d["inner_min_tol"] = p.inner_solver_parameters.min_tolerance; + d["diag_jacobi_precond"] = p.diag_jacobi_precond; return d; } @@ -381,6 +382,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p geti("inner_iter_limit", p->inner_solver_parameters.iteration_limit); getf("inner_init_tol", p->inner_solver_parameters.initial_tolerance); getf("inner_min_tol", p->inner_solver_parameters.min_tolerance); + getb("diag_jacobi_precond", p->diag_jacobi_precond); // presolve getb("presolve", p->presolve);