From b6e58dbd8b905b31ebd6a68f15ec7949d6ed67f1 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Sun, 20 Sep 2026 02:10:22 +0000 Subject: [PATCH 01/15] improve infeasibility detection --- CMakeLists.txt | 7 + README.md | 6 +- include/cupdlpx_types.h | 1 + internal/cuda_to_hip.h | 7 + internal/infeasibility.h | 51 ++ internal/internal_types.h | 88 ++- internal/reduce.h | 81 +++ internal/termination.h | 82 +++ internal/utils.h | 27 - python/README.md | 3 +- python/cupdlpx/PDLP.py | 1 + python/cupdlpx/model.py | 5 +- python_bindings/_core_bindings.cpp | 7 +- src/cli.c | 20 + src/feasibility_polish.cu | 271 +++++++--- src/infeasibility.cu | 589 ++++++++++++++++++++ src/reduce.cu | 51 ++ src/solver.cu | 211 +++----- src/termination.cu | 285 ++++++++++ src/utils.cu | 632 +--------------------- test/test_infeasibility_false_positive.py | 117 ++++ test/test_infeasible_unbounded.py | 28 +- test/test_interface.c | 2 +- 23 files changed, 1722 insertions(+), 850 deletions(-) create mode 100644 internal/infeasibility.h create mode 100644 internal/reduce.h create mode 100644 internal/termination.h create mode 100644 src/infeasibility.cu create mode 100644 src/reduce.cu create mode 100644 src/termination.cu create mode 100644 test/test_infeasibility_false_positive.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 44625f7..0c27288 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -70,6 +70,11 @@ endif() # ----------------------------------------------------------------------------- add_library(cupdlpx_compile_flags INTERFACE) +if(NOT USE_HIP) + # reduce.h passes device lambdas to thrust::transform_iterator + target_compile_options(cupdlpx_compile_flags INTERFACE $<$:--extended-lambda>) +endif() + if(MSVC) target_compile_options(cupdlpx_compile_flags INTERFACE $<$,$>:/O2 /W4 /Zi> @@ -126,6 +131,7 @@ if(USE_HIP) find_package(hipsparse REQUIRED) find_package(hipcub REQUIRED) find_package(rocprim REQUIRED) + find_package(rocthrust REQUIRED) endif() include(FetchContent) @@ -206,6 +212,7 @@ if(USE_HIP) PUBLIC roc::hipblas PUBLIC roc::hipsparse PUBLIC hip::hipcub + PUBLIC roc::rocthrust PUBLIC ZLIB::ZLIB PUBLIC PSLP ) diff --git a/README.md b/README.md index 517fdaa..98cdc4d 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,7 @@ This will create the solver binary at `./build/cupdlpx`. #### Building for AMD GPUs (ROCm/HIP) To target AMD GPUs, configure with `-DUSE_HIP=ON` and select the GPU architecture with `-DCMAKE_HIP_ARCHITECTURES`. The CUDA sources are compiled -as HIP and the cuBLAS/cuSPARSE/CUB calls are mapped to hipBLAS/hipSPARSE/hipCUB. +as HIP and the cuBLAS/cuSPARSE/CUB/Thrust calls are mapped to hipBLAS/hipSPARSE/hipCUB/rocThrust. ```bash cmake -B build -DUSE_HIP=ON -DCMAKE_HIP_ARCHITECTURES=gfx90a -DCMAKE_PREFIX_PATH=/opt/rocm cmake --build build --clean-first @@ -57,7 +57,7 @@ Set `CMAKE_HIP_ARCHITECTURES` to match your GPU (for example `gfx90a` for MI200, `gfx1100` for RDNA3 desktop, or `gfx1201` for RDNA4). If the ROCm install is not on CMake's default search path, point `-DCMAKE_PREFIX_PATH` at it (e.g. `/opt/rocm`) so `find_package` can locate hip, hipBLAS, hipSPARSE, -and hipCUB. The resulting `./build/cupdlpx` binary is used exactly as in the +hipCUB, and rocThrust. The resulting `./build/cupdlpx` binary is used exactly as in the CUDA build. #### Verifying the Installation @@ -117,6 +117,8 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` | | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | | `--no_presolve` | `flag` | Disable presolve | `enabled` | +| `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate | `1e-8` | +| `--no_infeasibility_detection` | `flag` | Disable infeasibility detection | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 16d3f04..25b9386 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -108,6 +108,7 @@ extern "C" bool feasibility_polishing; norm_type_t optimality_norm; bool presolve; + bool infeasibility_detection; double matrix_zero_tol; double infinite_bound; int geometric_mean_iterations; diff --git a/internal/cuda_to_hip.h b/internal/cuda_to_hip.h index fc9975e..4b230f7 100644 --- a/internal/cuda_to_hip.h +++ b/internal/cuda_to_hip.h @@ -154,6 +154,13 @@ static inline const char *cublasGetStatusName(hipblasStatus_t status) #define cusparseSpMV hipsparseSpMV #define cusparseSpMV_bufferSize hipsparseSpMV_bufferSize #define cusparseSpMV_preprocess hipsparseSpMV_preprocess +#define cusparseDnMatDescr_t hipsparseDnMatDescr_t +#define cusparseCreateDnMat hipsparseCreateDnMat +#define cusparseDestroyDnMat hipsparseDestroyDnMat +#define cusparseSpMM hipsparseSpMM +#define cusparseSpMM_bufferSize hipsparseSpMM_bufferSize +#define CUSPARSE_ORDER_COL HIPSPARSE_ORDER_COL +#define CUSPARSE_SPMM_ALG_DEFAULT HIPSPARSE_SPMM_ALG_DEFAULT #define cusparseCsr2cscEx2 hipsparseCsr2cscEx2 #define cusparseCsr2cscEx2_bufferSize hipsparseCsr2cscEx2_bufferSize diff --git a/internal/infeasibility.h b/internal/infeasibility.h new file mode 100644 index 0000000..acf2001 --- /dev/null +++ b/internal/infeasibility.h @@ -0,0 +1,51 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#pragma once + +#include "utils.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + + // Detection of primal and dual infeasibility from candidate rays; see the + // comment above compute_infeasibility_device in infeasibility.cu for the + // method. Same protocol as termination.h: compute_infeasibility_device + // only queues kernels (the candidates, A d, A^T y and their reductions) and + // never waits; the caller fetches the scalars; compute_infeasibility_host + // turns them into the certificate errors of the state. + // + // A check has a second round, because |A||d| and |A^T||y| cost two more + // SpMMs and are only needed when some candidate passes its objective + // margin. compute_infeasibility_host runs it when needed, with one more + // copy and synchronization; otherwise it is host arithmetic only. + + // The context owns every buffer of a check. It is created once before the + // main loop, so that a check never allocates, and destroyed with the solver + // state. + infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state); + + void infeasibility_context_destroy(infeas_ctx_t *ctx); + + void compute_infeasibility_device(pdhg_solver_state_t *state); + + void compute_infeasibility_host(pdhg_solver_state_t *state); + +#ifdef __cplusplus +} +#endif diff --git a/internal/internal_types.h b/internal/internal_types.h index bde5a4e..fe73a3a 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -25,6 +25,75 @@ limitations under the License. #include #include +// The scalars a termination check produces on the device, copied to the host +// in one go; see pdhg_solver_state_t::scalars_h. All of them are in the +// rescaled space; the host halves unscale. Name parts: res = residual, +// sq = sum of squares, max = infinity norm, dot = inner product. +// +// Each member is the destination of one reduction (reduce.h). + +#define INFEAS_NUM_CANDIDATES 3 // candidate rays of one infeasibility check + +// One candidate of the infeasibility check (infeasibility.cu): a primal ray d, +// which certifies dual infeasibility, and a dual ray y with reduced cost +// r = -A^T y, which certifies primal infeasibility. Every norm here is the +// infinity norm: the certificate criterion is defined with it, and the +// reported statistics follow PDLP's normalization. +// +// The violation exists twice because the two are different maxima over the +// rows. The scaled one is the numerator of the certificate error +// viol_max_scaled / product_abs_max, evaluated in the rescaled space because +// the ratio is not invariant to non-uniform row scaling and the rescaling has +// equilibrated the rows. The orig one is in the units of the original problem +// and is only reported (max_primal_ray_infeasibility and its dual twin). +typedef struct +{ + double max_abs; // ||d||_inf in the original space; normalizes the reported statistics + double obj; // c^T d, the objective term + double obj_abs; // |c|^T |d|, the same with absolute values + double viol_max_scaled; // max violation of the sign conditions on A d + double product_abs_max; // max_i (|A||d|)_i, scaled like the violation above (second round) + double viol_max_orig; // the violation in original units +} primal_ray_scalars_t; + +typedef struct +{ + double max_abs; // ||y||_inf in the original space + double obj_con_bounds; // l_c^T y^+ - u_c^T y^-, the constraint-bound part of the dual objective term + double obj_var_bounds; // l_v^T r^+ - u_v^T r^-, the variable-bound part + double obj_con_bounds_abs; // |l_c|^T y^+ + |u_c|^T y^-, the constraint-bound part with absolute values + double obj_var_bounds_abs; // |l_v|^T r^+ + |u_v|^T r^-, the variable-bound part with absolute values + double viol_max_scaled; // max violation of the sign conditions on r + double product_abs_max; // max_j (|A^T||y|)_j, scaled like the violation above (second round) + double viol_max_orig; // the violation in original units +} dual_ray_scalars_t; + +typedef struct +{ + primal_ray_scalars_t primal_ray; + dual_ray_scalars_t dual_ray; +} infeas_candidate_scalars_t; + +typedef struct +{ + double fixed_point_dual_sq; // ||y - T(z)_y||^2, dual part of the fixed point error + double fixed_point_primal_sq; // ||x - T(z)_x||^2, primal part + double fixed_point_cross; // (A^T dy)^T dx, cross term of the same error + double primal_res_sq; // sum of squares of the primal residual + double primal_res_max; // infinity norm of the primal residual + double primal_slack_sum; // sum of the primal slacks, a dual objective term + double dual_res_sq; // sum of squares of the dual residual + double dual_res_max; // infinity norm of the dual residual + double objective_dot; // c^T x, the primal objective + double dual_slack_dot; // (dual slack)^T x, the other dual objective term + double restart_primal_sq; // ||x - x_anchor||^2, the primal weight's PID input + double restart_dual_sq; // ||y - y_anchor||^2, the same on the dual side + infeas_candidate_scalars_t infeas[INFEAS_NUM_CANDIDATES]; +} device_scalars_t; + +// Private state of the infeasibility detection, defined in infeasibility.cu. +typedef struct infeas_ctx infeas_ctx_t; + typedef struct { int num_rows; @@ -106,6 +175,8 @@ typedef struct double max_dual_ray_infeasibility; double primal_ray_linear_objective; double dual_ray_objective; + double primal_ray_certificate_error; + double dual_ray_certificate_error; termination_reason_t termination_reason; double *delta_primal_solution; @@ -119,9 +190,20 @@ typedef struct cusparseHandle_t sparse_handle; cublasHandle_t blas_handle; void *spmv_ctx; - - double *ones_primal_d; - double *ones_dual_d; + infeas_ctx_t *infeas_ctx; // detection buffers, created before the main loop when detection is on + + // Scalars of a termination check (optimality, infeasibility) and of the + // restart; see device_scalars_t above. Every quantity is produced on the + // device by the reductions of reduce.h and the whole struct is copied to + // the host in one go (fetch_device_scalars in termination.cu), so a check + // costs one host round trip instead of one per scalar. Both norms of each residual are always computed, so that L2 + // or L-inf stays a host-side choice. + device_scalars_t *scalars_d; + device_scalars_t *scalars_h; // pinned, so that the copy is asynchronous + void *reduce_temp_d; // temporary storage of the reductions + size_t reduce_temp_bytes; + + double *delta_dual_product; // A^T delta_dual_solution, scratch of the fixed point error double feasibility_polishing_time; int feasibility_iteration; diff --git a/internal/reduce.h b/internal/reduce.h new file mode 100644 index 0000000..e1b0a49 --- /dev/null +++ b/internal/reduce.h @@ -0,0 +1,81 @@ +// Copyright 2025 Haihao Lu +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Reductions of a termination check. C++ only. +// +// Every scalar of a check is one cub::DeviceReduce call over a +// thrust::transform_iterator: a device lambda computes the quantity from the +// element index on the fly, and the result is written to a member of the +// device-side device_scalars_t (internal_types.h), never to the host. The +// calls only queue kernels and never wait; the whole struct is copied to the +// host once per check (termination.h). +// DeviceReduce uses no atomics ("run-to-run" determinism), so a given problem +// size always produces the same result. +// +// The lambdas capture device pointers by value from local variables (never +// the solver state itself) and are marked __host__ __device__ so that the +// iterator can deduce their result type on the host; nvcc needs +// --extended-lambda for that (set in CMakeLists.txt). + +#ifndef CUPDLPX_REDUCE_H +#define CUPDLPX_REDUCE_H + +#include "utils.h" + +#include +#include + +// Temporary storage of cub::DeviceReduce, sized once for the longest vector of +// the problem and shared by every call (reduce.cu). +void reduce_temp_storage_create(pdhg_solver_state_t *state); + +// Exits with a message when a call would need more temporary storage than +// reduce_temp_storage_create allocated. +void reduce_check_temp_storage(const pdhg_solver_state_t *state, size_t bytes); + +// out = sum over i < n of f(i), on the stream, into device memory. +template void reduce_sum(pdhg_solver_state_t *state, int n, F f, double *out) +{ + if (n <= 0) + { + CUDA_CHECK(cudaMemsetAsync(out, 0, sizeof(double), state->stream)); + return; + } + auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); + size_t bytes = 0; + CUDA_CHECK(cub::DeviceReduce::Sum(NULL, bytes, it, out, n, state->stream)); + reduce_check_temp_storage(state, bytes); + bytes = state->reduce_temp_bytes; + CUDA_CHECK(cub::DeviceReduce::Sum(state->reduce_temp_d, bytes, it, out, n, state->stream)); +} + +// out = max over i < n of f(i), on the stream, into device memory. Every +// maximum the solver takes is a norm or a violation, hence nonnegative, so an +// empty vector yields 0. +template void reduce_max(pdhg_solver_state_t *state, int n, F f, double *out) +{ + if (n <= 0) + { + CUDA_CHECK(cudaMemsetAsync(out, 0, sizeof(double), state->stream)); + return; + } + auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); + size_t bytes = 0; + CUDA_CHECK(cub::DeviceReduce::Max(NULL, bytes, it, out, n, state->stream)); + reduce_check_temp_storage(state, bytes); + bytes = state->reduce_temp_bytes; + CUDA_CHECK(cub::DeviceReduce::Max(state->reduce_temp_d, bytes, it, out, n, state->stream)); +} + +#endif // CUPDLPX_REDUCE_H diff --git a/internal/termination.h b/internal/termination.h new file mode 100644 index 0000000..b72a894 --- /dev/null +++ b/internal/termination.h @@ -0,0 +1,82 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#pragma once + +#include "utils.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + + // The termination check: fixed point error, residuals and objective values + // of the current iterate, and the criteria that decide whether the solver + // stops. The infeasibility check (infeasibility.h) follows the same + // protocol. + // + // Every scalar of a check is produced on the device by the reductions of + // reduce.h and lands in state->scalars_d. A check is computed in halves so + // that it synchronizes once: a _device half only queues kernels on the + // stream and never waits, so the host can queue a whole check while the + // GPU is still running the iterations before it; a _host half only does + // host arithmetic on state->scalars_h. In between, the caller calls + // fetch_device_scalars. The functions without a suffix run a whole check, + // synchronization included. + + // Copies the whole scalar struct to the host and waits for the stream: + // the one synchronization of a check. + void fetch_device_scalars(pdhg_solver_state_t *state); + + // state->delta_primal_solution = to_primal - from_primal, and the same for + // the dual vectors. + void compute_solution_difference(pdhg_solver_state_t *state, + const double *from_primal, + const double *to_primal, + const double *from_dual, + const double *to_dual); + + // Residual vectors and their reductions from the products already in the + // state (A x in primal_product, A^T y in dual_product): the primal + // residual, the dual residual, the primal slack (the contribution of y to + // the dual objective against the finite constraint bounds given), c^T x + // and (dual slack)^T x. The optimality check passes the state's own + // vectors; feasibility polishing, whose states carry a modified problem, + // passes the original problem's objective vector, primal solution and + // finite bounds where its formulas need them. + void launch_residual_check(pdhg_solver_state_t *state, + const double *objective_vector, + const double *primal_solution, + const double *constraint_lower_bound_finite, + const double *constraint_upper_bound_finite); + + void compute_residual_device(pdhg_solver_state_t *state); + + void compute_residual_host(pdhg_solver_state_t *state, norm_type_t optimality_norm); + + void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm); + + void compute_fixed_point_error_device(pdhg_solver_state_t *state); + + void compute_fixed_point_error_host(pdhg_solver_state_t *state); + + void compute_fixed_point_error(pdhg_solver_state_t *state); + + void check_termination_criteria(pdhg_solver_state_t *solver_state, const termination_criteria_t *criteria); + +#ifdef __cplusplus +} +#endif diff --git a/internal/utils.h b/internal/utils.h index dc04747..77e5837 100644 --- a/internal/utils.h +++ b/internal/utils.h @@ -123,8 +123,6 @@ extern "C" const restart_parameters_t *restart_params, int termination_evaluation_frequency); - void check_termination_criteria(pdhg_solver_state_t *solver_state, const termination_criteria_t *criteria); - void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem); void filter_constraint_matrix_entries(lp_problem_t *out, const lp_problem_t *in, const pdhg_parameters_t *params); @@ -143,10 +141,6 @@ extern "C" int get_print_frequency(int iter); - void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm); - - void compute_infeasibility_information(pdhg_solver_state_t *state); - void fill_or_copy(double **dest, int n, const double *src, double fill_value); int dense_to_csr(const matrix_desc_t *desc, int **row_ptr, int **col_ind, double **vals, int *nnz_out); @@ -155,27 +149,6 @@ extern "C" int coo_to_csr(const matrix_desc_t *desc, int **row_ptr, int **col_ind, double **vals); - void check_feas_polishing_termination_criteria(pdhg_solver_state_t *solver_state, - const pdhg_solver_state_t *ori_solver_state, - const termination_criteria_t *criteria, - bool is_primal_polish); - - void print_initial_feas_polish_info(bool is_primal_polish, const pdhg_parameters_t *params); - - void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state, bool verbose, bool is_primal_polish); - - void pdhg_feas_polish_final_log(const pdhg_solver_state_t *primal_state, - const pdhg_solver_state_t *dual_state, - bool verbose); - - void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, - const pdhg_solver_state_t *ori_state, - norm_type_t optimality_norm); - - void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, - const pdhg_solver_state_t *ori_state, - norm_type_t optimality_norm); - void set_default_parameters(pdhg_parameters_t *params); #ifdef __cplusplus diff --git a/python/README.md b/python/README.md index 0a84ed1..c7debc5 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. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-14` | Relative tolerance on the ray certificate for declaring primal/dual infeasibility. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix), and its objective term has the right sign with a fixed relative margin. `0` disables infeasibility detection. | +| `InfeasibilityDetection` | `infeasibility_detection` | bool | `True` | Whether to look for infeasibility certificates at all. Setting `InfeasibleTol` to `0` disables detection too. | | `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. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index fb0cb2f..412426f 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -42,6 +42,7 @@ "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", "InfeasibleTol": "eps_infeasible_relative", + "InfeasibilityDetection": "infeasibility_detection", # scaling / step size "GeoMeanIters": "geometric_mean_iterations", "RuizIters": "l_inf_ruiz_iterations", diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index dcb7f6c..a28b857 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -37,6 +37,7 @@ "bound_objective_rescaling", "feasibility_polishing", "presolve", + "infeasibility_detection", } ) _INT_PARAMS = frozenset( @@ -71,13 +72,13 @@ { "eps_optimal_relative", "eps_feasible_relative", - "eps_infeasible_relative", "eps_feas_polish_relative", "sv_tol", "infinite_bound", } ) -_NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol"}) +# eps_infeasible_relative = 0 also disables detection, like infeasibility_detection=False +_NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol", "eps_infeasible_relative"}) _STRING_PARAMS = frozenset({"optimality_norm"}) # int params are stored as C int32 on the backend diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index ffc6816..35f001f 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -297,6 +297,7 @@ static py::dict get_default_params_py() // feasiblity polishing d["feasibility_polishing"] = p.feasibility_polishing; + d["infeasibility_detection"] = p.infeasibility_detection; d["eps_feas_polish_relative"] = p.termination_criteria.eps_feas_polish_relative; // Termination criteria norm @@ -409,6 +410,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p // Feasibility Polishing getb("feasibility_polishing", p->feasibility_polishing); + getb("infeasibility_detection", p->infeasibility_detection); getf("eps_feas_polish_relative", p->termination_criteria.eps_feas_polish_relative); // Termination criteria norm @@ -437,8 +439,9 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p throw std::invalid_argument("eps_optimal_relative must be positive."); if (p->termination_criteria.eps_feasible_relative <= 0.0) throw std::invalid_argument("eps_feasible_relative must be positive."); - if (p->termination_criteria.eps_infeasible_relative <= 0.0) - throw std::invalid_argument("eps_infeasible_relative must be positive."); + if (p->termination_criteria.eps_infeasible_relative < 0.0) + throw std::invalid_argument( + "eps_infeasible_relative must be non-negative (0 disables infeasibility detection)."); if (p->termination_criteria.eps_feas_polish_relative <= 0.0) throw std::invalid_argument("eps_feas_polish_relative must be positive."); if (p->sv_tol <= 0.0) diff --git a/src/cli.c b/src/cli.c index 999c812..2ae78b5 100644 --- a/src/cli.c +++ b/src/cli.c @@ -123,6 +123,10 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir, // } } fprintf(outfile, "Precondition time (sec): %e\n", result->rescaling_time_sec); + fprintf(outfile, "Max Primal Ray Infeasibility: %e\n", result->max_primal_ray_infeasibility); + fprintf(outfile, "Primal Ray Linear Objective: %e\n", result->primal_ray_linear_objective); + fprintf(outfile, "Max Dual Ray Infeasibility: %e\n", result->max_dual_ray_infeasibility); + fprintf(outfile, "Dual Ray Objective: %e\n", result->dual_ray_objective); fprintf(outfile, "Runtime (sec): %e\n", result->cumulative_time_sec); fprintf(outfile, "Iterations Count: %d\n", result->total_count); fprintf(outfile, "Primal Objective Value: %e\n", result->primal_objective_value); @@ -206,9 +210,17 @@ void print_usage(const char *prog_name) fprintf(stderr, " --eps_feas_polish Relative feasibility " "polish tolerance (default: 1e-6).\n"); + fprintf(stderr, + " --eps_infeasible Infeasibility detection tolerance: " + "largest relative\n" + " constraint residual of an accepted ray " + "certificate (default: 1e-8; 0 also disables detection).\n"); fprintf(stderr, " --opt_norm " "Norm for optimality criteria: l2 or linf (default: l2).\n"); + fprintf(stderr, + " --no_infeasibility_detection " + "Disable infeasibility detection.\n"); fprintf(stderr, " --no_presolve " "Disable presolve (default: enabled).\n"); @@ -233,6 +245,7 @@ int main(int argc, char *argv[]) {"eps_opt", required_argument, 0, 1003}, {"eps_feas", required_argument, 0, 1004}, {"eps_feas_polish", required_argument, 0, 1006}, + {"eps_infeasible", required_argument, 0, 1019}, {"feasibility_polishing", no_argument, 0, 'f'}, {"geo_mean_iter", required_argument, 0, 1018}, {"l_inf_ruiz_iter", required_argument, 0, 1007}, @@ -243,6 +256,7 @@ int main(int argc, char *argv[]) {"sv_tol", required_argument, 0, 1012}, {"eval_freq", required_argument, 0, 1013}, {"opt_norm", required_argument, 0, 1014}, + {"no_infeasibility_detection", no_argument, 0, 1020}, {"no_presolve", no_argument, 0, 1015}, {"matrix_zero_tol", required_argument, 0, 1016}, {"infinite_bound", required_argument, 0, 1017}, @@ -277,6 +291,9 @@ int main(int argc, char *argv[]) case 1006: // --eps_feas_polish_relative params.termination_criteria.eps_feas_polish_relative = atof(optarg); break; + case 1019: // --eps_infeasible + params.termination_criteria.eps_infeasible_relative = atof(optarg); + break; case 'f': // --feasibility_polishing params.feasibility_polishing = true; break; @@ -322,6 +339,9 @@ int main(int argc, char *argv[]) } } break; + case 1020: // --no_infeasibility_detection + params.infeasibility_detection = false; + break; case 1015: // --no_presolve params.presolve = false; break; diff --git a/src/feasibility_polish.cu b/src/feasibility_polish.cu index 0584752..a023d87 100644 --- a/src/feasibility_polish.cu +++ b/src/feasibility_polish.cu @@ -15,6 +15,8 @@ limitations under the License. */ #include "feasibility_polish.h" +#include "reduce.h" +#include "termination.h" #include "utils.h" #include @@ -32,8 +34,13 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, int k_offset, double static void perform_primal_restart(pdhg_solver_state_t *state); static void perform_dual_restart(pdhg_solver_state_t *state); -static void compute_primal_fixed_point_error(pdhg_solver_state_t *state); -static void compute_dual_fixed_point_error(pdhg_solver_state_t *state); +static void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish); +static void compute_polish_fixed_point_error_host(pdhg_solver_state_t *state, bool is_primal_polish); +static void compute_polish_fixed_point_error(pdhg_solver_state_t *state, bool is_primal_polish); +static void +compute_polish_residual_device(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, bool is_primal_polish); +static void +compute_polish_residual_host(pdhg_solver_state_t *state, norm_type_t optimality_norm, bool is_primal_polish); static pdhg_solver_state_t *initialize_primal_feas_polish_state(const pdhg_solver_state_t *original_state); static pdhg_solver_state_t *initialize_dual_feas_polish_state(const pdhg_solver_state_t *original_state); void primal_feasibility_polish(const pdhg_parameters_t *params, @@ -45,14 +52,193 @@ void dual_feasibility_polish(const pdhg_parameters_t *params, void primal_feas_polish_state_free(pdhg_solver_state_t *state); void dual_feas_polish_state_free(pdhg_solver_state_t *state); __global__ void zero_finite_value_vectors_kernel(double *__restrict__ vec, int n); -__global__ void compute_delta_primal_solution_kernel(const double *__restrict__ initial_primal, - const double *__restrict__ pdhg_primal, - double *__restrict__ delta_primal, - int n_vars); -__global__ void compute_delta_dual_solution_kernel(const double *__restrict__ initial_dual, - const double *__restrict__ pdhg_dual, - double *__restrict__ delta_dual, - int n_cons); + +// Termination, logging and residuals of the polishing phases. + +static void check_feas_polishing_termination_criteria(pdhg_solver_state_t *solver_state, + const pdhg_solver_state_t *ori_solver_state, + const termination_criteria_t *criteria, + bool is_primal_polish) +{ + solver_state->cumulative_time_sec = (double)(clock() - solver_state->start_time) / CLOCKS_PER_SEC; + if (is_primal_polish) + { + if (solver_state->relative_primal_residual <= criteria->eps_feas_polish_relative) + { + solver_state->termination_reason = TERMINATION_REASON_FEAS_POLISH_SUCCESS; + return; + } + } + else + { + if (solver_state->relative_dual_residual <= criteria->eps_feas_polish_relative) + { + solver_state->termination_reason = TERMINATION_REASON_FEAS_POLISH_SUCCESS; + return; + } + } + if (solver_state->total_count >= criteria->iteration_limit) + { + solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; + return; + } + double total_time_sec = (double)(clock() - ori_solver_state->start_time) / CLOCKS_PER_SEC; + if (total_time_sec >= criteria->time_sec_limit) + { + solver_state->termination_reason = TERMINATION_REASON_TIME_LIMIT; + return; + } +} + +static void print_initial_feas_polish_info(bool is_primal_polish, const pdhg_parameters_t *params) +{ + if (!params->verbose) + { + return; + } + printf("---------------------------------------------------------------------------------------\n"); + printf("Starting %s Feasibility Polishing Phase with relative tolerance %.2e\n", + is_primal_polish ? "Primal" : "Dual", + params->termination_criteria.eps_feas_polish_relative); + printf("---------------------------------------------------------------------------------------\n"); + if (is_primal_polish) + printf("%s %s | %s | %s | %s \n", " iter", " time ", "pr obj", " abs pr res ", " rel pr res "); + else + printf("%s %s | %s | %s | %s \n", " iter", " time ", "du obj", " abs du res ", " rel du res "); + printf("---------------------------------------------------------------------------------------\n"); +} + +static void +pdhg_feas_polish_final_log(const pdhg_solver_state_t *primal_state, const pdhg_solver_state_t *dual_state, bool verbose) +{ + if (!verbose) + { + return; + } + printf("---------------------------------------------------------------------------------------\n"); + printf("Feasibility Polishing Summary\n"); + printf(" Primal Status : %s\n", termination_reason_to_string(primal_state->termination_reason)); + printf(" Primal Iterations : %d\n", primal_state->total_count); + printf(" Primal Time Usage : %.3g sec\n", primal_state->cumulative_time_sec); + printf(" Dual Status : %s\n", termination_reason_to_string(dual_state->termination_reason)); + printf(" Dual Iterations : %d\n", dual_state->total_count); + printf(" Dual Time Usage : %.3g sec\n", dual_state->cumulative_time_sec); + printf(" Primal Residual : %.3e\n", primal_state->relative_primal_residual); + printf(" Dual Residual : %.3e\n", dual_state->relative_dual_residual); + printf(" Primal Dual Gap : %.3e\n", + fabs(primal_state->primal_objective_value - dual_state->dual_objective_value) / + (1.0 + fabs(primal_state->primal_objective_value) + fabs(dual_state->dual_objective_value))); +} + +static void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state, bool verbose, bool is_primal_polish) +{ + if (!verbose) + { + return; + } + if (state->total_count % get_print_frequency(state->total_count) == 0) + { + if (is_primal_polish) + { + printf("%6d %.1e | %8.1e | %.1e | %.1e \n", + state->total_count, + state->cumulative_time_sec, + state->original_objective_sign * state->primal_objective_value, + state->absolute_primal_residual, + state->relative_primal_residual); + } + else + { + printf("%6d %.1e | %8.1e | %.1e | %.1e \n", + state->total_count, + state->cumulative_time_sec, + state->original_objective_sign * state->dual_objective_value, + state->absolute_dual_residual, + state->relative_dual_residual); + } + } +} + +// Residuals and objective values of a polishing iterate. The same reductions +// as the main termination check (termination.h) with one synchronization per +// check; the polish only needs its own side, so the other side's slots hold +// junk that is never read. The polishing states carry a modified problem, so +// the vectors the objective values need come from the original state: the +// primal polish zeroes its objective, the dual polish zeroes its primal +// iterate and its finite bounds. +static void +compute_polish_residual_device(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, bool is_primal_polish) +{ + if (is_primal_polish) + { + cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->pdhg_primal_solution, state->primal_product); + launch_residual_check(state, + ori_state->objective_vector, + state->pdhg_primal_solution, + state->constraint_lower_bound_finite_val, + state->constraint_upper_bound_finite_val); + } + else + { + cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->pdhg_dual_solution, state->dual_product); + launch_residual_check(state, + state->objective_vector, + ori_state->pdhg_primal_solution, + ori_state->constraint_lower_bound_finite_val, + ori_state->constraint_upper_bound_finite_val); + } +} + +static void compute_polish_residual_host(pdhg_solver_state_t *state, norm_type_t optimality_norm, bool is_primal_polish) +{ + const device_scalars_t *v = state->scalars_h; + bool linf = optimality_norm == NORM_TYPE_L_INF; + double scale = state->constraint_bound_rescaling * state->objective_vector_rescaling; + if (is_primal_polish) + { + state->absolute_primal_residual = + (linf ? v->primal_res_max : sqrt(v->primal_res_sq)) / state->constraint_bound_rescaling; + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); + state->primal_objective_value = v->objective_dot / scale + state->objective_constant; + } + else + { + state->absolute_dual_residual = + (linf ? v->dual_res_max : sqrt(v->dual_res_sq)) / state->objective_vector_rescaling; + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); + state->dual_objective_value = (v->dual_slack_dot + v->primal_slack_sum) / scale + state->objective_constant; + } +} + +// The fixed point error of a polishing iterate: only its own side of the +// movement term, without the square root (the restart test compares ratios). +static void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish) +{ + compute_solution_difference(state, + state->pdhg_primal_solution, + state->reflected_primal_solution, + state->pdhg_dual_solution, + state->reflected_dual_solution); + const double *delta = is_primal_polish ? state->delta_primal_solution : state->delta_dual_solution; + int len = is_primal_polish ? state->num_variables : state->num_constraints; + double *out = is_primal_polish ? &state->scalars_d->fixed_point_primal_sq : &state->scalars_d->fixed_point_dual_sq; + reduce_sum(state, len, [=] __host__ __device__(int i) { return delta[i] * delta[i]; }, out); +} + +static void compute_polish_fixed_point_error_host(pdhg_solver_state_t *state, bool is_primal_polish) +{ + const device_scalars_t *v = state->scalars_h; + state->fixed_point_error = is_primal_polish ? v->fixed_point_primal_sq * state->primal_weight + : v->fixed_point_dual_sq / state->primal_weight; +} + +static void compute_polish_fixed_point_error(pdhg_solver_state_t *state, bool is_primal_polish) +{ + compute_polish_fixed_point_error_device(state, is_primal_polish); + fetch_device_scalars(state); + compute_polish_fixed_point_error_host(state, is_primal_polish); +} + // Feasibility Polishing void feasibility_polish(const pdhg_parameters_t *params, pdhg_solver_state_t *state) { @@ -139,7 +325,7 @@ void primal_feasibility_polish(const pdhg_parameters_t *params, if (do_restart) { - compute_primal_fixed_point_error(state); + compute_polish_fixed_point_error(state, true); state->initial_fixed_point_error = state->fixed_point_error; do_restart = false; } @@ -169,8 +355,11 @@ void primal_feasibility_polish(const pdhg_parameters_t *params, } CUDA_CHECK(cudaGraphLaunch(graphExec, state->stream)); - compute_primal_fixed_point_error(state); - compute_primal_feas_polish_residual(state, ori_state, params->optimality_norm); + compute_polish_fixed_point_error_device(state, true); + compute_polish_residual_device(state, ori_state, true); + fetch_device_scalars(state); + compute_polish_fixed_point_error_host(state, true); + compute_polish_residual_host(state, params->optimality_norm, true); state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; @@ -214,7 +403,7 @@ void dual_feasibility_polish(const pdhg_parameters_t *params, if (do_restart) { - compute_dual_fixed_point_error(state); + compute_polish_fixed_point_error(state, false); state->initial_fixed_point_error = state->fixed_point_error; do_restart = false; } @@ -244,8 +433,11 @@ void dual_feasibility_polish(const pdhg_parameters_t *params, } CUDA_CHECK(cudaGraphLaunch(graphExec, state->stream)); - compute_dual_fixed_point_error(state); - compute_dual_feas_polish_residual(state, ori_state, params->optimality_norm); + compute_polish_fixed_point_error_device(state, false); + compute_polish_residual_device(state, ori_state, false); + fetch_device_scalars(state); + compute_polish_fixed_point_error_host(state, false); + compute_polish_residual_host(state, params->optimality_norm, false); state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; @@ -561,50 +753,3 @@ static void perform_dual_restart(pdhg_solver_state_t *state) state->inner_count = 0; state->last_trial_fixed_point_error = INFINITY; } - -__global__ void compute_delta_primal_solution_kernel(const double *__restrict__ initial_primal, - const double *__restrict__ pdhg_primal, - double *__restrict__ delta_primal, - int n_vars) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < n_vars) - { - delta_primal[i] = pdhg_primal[i] - initial_primal[i]; - } -} - -__global__ void compute_delta_dual_solution_kernel(const double *__restrict__ initial_dual, - const double *__restrict__ pdhg_dual, - double *__restrict__ delta_dual, - int n_cons) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < n_cons) - { - delta_dual[i] = pdhg_dual[i] - initial_dual[i]; - } -} - -static void compute_primal_fixed_point_error(pdhg_solver_state_t *state) -{ - compute_delta_primal_solution_kernel<<num_blocks_primal, THREADS_PER_BLOCK, 0, state->stream>>>( - state->pdhg_primal_solution, - state->reflected_primal_solution, - state->delta_primal_solution, - state->num_variables); - double primal_norm = 0.0; - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->delta_primal_solution, 1, &primal_norm)); - state->fixed_point_error = primal_norm * primal_norm * state->primal_weight; -} - -static void compute_dual_fixed_point_error(pdhg_solver_state_t *state) -{ - compute_delta_dual_solution_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->pdhg_dual_solution, state->reflected_dual_solution, state->delta_dual_solution, state->num_constraints); - double dual_norm = 0.0; - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->delta_dual_solution, 1, &dual_norm)); - state->fixed_point_error = dual_norm * dual_norm / state->primal_weight; -} diff --git a/src/infeasibility.cu b/src/infeasibility.cu new file mode 100644 index 0000000..a3f1e71 --- /dev/null +++ b/src/infeasibility.cu @@ -0,0 +1,589 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +// Detection of primal and dual infeasibility from ray certificates. See the +// comment above compute_infeasibility_device for the method, and +// docs/infeasibility_detection/report.tex for its derivation and evaluation. + +#include "infeasibility.h" +#include "reduce.h" +#include "termination.h" +#include "utils.h" + +#include + +static const double HOST_ONE = 1.0; +static const double HOST_ZERO = 0.0; + +// Statistics of one candidate ray, in the original problem space. +typedef struct +{ + double max_primal_ray_infeasibility; + double primal_ray_linear_objective; + double max_dual_ray_infeasibility; + double dual_ray_objective; + // the objective terms with absolute values, |c|^T|d| and the absolute-value + // sum of the dual objective terms, in the same scale as the objective terms + // above; the objective term relative to this is the margin t3 + double primal_objective_abs; + double dual_objective_abs; +} ray_candidate_stats_t; + +__global__ void abs_kernel(double *__restrict__ out, const double *__restrict__ x, int n) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n) + out[i] = fabs(x[i]); +} + +// --------------------------------------------------------------------------- +// Batched infeasibility detection. The three candidate directions are stacked +// column by column into dense n x K and m x K matrices and multiplied by A and +// A^T with one SpMM each, so that the matrix is read once per product instead +// of once per candidate; the products |A||d| and |A^T||y| use the absolute +// values of the matrix entries stored on the same CSR indices. Every scalar of +// a check comes out of the fused passes below (one over the variables and one +// over the constraints per candidate, all quantities at once, deterministic +// two-stage reductions into the solver state's scalar array), so a check +// synchronizes twice at most: once before the objective margin is tested and +// once after the absolute-value products. +// +// Memory: a check keeps four dense K-column matrices, the candidates D and Y +// and the products prod_m and prod_n. The first round fills prod_m with A D +// and prod_n with A^T Y and reduces everything that depends on them, including +// the violations that the second round's ratios need. The second round then no longer +// needs D, Y or the products, so it overwrites D with |D| in place and prod_m +// with |A||D| (likewise on the dual side) and reduces only their maxima. The +// next check rewrites all four matrices. +// --------------------------------------------------------------------------- + +// The scalars of a check are the infeas[] members of device_scalars_t +// (internal_types.h), one infeas_candidate_scalars_t per candidate. + +struct infeas_ctx +{ + double *abs_val_A; + double *abs_val_AT; + cusparseSpMatDescr_t matA, matAT, matAbsA, matAbsAT; + // column-major, one column per candidate; see "Memory" above for what + // each holds between the two rounds of a check + double *D; // n x K: primal candidates, then |D| + double *Y; // m x K: dual candidates, then |Y| + double *prod_m; // m x K: A D, then |A||D| + double *prod_n; // n x K: A^T Y, then |A^T||Y| + cusparseDnMatDescr_t dnD, dnY, dnProdM, dnProdN; + void *buf_A, *buf_AT, *buf_absA, *buf_absAT; + // carried from the margin test to the finish of a check + ray_candidate_stats_t stats[INFEAS_NUM_CANDIDATES]; + bool need_primal[INFEAS_NUM_CANDIDATES]; + bool need_dual[INFEAS_NUM_CANDIDATES]; +}; + +// The candidate directions, each x + beta * y of two device vectors (y may be +// NULL for x alone); passed to the kernels by value. +typedef struct +{ + const double *x[INFEAS_NUM_CANDIDATES]; + const double *y[INFEAS_NUM_CANDIDATES]; + double beta[INFEAS_NUM_CANDIDATES]; +} candidate_sources_t; + +__device__ __forceinline__ double candidate_entry(const candidate_sources_t &src, int c, int i) +{ + return src.y[c] ? src.x[c][i] + src.beta[c] * src.y[c][i] : src.x[c][i]; +} + +// Violation of the sign conditions of a primal ray: (A d)_i must be >= 0 when +// the constraint has a finite lower bound and <= 0 when it has a finite upper +// bound. +__host__ __device__ __forceinline__ double primal_ray_row_violation(double ad, bool lb_finite, bool ub_finite) +{ + return (lb_finite ? fmax(-ad, 0.0) : 0.0) + (ub_finite ? fmax(ad, 0.0) : 0.0); +} + +// Violation of the sign conditions of a dual ray on its reduced cost +// r_j = -(A^T y)_j: r_j > 0 needs a finite lower bound on x_j and r_j < 0 a +// finite upper bound. +__host__ __device__ __forceinline__ double dual_ray_reduced_cost_violation(double r, bool lb_finite, bool ub_finite) +{ + return (lb_finite ? 0.0 : fmax(r, 0.0)) + (ub_finite ? 0.0 : fmax(-r, 0.0)); +} + +// The two summands of the contribution of a dual ray entry v (y_i, or the +// reduced cost r_j) to the dual objective: v^+ times the finite lower bound +// and v^- times the finite upper bound (0 when the bound is infinite). +__host__ __device__ __forceinline__ double lower_bound_term(double v, double lower_bound_finite) +{ + return fmax(v, 0.0) * lower_bound_finite; +} + +__host__ __device__ __forceinline__ double upper_bound_term(double v, double upper_bound_finite) +{ + return fmin(v, 0.0) * upper_bound_finite; +} + +// Forms candidate blockIdx.y of a primal ray and projects it onto the sign +// conditions of a certificate. +__global__ void form_primal_candidates_kernel(candidate_sources_t src, + const double *__restrict__ var_lb, + const double *__restrict__ var_ub, + double *__restrict__ D, + int n) +{ + int j = blockIdx.x * blockDim.x + threadIdx.x; + int c = blockIdx.y; + if (j >= n) + return; + double v = candidate_entry(src, c, j); + if (isfinite(var_lb[j])) + v = fmax(v, 0.0); + if (isfinite(var_ub[j])) + v = fmin(v, 0.0); + D[(size_t)c * n + j] = v; +} + +__global__ void form_dual_candidates_kernel(candidate_sources_t src, + const double *__restrict__ con_lb, + const double *__restrict__ con_ub, + double *__restrict__ Y, + int m) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + int c = blockIdx.y; + if (i >= m) + return; + double v = candidate_entry(src, c, i); + if (!isfinite(con_lb[i])) + v = fmin(v, 0.0); + if (!isfinite(con_ub[i])) + v = fmax(v, 0.0); + Y[(size_t)c * m + i] = v; +} + +// Replaces every entry of column blockIdx.y of a column-major matrix with its +// absolute value. +__global__ void abs_columns_kernel(double *__restrict__ X, int rows) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < rows) + X[(size_t)blockIdx.y * rows + i] = fabs(X[(size_t)blockIdx.y * rows + i]); +} + +// First-round scalars of candidate c from its columns of D, Y and the products +// A D and A^T Y. Violations come in the scaled space (they decide the +// certificate) and in the original space (they are reported). +static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + const int n = state->num_variables, m = state->num_constraints; + infeas_candidate_scalars_t *out = &state->scalars_d->infeas[c]; + const double *d = ctx->D + (size_t)c * n; // the primal ray + const double *y = ctx->Y + (size_t)c * m; // the dual ray + const double *Ad = ctx->prod_m + (size_t)c * m; + const double *ATy = ctx->prod_n + (size_t)c * n; + const double *objective = state->objective_vector; + const double *variable_rescaling = state->variable_rescaling; + const double *constraint_rescaling = state->constraint_rescaling; + const double bound_rescaling = state->constraint_bound_rescaling; + const double objective_rescaling = state->objective_vector_rescaling; + const double *variable_lower = state->variable_lower_bound, *variable_upper = state->variable_upper_bound; + const double *variable_lower_finite = state->variable_lower_bound_finite_val; + const double *variable_upper_finite = state->variable_upper_bound_finite_val; + const double *constraint_lower = state->constraint_lower_bound, *constraint_upper = state->constraint_upper_bound; + const double *constraint_lower_finite = state->constraint_lower_bound_finite_val; + const double *constraint_upper_finite = state->constraint_upper_bound_finite_val; + + // the primal ray d: its norm in the original space, the objective term + // c^T d with and without absolute values, and the violation of the sign + // conditions on A d + auto row_violation = [=] __host__ __device__(int i) + { return primal_ray_row_violation(Ad[i], isfinite(constraint_lower[i]), isfinite(constraint_upper[i])); }; + reduce_max( + state, + n, + [=] __host__ __device__(int j) { return fabs(d[j]) / (variable_rescaling[j] * bound_rescaling); }, + &out->primal_ray.max_abs); + reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * d[j]; }, &out->primal_ray.obj); + reduce_sum( + state, n, [=] __host__ __device__(int j) { return fabs(objective[j] * d[j]); }, &out->primal_ray.obj_abs); + reduce_max(state, m, row_violation, &out->primal_ray.viol_max_scaled); + reduce_max( + state, + m, + [=] __host__ __device__(int i) { return row_violation(i) * constraint_rescaling[i]; }, + &out->primal_ray.viol_max_orig); + + // the dual ray y and its reduced cost r = -A^T y: the norm of y in the + // original space, the two parts of the dual objective term with and + // without absolute values, and the violation of the sign conditions on r + auto reduced_cost = [=] __host__ __device__(int j) { return -ATy[j]; }; + auto y_lower_term = [=] __host__ __device__(int i) { return lower_bound_term(y[i], constraint_lower_finite[i]); }; + auto y_upper_term = [=] __host__ __device__(int i) { return upper_bound_term(y[i], constraint_upper_finite[i]); }; + auto r_lower_term = [=] __host__ __device__(int j) + { return lower_bound_term(reduced_cost(j), variable_lower_finite[j]); }; + auto r_upper_term = [=] __host__ __device__(int j) + { return upper_bound_term(reduced_cost(j), variable_upper_finite[j]); }; + auto reduced_cost_violation = [=] __host__ __device__(int j) + { + return dual_ray_reduced_cost_violation( + reduced_cost(j), isfinite(variable_lower[j]), isfinite(variable_upper[j])); + }; + reduce_max( + state, + m, + [=] __host__ __device__(int i) { return fabs(y[i]) / (constraint_rescaling[i] * objective_rescaling); }, + &out->dual_ray.max_abs); + reduce_sum( + state, + m, + [=] __host__ __device__(int i) { return y_lower_term(i) + y_upper_term(i); }, + &out->dual_ray.obj_con_bounds); + reduce_sum( + state, + n, + [=] __host__ __device__(int j) { return r_lower_term(j) + r_upper_term(j); }, + &out->dual_ray.obj_var_bounds); + reduce_sum( + state, + m, + [=] __host__ __device__(int i) { return fabs(y_lower_term(i)) + fabs(y_upper_term(i)); }, + &out->dual_ray.obj_con_bounds_abs); + reduce_sum( + state, + n, + [=] __host__ __device__(int j) { return fabs(r_lower_term(j)) + fabs(r_upper_term(j)); }, + &out->dual_ray.obj_var_bounds_abs); + reduce_max(state, n, reduced_cost_violation, &out->dual_ray.viol_max_scaled); + reduce_max( + state, + n, + [=] __host__ __device__(int j) { return reduced_cost_violation(j) * variable_rescaling[j]; }, + &out->dual_ray.viol_max_orig); +} + +static void spmm_buffer( + cusparseHandle_t handle, cusparseSpMatDescr_t mat, cusparseDnMatDescr_t in, cusparseDnMatDescr_t out, void **buffer) +{ + size_t size = 0; + CUSPARSE_CHECK(cusparseSpMM_bufferSize(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + CUSPARSE_OPERATION_NON_TRANSPOSE, + &HOST_ONE, + mat, + in, + &HOST_ZERO, + out, + CUDA_R_64F, + CUSPARSE_SPMM_ALG_DEFAULT, + &size)); + CUDA_CHECK(cudaMalloc(buffer, size > 0 ? size : 16)); +} + +static void +spmm(cusparseHandle_t handle, cusparseSpMatDescr_t mat, cusparseDnMatDescr_t in, cusparseDnMatDescr_t out, void *buffer) +{ + CUSPARSE_CHECK(cusparseSpMM(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + CUSPARSE_OPERATION_NON_TRANSPOSE, + &HOST_ONE, + mat, + in, + &HOST_ZERO, + out, + CUDA_R_64F, + CUSPARSE_SPMM_ALG_DEFAULT, + buffer)); +} + +infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state) +{ + infeas_ctx_t *ctx = (infeas_ctx_t *)safe_calloc(1, sizeof(infeas_ctx_t)); + const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; + const cu_sparse_matrix_csr_t *A = state->constraint_matrix, *AT = state->constraint_matrix_t; + const int nnz = A->num_nonzeros; + + CUDA_CHECK(cudaMalloc(&ctx->abs_val_A, (size_t)nnz * sizeof(double))); + CUDA_CHECK(cudaMalloc(&ctx->abs_val_AT, (size_t)nnz * sizeof(double))); + abs_kernel<<num_blocks_nnz, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->abs_val_A, A->val, nnz); + abs_kernel<<num_blocks_nnz, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->abs_val_AT, AT->val, nnz); + +#define CREATE_CSR(descr, M, vals) \ + CUSPARSE_CHECK(cusparseCreateCsr(&(descr), \ + (M)->num_rows, \ + (M)->num_cols, \ + (M)->num_nonzeros, \ + (M)->row_ptr, \ + (M)->col_ind, \ + (vals), \ + CUSPARSE_INDEX_32I, \ + CUSPARSE_INDEX_32I, \ + CUSPARSE_INDEX_BASE_ZERO, \ + CUDA_R_64F)) + CREATE_CSR(ctx->matA, A, A->val); + CREATE_CSR(ctx->matAT, AT, AT->val); + CREATE_CSR(ctx->matAbsA, A, ctx->abs_val_A); + CREATE_CSR(ctx->matAbsAT, AT, ctx->abs_val_AT); +#undef CREATE_CSR + +#define CREATE_DN(ptr, descr, rows) \ + CUDA_CHECK(cudaMalloc(&(ptr), (size_t)(rows) * K * sizeof(double))); \ + CUSPARSE_CHECK(cusparseCreateDnMat(&(descr), (rows), K, (rows), (ptr), CUDA_R_64F, CUSPARSE_ORDER_COL)) + CREATE_DN(ctx->D, ctx->dnD, n); + CREATE_DN(ctx->prod_n, ctx->dnProdN, n); + CREATE_DN(ctx->Y, ctx->dnY, m); + CREATE_DN(ctx->prod_m, ctx->dnProdM, m); +#undef CREATE_DN + + spmm_buffer(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, &ctx->buf_A); + spmm_buffer(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, &ctx->buf_AT); + spmm_buffer(state->sparse_handle, ctx->matAbsA, ctx->dnD, ctx->dnProdM, &ctx->buf_absA); + spmm_buffer(state->sparse_handle, ctx->matAbsAT, ctx->dnY, ctx->dnProdN, &ctx->buf_absAT); + return ctx; +} + +void infeasibility_context_destroy(infeas_ctx_t *ctx) +{ + if (!ctx) + return; + cusparseSpMatDescr_t mats[] = {ctx->matA, ctx->matAT, ctx->matAbsA, ctx->matAbsAT}; + for (size_t i = 0; i < sizeof(mats) / sizeof(mats[0]); i++) + if (mats[i]) + CUSPARSE_CHECK(cusparseDestroySpMat(mats[i])); + cusparseDnMatDescr_t dns[] = {ctx->dnD, ctx->dnY, ctx->dnProdM, ctx->dnProdN}; + for (size_t i = 0; i < sizeof(dns) / sizeof(dns[0]); i++) + if (dns[i]) + CUSPARSE_CHECK(cusparseDestroyDnMat(dns[i])); + void *bufs[] = {ctx->abs_val_A, + ctx->abs_val_AT, + ctx->D, + ctx->Y, + ctx->prod_m, + ctx->prod_n, + ctx->buf_A, + ctx->buf_AT, + ctx->buf_absA, + ctx->buf_absAT}; + for (size_t i = 0; i < sizeof(bufs) / sizeof(bufs[0]); i++) + if (bufs[i]) + CUDA_CHECK(cudaFree(bufs[i])); + free(ctx); +} + +// Threshold on the objective-side margin of a certificate (t3 below). It only +// has to clear rounding noise and is deliberately independent of +// eps_infeasible, so that loosening the tolerance can only make detection +// easier. +#define INFEAS_OBJECTIVE_MARGIN 1e-8 + +// On an infeasible problem the PDHG iterates diverge along the infimal +// displacement vector v, whose primal or dual part is an infeasibility +// certificate (Applegate et al. 2021; Lu and Yang 2024 for Halpern PDHG). Three +// candidate directions are evaluated at every termination check: +// 0. difference of iterates T(z^k) - z^k, which converges to v at rate +// sqrt(log k)/k within a restart epoch and linearly at the epoch starts; +// 1. the iterate T(z^k) itself, the anchor plus the displacement below; +// 2. displacement since the last restart T(z^k) - z^0, the normalized +// Halpern iterate of Lu and Yang, rate 1/k within an epoch. +// The three are complementary in when they are close to v, and no pair matches +// all three on the benchmarks. +// +// Each candidate is projected onto the sign conditions of a certificate and +// normalized by its infinity norm in the original space; its primal part d is +// a certificate of dual infeasibility and its dual part y one of primal +// infeasibility when +// the objective term has the right sign with a relative margin, +// -c^T d >= INFEAS_OBJECTIVE_MARGIN |c|^T |d| (t3, primal side) +// Omega(y) >= INFEAS_OBJECTIVE_MARGIN sum |terms| (t3, dual side) +// and the violation of the homogeneous constraints, relative to the +// activity the ray produces, is at most eps_infeasible: +// viol(A d) / max_i (|A||d|)_i <= eps (t2, primal side) +// viol(A^T y) / max_j (|A^T||y|)_j <= eps (t2, dual side). +// Both conditions are relative to the data: a ray with t2 <= eps is an exact +// certificate for a problem whose matrix differs from A by a relative amount +// eps, and its objective condition survives relative perturbations of c below +// the margin. The plain PDLP ratio viol / (-c^T d) is not used: it is not scale +// invariant, and a large objective vector lets junk directions pass. The +// certificate error stored in the state is t2 when the margin holds and +// INFINITY otherwise; separately for the primal and the dual ray, the candidate +// with the smallest error is kept and its statistics are reported. +// +// First round, kernels only: the candidates, A D and A^T Y, and every scalar +// that depends on them. infeasibility_context_create must have run. +void compute_infeasibility_device(pdhg_solver_state_t *state) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; + + // 1. the candidates, one per column, projected + candidate_sources_t ps, ds; + // T(z) - z = reflected - pdhg, since reflected = 2 T(z) - z + ps.x[0] = state->reflected_primal_solution, ps.y[0] = state->pdhg_primal_solution, ps.beta[0] = -1.0; + ds.x[0] = state->reflected_dual_solution, ds.y[0] = state->pdhg_dual_solution, ds.beta[0] = -1.0; + // T(z) + ps.x[1] = state->pdhg_primal_solution, ps.y[1] = NULL, ps.beta[1] = 0.0; + ds.x[1] = state->pdhg_dual_solution, ds.y[1] = NULL, ds.beta[1] = 0.0; + // T(z) - z^0 + ps.x[2] = state->pdhg_primal_solution, ps.y[2] = state->initial_primal_solution, ps.beta[2] = -1.0; + ds.x[2] = state->pdhg_dual_solution, ds.y[2] = state->initial_dual_solution, ds.beta[2] = -1.0; + form_primal_candidates_kernel<<num_blocks_primal, K), THREADS_PER_BLOCK, 0, state->stream>>>( + ps, state->variable_lower_bound, state->variable_upper_bound, ctx->D, n); + form_dual_candidates_kernel<<num_blocks_dual, K), THREADS_PER_BLOCK, 0, state->stream>>>( + ds, state->constraint_lower_bound, state->constraint_upper_bound, ctx->Y, m); + + // 2. A D and A^T Y in one product each + spmm(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, ctx->buf_A); + spmm(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, ctx->buf_AT); + + // 3. every scalar that depends on them, including the violations the + // certificate ratios need, so that the second round may overwrite D, Y + // and the products + for (int c = 0; c < K; c++) + { + launch_candidate_reductions(state, c); + } +} + +// Statistics of every candidate in the original space and the objective-side +// margin test, from the fetched first-round scalars. Returns whether any +// candidate passed the test, in which case the second round is needed. +static bool evaluate_objective_margins(pdhg_solver_state_t *state) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + const double bc = state->constraint_bound_rescaling, oc = state->objective_vector_rescaling; + + bool any = false; + for (int c = 0; c < INFEAS_NUM_CANDIDATES; c++) + { + const infeas_candidate_scalars_t *v = &state->scalars_h->infeas[c]; + ray_candidate_stats_t *s = &ctx->stats[c]; + double primal_norm = v->primal_ray.max_abs, dual_norm = v->dual_ray.max_abs; + s->primal_ray_linear_objective = primal_norm > 0.0 ? v->primal_ray.obj / (bc * oc) / primal_norm : 0.0; + s->primal_objective_abs = primal_norm > 0.0 ? v->primal_ray.obj_abs / (bc * oc) / primal_norm : 0.0; + s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max_orig / bc / primal_norm : 0.0; + s->dual_ray_objective = + dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds + v->dual_ray.obj_var_bounds) / (bc * oc) / dual_norm : 0.0; + s->dual_objective_abs = dual_norm > 0.0 + ? (v->dual_ray.obj_con_bounds_abs + v->dual_ray.obj_var_bounds_abs) / (bc * oc) / dual_norm + : 0.0; + s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max_orig / oc / dual_norm : 0.0; + ctx->need_primal[c] = s->primal_ray_linear_objective < 0.0 && s->primal_objective_abs > 0.0 && + -s->primal_ray_linear_objective >= INFEAS_OBJECTIVE_MARGIN * s->primal_objective_abs; + ctx->need_dual[c] = s->dual_ray_objective > 0.0 && s->dual_objective_abs > 0.0 && + s->dual_ray_objective >= INFEAS_OBJECTIVE_MARGIN * s->dual_objective_abs; + any = any || ctx->need_primal[c] || ctx->need_dual[c]; + } + return any; +} + +// Second round, kernels only: |A||d| or |A^T||y| and their maxima for each side on +// which some candidate passed the margin (all candidates of that side at +// once). D becomes |D| in place and prod_m receives |A||D|; likewise Y and +// prod_n. +static void launch_abs_products(pdhg_solver_state_t *state) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; + + bool any_primal = false, any_dual = false; + for (int c = 0; c < K; c++) + { + any_primal = any_primal || ctx->need_primal[c]; + any_dual = any_dual || ctx->need_dual[c]; + } + if (any_primal) + { + abs_columns_kernel<<num_blocks_primal, K), THREADS_PER_BLOCK, 0, state->stream>>>(ctx->D, n); + spmm(state->sparse_handle, ctx->matAbsA, ctx->dnD, ctx->dnProdM, ctx->buf_absA); + for (int c = 0; c < K; c++) + { + const double *p = ctx->prod_m + (size_t)c * m; + reduce_max( + state, + m, + [=] __host__ __device__(int i) { return p[i]; }, + &state->scalars_d->infeas[c].primal_ray.product_abs_max); + } + } + if (any_dual) + { + abs_columns_kernel<<num_blocks_dual, K), THREADS_PER_BLOCK, 0, state->stream>>>(ctx->Y, m); + spmm(state->sparse_handle, ctx->matAbsAT, ctx->dnY, ctx->dnProdN, ctx->buf_absAT); + for (int c = 0; c < K; c++) + { + const double *p = ctx->prod_n + (size_t)c * n; + reduce_max( + state, + n, + [=] __host__ __device__(int j) { return p[j]; }, + &state->scalars_d->infeas[c].dual_ray.product_abs_max); + } + } +} + +// Certificate errors and the best candidate per side. When no candidate of a +// side passed the margin, that side reports the statistics of candidate 0 +// with an infinite error, so the termination test never accepts it. +static void select_best_candidates(pdhg_solver_state_t *state) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + int best_p = 0, best_d = 0; + double best_perr = INFINITY, best_derr = INFINITY; + for (int c = 0; c < INFEAS_NUM_CANDIDATES; c++) + { + const infeas_candidate_scalars_t *v = &state->scalars_h->infeas[c]; + double perr = INFINITY, derr = INFINITY; + if (ctx->need_primal[c]) + { + double viol = v->primal_ray.viol_max_scaled, act = v->primal_ray.product_abs_max; + perr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); + } + if (ctx->need_dual[c]) + { + double viol = v->dual_ray.viol_max_scaled, act = v->dual_ray.product_abs_max; + derr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); + } + if (perr < best_perr) + { + best_perr = perr; + best_p = c; + } + if (derr < best_derr) + { + best_derr = derr; + best_d = c; + } + } + state->max_primal_ray_infeasibility = ctx->stats[best_p].max_primal_ray_infeasibility; + state->primal_ray_linear_objective = ctx->stats[best_p].primal_ray_linear_objective; + state->primal_ray_certificate_error = best_perr; + state->max_dual_ray_infeasibility = ctx->stats[best_d].max_dual_ray_infeasibility; + state->dual_ray_objective = ctx->stats[best_d].dual_ray_objective; + state->dual_ray_certificate_error = best_derr; +} + +// Turns the fetched first-round scalars into the certificate errors and ray +// statistics of the state. Host arithmetic only, unless some candidate passes +// its objective margin: then the second round runs here, with one more copy +// and synchronization. +void compute_infeasibility_host(pdhg_solver_state_t *state) +{ + if (evaluate_objective_margins(state)) + { + launch_abs_products(state); + fetch_device_scalars(state); + } + select_best_candidates(state); +} diff --git a/src/reduce.cu b/src/reduce.cu new file mode 100644 index 0000000..7e472bd --- /dev/null +++ b/src/reduce.cu @@ -0,0 +1,51 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +// Temporary storage of the reductions; see reduce.h. + +#include "reduce.h" + +#include +#include + +void reduce_temp_storage_create(pdhg_solver_state_t *state) +{ + int n_max = state->num_variables > state->num_constraints ? state->num_variables : state->num_constraints; + if (n_max < 1) + n_max = 1; + const double *in = NULL; + double *out = NULL; + size_t sum_bytes = 0, max_bytes = 0; + CUDA_CHECK(cub::DeviceReduce::Sum(NULL, sum_bytes, in, out, n_max, state->stream)); + CUDA_CHECK(cub::DeviceReduce::Max(NULL, max_bytes, in, out, n_max, state->stream)); + // twice the queried size and a margin: the size depends only on the vector + // length and the value type, but reduce_check_temp_storage catches a + // transform iterator that needs more anyway + state->reduce_temp_bytes = 2 * (sum_bytes > max_bytes ? sum_bytes : max_bytes) + 256; + CUDA_CHECK(cudaMalloc(&state->reduce_temp_d, state->reduce_temp_bytes)); +} + +void reduce_check_temp_storage(const pdhg_solver_state_t *state, size_t bytes) +{ + if (bytes > state->reduce_temp_bytes) + { + fprintf(stderr, + "Error: a reduction needs %zu bytes of temporary storage, %zu allocated.\n", + bytes, + state->reduce_temp_bytes); + exit(EXIT_FAILURE); + } +} diff --git a/src/solver.cu b/src/solver.cu index 636c7b5..48bfca1 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -16,10 +16,13 @@ limitations under the License. #include "cupdlpx.h" #include "feasibility_polish.h" +#include "infeasibility.h" #include "internal_types.h" #include "preconditioner.h" #include "presolve.h" +#include "reduce.h" #include "solver.h" +#include "termination.h" #include "utils.h" #include @@ -48,21 +51,12 @@ __global__ void rescale_solution_kernel(double *__restrict__ primal_solution, const double constraint_bound_rescaling, int n_vars, int n_cons); -__global__ void compute_delta_solution_kernel(const double *__restrict__ initial_primal, - const double *__restrict__ pdhg_primal, - double *__restrict__ delta_primal, - const double *__restrict__ initial_dual, - const double *__restrict__ pdhg_dual, - double *__restrict__ delta_dual, - int n_vars, - int n_cons); static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, const lp_problem_t *original_problem); static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params); static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, const pdhg_parameters_t *params); static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_problem, const pdhg_parameters_t *params, objective_sense_t original_objective_sense); -static void compute_fixed_point_error(pdhg_solver_state_t *state); void pdhg_solver_state_free(pdhg_solver_state_t *state); void rescale_info_free(rescale_info_t *info); @@ -163,6 +157,13 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * cudaGraphExec_t graphExec = NULL; bool graph_created = false; + bool detect_infeasibility = + params->infeasibility_detection && params->termination_criteria.eps_infeasible_relative > 0.0; + if (detect_infeasibility) + { + state->infeas_ctx = infeasibility_context_create(state); + } + while (state->total_count < params->termination_criteria.iteration_limit) { sync_inner_count_to_gpu(state); @@ -191,6 +192,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * state, params->termination_evaluation_frequency, params->reflection_coefficient, true); compute_next_dual_solution( state, params->termination_evaluation_frequency, params->reflection_coefficient, true); + // end CUDA graph capture cudaGraph_t graph; @@ -200,10 +202,33 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * graph_created = true; } CUDA_CHECK(cudaGraphLaunch(graphExec, state->stream)); - compute_fixed_point_error(state); - compute_residual(state, params->optimality_norm); - compute_infeasibility_information(state); + // Termination check. The graph launch returns at once, so the host + // queues the kernels of the check (fixed point error, residuals and, if + // enabled, the infeasibility candidates) behind it while the GPU is + // still iterating; every scalar stays on the device until the one copy + // and synchronization of fetch_device_scalars. Only the infeasibility + // check may synchronize a second time, inside its host half, when some + // candidate passes its objective margin. + compute_fixed_point_error_device(state); + compute_residual_device(state); + if (detect_infeasibility) + { + compute_infeasibility_device(state); + } + fetch_device_scalars(state); + compute_fixed_point_error_host(state); + compute_residual_host(state, params->optimality_norm); + if (detect_infeasibility) + { + compute_infeasibility_host(state); + } + else + { + // detection off: no candidate evaluation, no |A| copies + state->primal_ray_certificate_error = INFINITY; + state->dual_ray_certificate_error = INFINITY; + } state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; @@ -454,6 +479,7 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ ALLOC_ZERO(state->dual_slack, var_bytes); ALLOC_ZERO(state->dual_residual, var_bytes); ALLOC_ZERO(state->delta_primal_solution, var_bytes); + ALLOC_ZERO(state->delta_dual_product, var_bytes); ALLOC_ZERO(state->initial_dual_solution, con_bytes); ALLOC_ZERO(state->current_dual_solution, con_bytes); @@ -603,26 +629,14 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ CUDA_CHECK(cudaMemset(state->d_dual_step_size, 0, sizeof(double))); CUDA_CHECK(cudaMemset(state->d_inner_count, 0, sizeof(int))); - CUDA_CHECK(cudaMalloc(&state->ones_primal_d, state->num_variables * sizeof(double))); - CUDA_CHECK(cudaMalloc(&state->ones_dual_d, state->num_constraints * sizeof(double))); - - double *ones_primal_h = (double *)safe_malloc(state->num_variables * sizeof(double)); - for (int i = 0; i < state->num_variables; ++i) - ones_primal_h[i] = 1.0; - CUDA_CHECK( - cudaMemcpy(state->ones_primal_d, ones_primal_h, state->num_variables * sizeof(double), cudaMemcpyHostToDevice)); - free(ones_primal_h); - - double *ones_dual_h = (double *)safe_malloc(state->num_constraints * sizeof(double)); - for (int i = 0; i < state->num_constraints; ++i) - ones_dual_h[i] = 1.0; - CUDA_CHECK( - cudaMemcpy(state->ones_dual_d, ones_dual_h, state->num_constraints * sizeof(double), cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMalloc(&state->scalars_d, sizeof(device_scalars_t))); + CUDA_CHECK(cudaHostAlloc(&state->scalars_h, sizeof(device_scalars_t), cudaHostAllocDefault)); // --- CUDA Graph Initialization --- CUDA_CHECK(cudaStreamCreate(&state->stream)); CUSPARSE_CHECK(cusparseSetStream(state->sparse_handle, state->stream)); CUBLAS_CHECK(cublasSetStream(state->blas_handle, state->stream)); + reduce_temp_storage_create(state); state->spmv_ctx = cupdlpx_spmv_ctx_create(state->sparse_handle, state->constraint_matrix, @@ -632,7 +646,6 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ state->pdhg_dual_solution, state->dual_product); - free(ones_dual_h); if (params->verbose) { printf("---------------------------------------------------------------------" @@ -856,27 +869,6 @@ __global__ void rescale_solution_kernel(double *__restrict__ primal_solution, } } -__global__ void compute_delta_solution_kernel(const double *__restrict__ initial_primal, - const double *__restrict__ pdhg_primal, - double *__restrict__ delta_primal, - const double *__restrict__ initial_dual, - const double *__restrict__ pdhg_dual, - double *__restrict__ delta_dual, - int n_vars, - int n_cons) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < n_vars) - { - delta_primal[i] = (pdhg_primal[i] - initial_primal[i]); - } - else if (i < n_vars + n_cons) - { - int idx = i - n_vars; - delta_dual[idx] = (pdhg_dual[idx] - initial_dual[idx]); - } -} - void compute_next_primal_solution(pdhg_solver_state_t *state, const int k_offset, const double reflection_coefficient, @@ -965,21 +957,26 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params) { - compute_delta_solution_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->initial_primal_solution, - state->pdhg_primal_solution, - state->delta_primal_solution, - state->initial_dual_solution, - state->pdhg_dual_solution, - state->delta_dual_solution, + compute_solution_difference(state, + state->initial_primal_solution, + state->pdhg_primal_solution, + state->initial_dual_solution, + state->pdhg_dual_solution); + + const double *delta_primal = state->delta_primal_solution, *delta_dual = state->delta_dual_solution; + reduce_sum( + state, state->num_variables, - state->num_constraints); - - double primal_dist, dual_dist; - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->delta_primal_solution, 1, &primal_dist)); - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->delta_dual_solution, 1, &dual_dist)); + [=] __host__ __device__(int j) { return delta_primal[j] * delta_primal[j]; }, + &state->scalars_d->restart_primal_sq); + reduce_sum( + state, + state->num_constraints, + [=] __host__ __device__(int i) { return delta_dual[i] * delta_dual[i]; }, + &state->scalars_d->restart_dual_sq); + fetch_device_scalars(state); + double primal_dist = sqrt(state->scalars_h->restart_primal_sq); + double dual_dist = sqrt(state->scalars_h->restart_dual_sq); double ratio_infeas = state->relative_dual_residual / state->relative_primal_residual; @@ -1009,22 +1006,26 @@ static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t state->best_primal_weight = state->primal_weight; } - CUDA_CHECK(cudaMemcpy(state->initial_primal_solution, - state->pdhg_primal_solution, - state->num_variables * sizeof(double), - cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(state->current_primal_solution, - state->pdhg_primal_solution, - state->num_variables * sizeof(double), - cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(state->initial_dual_solution, - state->pdhg_dual_solution, - state->num_constraints * sizeof(double), - cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(state->current_dual_solution, - state->pdhg_dual_solution, - state->num_constraints * sizeof(double), - cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpyAsync(state->initial_primal_solution, + state->pdhg_primal_solution, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaMemcpyAsync(state->current_primal_solution, + state->pdhg_primal_solution, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaMemcpyAsync(state->initial_dual_solution, + state->pdhg_dual_solution, + state->num_constraints * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); + CUDA_CHECK(cudaMemcpyAsync(state->current_dual_solution, + state->pdhg_dual_solution, + state->num_constraints * sizeof(double), + cudaMemcpyDeviceToDevice, + state->stream)); state->inner_count = 0; state->last_trial_fixed_point_error = INFINITY; @@ -1058,44 +1059,6 @@ static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, c state->best_primal_weight = state->primal_weight; } -static void compute_fixed_point_error(pdhg_solver_state_t *state) -{ - compute_delta_solution_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->pdhg_primal_solution, - state->reflected_primal_solution, - state->delta_primal_solution, - state->pdhg_dual_solution, - state->reflected_dual_solution, - state->delta_dual_solution, - state->num_variables, - state->num_constraints); - - cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->delta_dual_solution, state->dual_product); - - double interaction, movement; - - double primal_norm = 0.0; - double dual_norm = 0.0; - double cross_term = 0.0; - - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->delta_dual_solution, 1, &dual_norm)); - CUBLAS_CHECK( - cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->delta_primal_solution, 1, &primal_norm)); - movement = primal_norm * primal_norm * state->primal_weight + dual_norm * dual_norm / state->primal_weight; - - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->dual_product, - 1, - state->delta_primal_solution, - 1, - &cross_term)); - interaction = 2 * state->step_size * cross_term; - - state->fixed_point_error = sqrt(movement + interaction); -} - void pdhg_solver_state_free(pdhg_solver_state_t *state) { if (state == NULL) @@ -1105,6 +1068,8 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) if (state->spmv_ctx) cupdlpx_spmv_ctx_destroy(state->spmv_ctx); + if (state->infeas_ctx) + infeasibility_context_destroy(state->infeas_ctx); if (state->sparse_handle) CUSPARSE_CHECK(cusparseDestroy(state->sparse_handle)); if (state->blas_handle) @@ -1186,10 +1151,14 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) CUDA_CHECK(cudaFree(state->delta_primal_solution)); if (state->delta_dual_solution) CUDA_CHECK(cudaFree(state->delta_dual_solution)); - if (state->ones_primal_d) - CUDA_CHECK(cudaFree(state->ones_primal_d)); - if (state->ones_dual_d) - CUDA_CHECK(cudaFree(state->ones_dual_d)); + if (state->delta_dual_product) + CUDA_CHECK(cudaFree(state->delta_dual_product)); + if (state->scalars_d) + CUDA_CHECK(cudaFree(state->scalars_d)); + if (state->scalars_h) + CUDA_CHECK(cudaFreeHost(state->scalars_h)); + if (state->reduce_temp_d) + CUDA_CHECK(cudaFree(state->reduce_temp_d)); if (state->d_primal_step_size) CUDA_CHECK(cudaFree(state->d_primal_step_size)); if (state->d_dual_step_size) diff --git a/src/termination.cu b/src/termination.cu new file mode 100644 index 0000000..592179f --- /dev/null +++ b/src/termination.cu @@ -0,0 +1,285 @@ +/* +Copyright 2025 Haihao Lu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +// The termination check: fixed point error, residuals and objective values of +// the current iterate, and the criteria that decide whether the solver stops. +// See termination.h for the _device / _host protocol. + +#include "reduce.h" +#include "termination.h" +#include "utils.h" + +#include +#include + +__global__ void compute_delta_solution_kernel(const double *__restrict__ initial_primal, + const double *__restrict__ pdhg_primal, + double *__restrict__ delta_primal, + const double *__restrict__ initial_dual, + const double *__restrict__ pdhg_dual, + double *__restrict__ delta_dual, + int n_vars, + int n_cons) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n_vars) + { + delta_primal[i] = (pdhg_primal[i] - initial_primal[i]); + } + else if (i < n_vars + n_cons) + { + int idx = i - n_vars; + delta_dual[idx] = (pdhg_dual[idx] - initial_dual[idx]); + } +} + +void compute_solution_difference(pdhg_solver_state_t *state, + const double *from_primal, + const double *to_primal, + const double *from_dual, + const double *to_dual) +{ + compute_delta_solution_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( + from_primal, + to_primal, + state->delta_primal_solution, + from_dual, + to_dual, + state->delta_dual_solution, + state->num_variables, + state->num_constraints); +} + +__global__ void compute_residual_kernel(double *__restrict__ primal_residual, + const double *__restrict__ primal_product, + const double *__restrict__ constraint_lower_bound, + const double *__restrict__ constraint_upper_bound, + const double *__restrict__ dual_solution, + double *__restrict__ dual_residual, + const double *__restrict__ dual_product, + const double *__restrict__ dual_slack, + const double *__restrict__ objective_vector, + const double *__restrict__ constraint_rescaling, + const double *__restrict__ variable_rescaling, + double *__restrict__ primal_slack, + const double *__restrict__ const_lb_finite, + const double *__restrict__ const_ub_finite, + int num_constraints, + int num_variables) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + + if (i < num_constraints) + { + + double clamped_val = fmax(constraint_lower_bound[i], fmin(primal_product[i], constraint_upper_bound[i])); + primal_residual[i] = (primal_product[i] - clamped_val) * constraint_rescaling[i]; + + // the contribution of y_i to the dual objective, against the finite bounds + primal_slack[i] = + fmax(dual_solution[i], 0.0) * const_lb_finite[i] + fmin(dual_solution[i], 0.0) * const_ub_finite[i]; + } + else if (i < num_constraints + num_variables) + { + int idx = i - num_constraints; + dual_residual[idx] = (objective_vector[idx] - dual_product[idx] - dual_slack[idx]) * variable_rescaling[idx]; + } +} + +void fetch_device_scalars(pdhg_solver_state_t *state) +{ + CUDA_CHECK(cudaMemcpyAsync( + state->scalars_h, state->scalars_d, sizeof(device_scalars_t), cudaMemcpyDeviceToHost, state->stream)); + CUDA_CHECK(cudaStreamSynchronize(state->stream)); +} + +void launch_residual_check(pdhg_solver_state_t *state, + const double *objective_vector, + const double *primal_solution, + const double *constraint_lower_bound_finite, + const double *constraint_upper_bound_finite) +{ + compute_residual_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( + state->primal_residual, + state->primal_product, + state->constraint_lower_bound, + state->constraint_upper_bound, + state->pdhg_dual_solution, + state->dual_residual, + state->dual_product, + state->dual_slack, + objective_vector, + state->constraint_rescaling, + state->variable_rescaling, + state->primal_slack, + constraint_lower_bound_finite, + constraint_upper_bound_finite, + state->num_constraints, + state->num_variables); + + // the two norms of each residual (so that the norm type is a host-side + // choice), the primal objective c^T x and the two dual objective terms + device_scalars_t *out = state->scalars_d; + const int m = state->num_constraints, n = state->num_variables; + const double *primal_residual = state->primal_residual, *primal_slack = state->primal_slack; + const double *dual_residual = state->dual_residual, *dual_slack = state->dual_slack; + const double *objective = objective_vector, *x = primal_solution; + reduce_sum( + state, + m, + [=] __host__ __device__(int i) { return primal_residual[i] * primal_residual[i]; }, + &out->primal_res_sq); + reduce_max(state, m, [=] __host__ __device__(int i) { return fabs(primal_residual[i]); }, &out->primal_res_max); + reduce_sum(state, m, [=] __host__ __device__(int i) { return primal_slack[i]; }, &out->primal_slack_sum); + reduce_sum( + state, n, [=] __host__ __device__(int j) { return dual_residual[j] * dual_residual[j]; }, &out->dual_res_sq); + reduce_max(state, n, [=] __host__ __device__(int j) { return fabs(dual_residual[j]); }, &out->dual_res_max); + reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * x[j]; }, &out->objective_dot); + reduce_sum(state, n, [=] __host__ __device__(int j) { return dual_slack[j] * x[j]; }, &out->dual_slack_dot); +} + +// Queues the residuals and objective values of the current iterate; the +// results land in state->scalars_d. +void compute_residual_device(pdhg_solver_state_t *state) +{ + cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->pdhg_primal_solution, state->primal_product); + cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->pdhg_dual_solution, state->dual_product); + launch_residual_check(state, + state->objective_vector, + state->pdhg_primal_solution, + state->constraint_lower_bound_finite_val, + state->constraint_upper_bound_finite_val); +} + +// Turns the fetched scalars into the reported residuals and objective values. +void compute_residual_host(pdhg_solver_state_t *state, norm_type_t optimality_norm) +{ + const device_scalars_t *v = state->scalars_h; + bool linf = optimality_norm == NORM_TYPE_L_INF; + + state->absolute_primal_residual = + (linf ? v->primal_res_max : sqrt(v->primal_res_sq)) / state->constraint_bound_rescaling; + state->absolute_dual_residual = (linf ? v->dual_res_max : sqrt(v->dual_res_sq)) / state->objective_vector_rescaling; + + double scale = state->constraint_bound_rescaling * state->objective_vector_rescaling; + state->primal_objective_value = v->objective_dot / scale + state->objective_constant; + state->dual_objective_value = (v->dual_slack_dot + v->primal_slack_sum) / scale + state->objective_constant; + + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); + state->objective_gap = fabs(state->primal_objective_value - state->dual_objective_value); + state->relative_objective_gap = + state->objective_gap / (1.0 + fabs(state->primal_objective_value) + fabs(state->dual_objective_value)); +} + +void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm) +{ + compute_residual_device(state); + fetch_device_scalars(state); + compute_residual_host(state, optimality_norm); +} + +// Queues the fixed point error of the current iterate. +void compute_fixed_point_error_device(pdhg_solver_state_t *state) +{ + compute_solution_difference(state, + state->pdhg_primal_solution, + state->reflected_primal_solution, + state->pdhg_dual_solution, + state->reflected_dual_solution); + + cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->delta_dual_solution, state->delta_dual_product); + + // the squared norms of the two steps and the interaction term of the + // primal step with A^T times the dual step + device_scalars_t *out = state->scalars_d; + const int m = state->num_constraints, n = state->num_variables; + const double *delta_primal = state->delta_primal_solution, *delta_dual = state->delta_dual_solution; + const double *delta_dual_product = state->delta_dual_product; + reduce_sum( + state, m, [=] __host__ __device__(int i) { return delta_dual[i] * delta_dual[i]; }, &out->fixed_point_dual_sq); + reduce_sum( + state, + n, + [=] __host__ __device__(int j) { return delta_primal[j] * delta_primal[j]; }, + &out->fixed_point_primal_sq); + reduce_sum( + state, + n, + [=] __host__ __device__(int j) { return delta_dual_product[j] * delta_primal[j]; }, + &out->fixed_point_cross); +} + +void compute_fixed_point_error_host(pdhg_solver_state_t *state) +{ + const device_scalars_t *v = state->scalars_h; + double movement = v->fixed_point_primal_sq * state->primal_weight + v->fixed_point_dual_sq / state->primal_weight; + double interaction = 2.0 * state->step_size * v->fixed_point_cross; + state->fixed_point_error = sqrt(movement + interaction); +} + +void compute_fixed_point_error(pdhg_solver_state_t *state) +{ + compute_fixed_point_error_device(state); + fetch_device_scalars(state); + compute_fixed_point_error_host(state); +} + +static bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol) +{ + return state->relative_dual_residual < rel_feas_tol && state->relative_primal_residual < rel_feas_tol && + state->relative_objective_gap < rel_opt_tol; +} + +static bool primal_infeasibility_criteria_met(const pdhg_solver_state_t *state, double eps) +{ + return state->dual_ray_objective > 0.0 && state->dual_ray_certificate_error <= eps; +} + +static bool dual_infeasibility_criteria_met(const pdhg_solver_state_t *state, double eps) +{ + return state->primal_ray_linear_objective < 0.0 && state->primal_ray_certificate_error <= eps; +} + +void check_termination_criteria(pdhg_solver_state_t *solver_state, const termination_criteria_t *criteria) +{ + solver_state->cumulative_time_sec = (double)(clock() - solver_state->start_time) / CLOCKS_PER_SEC; + if (optimality_criteria_met(solver_state, criteria->eps_optimal_relative, criteria->eps_feasible_relative)) + { + solver_state->termination_reason = TERMINATION_REASON_OPTIMAL; + return; + } + if (primal_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) + { + solver_state->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE; + return; + } + if (dual_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) + { + solver_state->termination_reason = TERMINATION_REASON_DUAL_INFEASIBLE; + return; + } + if (solver_state->total_count >= criteria->iteration_limit) + { + solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; + return; + } + if (solver_state->cumulative_time_sec >= criteria->time_sec_limit) + { + solver_state->termination_reason = TERMINATION_REASON_TIME_LIMIT; + return; + } +} diff --git a/src/utils.cu b/src/utils.cu index 7ff69ce..059f8ba 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -18,6 +18,8 @@ limitations under the License. #include #include +#include + #ifndef CUPDLPX_VERSION #define CUPDLPX_VERSION "unknown" #endif @@ -222,60 +224,6 @@ const char *termination_reason_to_string(termination_reason_t reason) } } -bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol) -{ - return state->relative_dual_residual < rel_feas_tol && state->relative_primal_residual < rel_feas_tol && - state->relative_objective_gap < rel_opt_tol; -} - -bool primal_infeasibility_criteria_met(const pdhg_solver_state_t *state, double eps) -{ - if (state->dual_ray_objective <= 0.0) - { - return false; - } - return state->max_dual_ray_infeasibility / state->dual_ray_objective <= eps; -} - -bool dual_infeasibility_criteria_met(const pdhg_solver_state_t *state, double eps) -{ - if (state->primal_ray_linear_objective >= 0.0) - { - return false; - } - return state->max_primal_ray_infeasibility / (-state->primal_ray_linear_objective) <= eps; -} - -void check_termination_criteria(pdhg_solver_state_t *solver_state, const termination_criteria_t *criteria) -{ - solver_state->cumulative_time_sec = (double)(clock() - solver_state->start_time) / CLOCKS_PER_SEC; - if (optimality_criteria_met(solver_state, criteria->eps_optimal_relative, criteria->eps_feasible_relative)) - { - solver_state->termination_reason = TERMINATION_REASON_OPTIMAL; - return; - } - if (primal_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) - { - solver_state->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE; - return; - } - if (dual_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) - { - solver_state->termination_reason = TERMINATION_REASON_DUAL_INFEASIBLE; - return; - } - if (solver_state->total_count >= criteria->iteration_limit) - { - solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; - return; - } - if (solver_state->cumulative_time_sec >= criteria->time_sec_limit) - { - solver_state->termination_reason = TERMINATION_REASON_TIME_LIMIT; - return; - } -} - bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, const restart_parameters_t *restart_params, int termination_evaluation_frequency) @@ -328,7 +276,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->termination_criteria.time_sec_limit = 3600.0; params->termination_criteria.iteration_limit = INT32_MAX; params->termination_criteria.eps_feas_polish_relative = 1e-6; - params->termination_criteria.eps_infeasible_relative = 1e-14; + params->termination_criteria.eps_infeasible_relative = 1e-8; params->restart_params.artificial_restart_threshold = 0.36; params->restart_params.sufficient_reduction_for_restart = 0.2; @@ -340,6 +288,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->optimality_norm = NORM_TYPE_L2; params->presolve = true; + params->infeasibility_detection = true; params->matrix_zero_tol = 1e-9; params->infinite_bound = 1e20; } @@ -640,6 +589,7 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro params->termination_criteria.eps_infeasible_relative, default_params.termination_criteria.eps_infeasible_relative); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); + PRINT_DIFF_BOOL("infeasibility_detection", params->infeasibility_detection, default_params.infeasibility_detection); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); PRINT_DIFF_DBL("infinite_bound", params->infinite_bound, default_params.infinite_bound); } @@ -732,340 +682,6 @@ int get_print_frequency(int iter) return step; } -__global__ void compute_residual_kernel(double *__restrict__ primal_residual, - const double *__restrict__ primal_product, - const double *__restrict__ constraint_lower_bound, - const double *__restrict__ constraint_upper_bound, - const double *__restrict__ dual_solution, - double *__restrict__ dual_residual, - const double *__restrict__ dual_product, - const double *__restrict__ dual_slack, - const double *__restrict__ objective_vector, - const double *__restrict__ constraint_rescaling, - const double *__restrict__ variable_rescaling, - double *__restrict__ dual_obj_contribution, - const double *__restrict__ const_lb_finite, - const double *__restrict__ const_ub_finite, - int num_constraints, - int num_variables) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - - if (i < num_constraints) - { - - double clamped_val = fmax(constraint_lower_bound[i], fmin(primal_product[i], constraint_upper_bound[i])); - primal_residual[i] = (primal_product[i] - clamped_val) * constraint_rescaling[i]; - - dual_obj_contribution[i] = - fmax(dual_solution[i], 0.0) * const_lb_finite[i] + fmin(dual_solution[i], 0.0) * const_ub_finite[i]; - } - else if (i < num_constraints + num_variables) - { - int idx = i - num_constraints; - dual_residual[idx] = (objective_vector[idx] - dual_product[idx] - dual_slack[idx]) * variable_rescaling[idx]; - } -} - -__global__ void primal_infeasibility_project_kernel(double *__restrict__ primal_ray_estimate, - const double *__restrict__ variable_lower_bound, - const double *__restrict__ variable_upper_bound, - int num_variables) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < num_variables) - { - if (isfinite(variable_lower_bound[i])) - { - primal_ray_estimate[i] = fmax(primal_ray_estimate[i], 0.0); - } - if (isfinite(variable_upper_bound[i])) - { - primal_ray_estimate[i] = fmin(primal_ray_estimate[i], 0.0); - } - } -} - -__global__ void dual_infeasibility_project_kernel(double *__restrict__ dual_ray_estimate, - const double *__restrict__ constraint_lower_bound, - const double *__restrict__ constraint_upper_bound, - int num_constraints) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < num_constraints) - { - if (!isfinite(constraint_lower_bound[i])) - { - dual_ray_estimate[i] = fmin(dual_ray_estimate[i], 0.0); - } - if (!isfinite(constraint_upper_bound[i])) - { - dual_ray_estimate[i] = fmax(dual_ray_estimate[i], 0.0); - } - } -} - -__global__ void compute_primal_infeasibility_kernel(const double *__restrict__ primal_product, - const double *__restrict__ const_lb, - const double *__restrict__ const_ub, - int num_constraints, - double *__restrict__ primal_infeasibility, - const double *__restrict__ constraint_rescaling) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < num_constraints) - { - double pp_val = primal_product[i]; - primal_infeasibility[i] = - (fmax(0.0, -pp_val) * isfinite(const_lb[i]) + fmax(0.0, pp_val) * isfinite(const_ub[i])) * - constraint_rescaling[i]; - } -} - -__global__ void compute_dual_infeasibility_kernel(const double *dual_product, - const double *__restrict__ var_lb, - const double *__restrict__ var_ub, - int num_variables, - double *__restrict__ dual_infeasibility, - const double *__restrict__ variable_rescaling) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < num_variables) - { - double dp_val = -dual_product[i]; - dual_infeasibility[i] = (fmax(0.0, dp_val) * !isfinite(var_lb[i]) - fmin(0.0, dp_val) * !isfinite(var_ub[i])) * - variable_rescaling[i]; - } -} - -__global__ void -dual_solution_dual_objective_contribution_kernel(const double *__restrict__ constraint_lower_bound_finite_val, - const double *__restrict__ constraint_upper_bound_finite_val, - const double *__restrict__ dual_solution, - int num_constraints, - double *__restrict__ dual_objective_dual_solution_contribution_array) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - - if (i < num_constraints) - { - dual_objective_dual_solution_contribution_array[i] = - fmax(dual_solution[i], 0.0) * constraint_lower_bound_finite_val[i] + - fmin(dual_solution[i], 0.0) * constraint_upper_bound_finite_val[i]; - } -} - -__global__ void -dual_objective_dual_slack_contribution_array_kernel(const double *__restrict__ dual_slack, - double *__restrict__ dual_objective_dual_slack_contribution_array, - const double *__restrict__ variable_lower_bound_finite_val, - const double *__restrict__ variable_upper_bound_finite_val, - int num_variables) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - - if (i < num_variables) - { - dual_objective_dual_slack_contribution_array[i] = - fmax(-dual_slack[i], 0.0) * variable_lower_bound_finite_val[i] + - fmin(-dual_slack[i], 0.0) * variable_upper_bound_finite_val[i]; - } -} - -static double get_vector_inf_norm(cublasHandle_t handle, int n, const double *x_d) -{ - if (n <= 0) - return 0.0; - int index; - - cublasIdamax(handle, n, x_d, 1, &index); - double max_val; - - CUDA_CHECK(cudaMemcpy(&max_val, x_d + (index - 1), sizeof(double), cudaMemcpyDeviceToHost)); - return fabs(max_val); -} - -static double get_vector_sum(cublasHandle_t handle, int n, double *ones_d, const double *x_d) -{ - if (n <= 0) - return 0.0; - - double sum; - CUBLAS_CHECK(cublasDdot(handle, n, x_d, 1, ones_d, 1, &sum)); - return sum; -} - -void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm) -{ - cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->pdhg_primal_solution, state->primal_product); - cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->pdhg_dual_solution, state->dual_product); - - compute_residual_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->primal_residual, - state->primal_product, - state->constraint_lower_bound, - state->constraint_upper_bound, - state->pdhg_dual_solution, - state->dual_residual, - state->dual_product, - state->dual_slack, - state->objective_vector, - state->constraint_rescaling, - state->variable_rescaling, - state->primal_slack, - state->constraint_lower_bound_finite_val, - state->constraint_upper_bound_finite_val, - state->num_constraints, - state->num_variables); - - if (optimality_norm == NORM_TYPE_L_INF) - { - state->absolute_primal_residual = - get_vector_inf_norm(state->blas_handle, state->num_constraints, state->primal_residual); - } - else - { - CUBLAS_CHECK(cublasDnrm2_v2_64( - state->blas_handle, state->num_constraints, state->primal_residual, 1, &state->absolute_primal_residual)); - } - state->absolute_primal_residual /= state->constraint_bound_rescaling; - - if (optimality_norm == NORM_TYPE_L_INF) - { - state->absolute_dual_residual = - get_vector_inf_norm(state->blas_handle, state->num_variables, state->dual_residual); - } - else - { - CUBLAS_CHECK(cublasDnrm2_v2_64( - state->blas_handle, state->num_variables, state->dual_residual, 1, &state->absolute_dual_residual)); - } - - state->absolute_dual_residual /= state->objective_vector_rescaling; - - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->objective_vector, - 1, - state->pdhg_primal_solution, - 1, - &state->primal_objective_value)); - state->primal_objective_value = - state->primal_objective_value / (state->constraint_bound_rescaling * state->objective_vector_rescaling) + - state->objective_constant; - - double base_dual_objective; - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->dual_slack, - 1, - state->pdhg_primal_solution, - 1, - &base_dual_objective)); - double dual_slack_sum = - get_vector_sum(state->blas_handle, state->num_constraints, state->ones_dual_d, state->primal_slack); - state->dual_objective_value = (base_dual_objective + dual_slack_sum) / - (state->constraint_bound_rescaling * state->objective_vector_rescaling) + - state->objective_constant; - - state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); - - state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); - - state->objective_gap = fabs(state->primal_objective_value - state->dual_objective_value); - - state->relative_objective_gap = - state->objective_gap / (1.0 + fabs(state->primal_objective_value) + fabs(state->dual_objective_value)); -} - -void compute_infeasibility_information(pdhg_solver_state_t *state) -{ - primal_infeasibility_project_kernel<<num_blocks_primal, THREADS_PER_BLOCK, 0, state->stream>>>( - state->delta_primal_solution, state->variable_lower_bound, state->variable_upper_bound, state->num_variables); - dual_infeasibility_project_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->delta_dual_solution, - state->constraint_lower_bound, - state->constraint_upper_bound, - state->num_constraints); - - double primal_ray_inf_norm = - get_vector_inf_norm(state->blas_handle, state->num_variables, state->delta_primal_solution); - if (primal_ray_inf_norm > 0.0) - { - double scale = 1.0 / primal_ray_inf_norm; - cublasDscal(state->blas_handle, state->num_variables, &scale, state->delta_primal_solution, 1); - } - double dual_ray_inf_norm = - get_vector_inf_norm(state->blas_handle, state->num_constraints, state->delta_dual_solution); - - cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->delta_primal_solution, state->primal_product); - cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->delta_dual_solution, state->dual_product); - - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->objective_vector, - 1, - state->delta_primal_solution, - 1, - &state->primal_ray_linear_objective)); - state->primal_ray_linear_objective /= (state->constraint_bound_rescaling * state->objective_vector_rescaling); - - dual_solution_dual_objective_contribution_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->constraint_lower_bound_finite_val, - state->constraint_upper_bound_finite_val, - state->delta_dual_solution, - state->num_constraints, - state->primal_slack); - - dual_objective_dual_slack_contribution_array_kernel<<num_blocks_primal, - THREADS_PER_BLOCK, - 0, - state->stream>>>(state->dual_product, - state->dual_slack, - state->variable_lower_bound_finite_val, - state->variable_upper_bound_finite_val, - state->num_variables); - - double sum_primal_slack = - get_vector_sum(state->blas_handle, state->num_constraints, state->ones_dual_d, state->primal_slack); - double sum_dual_slack = - get_vector_sum(state->blas_handle, state->num_variables, state->ones_primal_d, state->dual_slack); - state->dual_ray_objective = - (sum_primal_slack + sum_dual_slack) / (state->constraint_bound_rescaling * state->objective_vector_rescaling); - - compute_primal_infeasibility_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->primal_product, - state->constraint_lower_bound, - state->constraint_upper_bound, - state->num_constraints, - state->primal_slack, - state->constraint_rescaling); - compute_dual_infeasibility_kernel<<num_blocks_primal, THREADS_PER_BLOCK, 0, state->stream>>>( - state->dual_product, - state->variable_lower_bound, - state->variable_upper_bound, - state->num_variables, - state->dual_slack, - state->variable_rescaling); - - state->max_primal_ray_infeasibility = - get_vector_inf_norm(state->blas_handle, state->num_constraints, state->primal_slack); - double dual_slack_norm = get_vector_inf_norm(state->blas_handle, state->num_variables, state->dual_slack); - state->max_dual_ray_infeasibility = dual_slack_norm; - - double scaling_factor = fmax(dual_ray_inf_norm, dual_slack_norm); - if (scaling_factor > 0.0) - { - state->max_dual_ray_infeasibility /= scaling_factor; - state->dual_ray_objective /= scaling_factor; - } - else - { - state->max_dual_ray_infeasibility = 0.0; - state->dual_ray_objective = 0.0; - } -} - // helper function to allocate and fill or copy an array void fill_or_copy(double **dst, int n, const double *src, double fill_val) { @@ -1228,241 +844,3 @@ int coo_to_csr(const matrix_desc_t *desc, int **row_ptr, int **col_ind, double * free(next); return 0; } - -void check_feas_polishing_termination_criteria(pdhg_solver_state_t *solver_state, - const pdhg_solver_state_t *ori_solver_state, - const termination_criteria_t *criteria, - bool is_primal_polish) -{ - solver_state->cumulative_time_sec = (double)(clock() - solver_state->start_time) / CLOCKS_PER_SEC; - if (is_primal_polish) - { - if (solver_state->relative_primal_residual <= criteria->eps_feas_polish_relative) - { - solver_state->termination_reason = TERMINATION_REASON_FEAS_POLISH_SUCCESS; - return; - } - } - else - { - if (solver_state->relative_dual_residual <= criteria->eps_feas_polish_relative) - { - solver_state->termination_reason = TERMINATION_REASON_FEAS_POLISH_SUCCESS; - return; - } - } - if (solver_state->total_count >= criteria->iteration_limit) - { - solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; - return; - } - double total_time_sec = (double)(clock() - ori_solver_state->start_time) / CLOCKS_PER_SEC; - if (total_time_sec >= criteria->time_sec_limit) - { - solver_state->termination_reason = TERMINATION_REASON_TIME_LIMIT; - return; - } -} - -void print_initial_feas_polish_info(bool is_primal_polish, const pdhg_parameters_t *params) -{ - if (!params->verbose) - { - return; - } - printf("---------------------------------------------------------------------------------------\n"); - printf("Starting %s Feasibility Polishing Phase with relative tolerance %.2e\n", - is_primal_polish ? "Primal" : "Dual", - params->termination_criteria.eps_feas_polish_relative); - printf("---------------------------------------------------------------------------------------\n"); - if (is_primal_polish) - printf("%s %s | %s | %s | %s \n", " iter", " time ", "pr obj", " abs pr res ", " rel pr res "); - else - printf("%s %s | %s | %s | %s \n", " iter", " time ", "du obj", " abs du res ", " rel du res "); - printf("---------------------------------------------------------------------------------------\n"); -} - -void pdhg_feas_polish_final_log(const pdhg_solver_state_t *primal_state, - const pdhg_solver_state_t *dual_state, - bool verbose) -{ - if (!verbose) - { - return; - } - printf("---------------------------------------------------------------------------------------\n"); - printf("Feasibility Polishing Summary\n"); - printf(" Primal Status : %s\n", termination_reason_to_string(primal_state->termination_reason)); - printf(" Primal Iterations : %d\n", primal_state->total_count); - printf(" Primal Time Usage : %.3g sec\n", primal_state->cumulative_time_sec); - printf(" Dual Status : %s\n", termination_reason_to_string(dual_state->termination_reason)); - printf(" Dual Iterations : %d\n", dual_state->total_count); - printf(" Dual Time Usage : %.3g sec\n", dual_state->cumulative_time_sec); - printf(" Primal Residual : %.3e\n", primal_state->relative_primal_residual); - printf(" Dual Residual : %.3e\n", dual_state->relative_dual_residual); - printf(" Primal Dual Gap : %.3e\n", - fabs(primal_state->primal_objective_value - dual_state->dual_objective_value) / - (1.0 + fabs(primal_state->primal_objective_value) + fabs(dual_state->dual_objective_value))); -} - -void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state, bool verbose, bool is_primal_polish) -{ - if (!verbose) - { - return; - } - if (state->total_count % get_print_frequency(state->total_count) == 0) - { - if (is_primal_polish) - { - printf("%6d %.1e | %8.1e | %.1e | %.1e \n", - state->total_count, - state->cumulative_time_sec, - state->original_objective_sign * state->primal_objective_value, - state->absolute_primal_residual, - state->relative_primal_residual); - } - else - { - printf("%6d %.1e | %8.1e | %.1e | %.1e \n", - state->total_count, - state->cumulative_time_sec, - state->original_objective_sign * state->dual_objective_value, - state->absolute_dual_residual, - state->relative_dual_residual); - } - } -} - -__global__ void compute_primal_feas_polish_residual_kernel(double *__restrict__ primal_residual, - const double *__restrict__ primal_product, - const double *__restrict__ constraint_lower_bound, - const double *__restrict__ constraint_upper_bound, - const double *__restrict__ constraint_rescaling, - int num_constraints) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - - if (i < num_constraints) - { - - double clamped_val = fmax(constraint_lower_bound[i], fmin(primal_product[i], constraint_upper_bound[i])); - primal_residual[i] = (primal_product[i] - clamped_val) * constraint_rescaling[i]; - } -} - -__global__ void compute_dual_feas_polish_residual_kernel(double *__restrict__ dual_residual, - const double *__restrict__ dual_solution, - const double *__restrict__ dual_product, - const double *__restrict__ dual_slack, - const double *__restrict__ objective_vector, - const double *__restrict__ variable_rescaling, - double *__restrict__ dual_obj_contribution, - const double *__restrict__ const_lb_finite, - const double *__restrict__ const_ub_finite, - int num_variables, - int num_constraints) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - - if (i < num_variables) - { - dual_residual[i] = (objective_vector[i] - dual_product[i] - dual_slack[i]) * variable_rescaling[i]; - } - else if (i < num_constraints + num_variables) - { - int idx = i - num_variables; - dual_obj_contribution[idx] = - fmax(dual_solution[idx], 0.0) * const_lb_finite[idx] + fmin(dual_solution[idx], 0.0) * const_ub_finite[idx]; - } -} - -void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, - const pdhg_solver_state_t *ori_state, - norm_type_t optimality_norm) -{ - cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->pdhg_primal_solution, state->primal_product); - - compute_primal_feas_polish_residual_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( - state->primal_residual, - state->primal_product, - state->constraint_lower_bound, - state->constraint_upper_bound, - state->constraint_rescaling, - state->num_constraints); - - if (optimality_norm == NORM_TYPE_L_INF) - { - state->absolute_primal_residual = - get_vector_inf_norm(state->blas_handle, state->num_constraints, state->primal_residual); - } - else - { - CUBLAS_CHECK(cublasDnrm2_v2_64( - state->blas_handle, state->num_constraints, state->primal_residual, 1, &state->absolute_primal_residual)); - } - - state->absolute_primal_residual /= state->constraint_bound_rescaling; - - state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); - - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - ori_state->objective_vector, - 1, - state->pdhg_primal_solution, - 1, - &state->primal_objective_value)); - state->primal_objective_value = - state->primal_objective_value / (state->constraint_bound_rescaling * state->objective_vector_rescaling) + - state->objective_constant; -} - -void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, - const pdhg_solver_state_t *ori_state, - norm_type_t optimality_norm) -{ - cupdlpx_spmv_ATx(state->sparse_handle, state->spmv_ctx, state->pdhg_dual_solution, state->dual_product); - - compute_dual_feas_polish_residual_kernel<<num_blocks_primal_dual, THREADS_PER_BLOCK>>>( - state->dual_residual, - state->pdhg_dual_solution, - state->dual_product, - state->dual_slack, - state->objective_vector, - state->variable_rescaling, - state->primal_slack, - ori_state->constraint_lower_bound_finite_val, - ori_state->constraint_upper_bound_finite_val, - state->num_variables, - state->num_constraints); - - if (optimality_norm == NORM_TYPE_L_INF) - { - state->absolute_dual_residual = - get_vector_inf_norm(state->blas_handle, state->num_variables, state->dual_residual); - } - else - { - CUBLAS_CHECK(cublasDnrm2_v2_64( - state->blas_handle, state->num_variables, state->dual_residual, 1, &state->absolute_dual_residual)); - } - - state->absolute_dual_residual /= state->objective_vector_rescaling; - - state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); - - double base_dual_objective; - CUBLAS_CHECK(cublasDdot(state->blas_handle, - state->num_variables, - state->dual_slack, - 1, - ori_state->pdhg_primal_solution, - 1, - &base_dual_objective)); - double dual_slack_sum = - get_vector_sum(state->blas_handle, state->num_constraints, state->ones_dual_d, state->primal_slack); - state->dual_objective_value = (base_dual_objective + dual_slack_sum) / - (state->constraint_bound_rescaling * state->objective_vector_rescaling) + - state->objective_constant; -} diff --git a/test/test_infeasibility_false_positive.py b/test/test_infeasibility_false_positive.py new file mode 100644 index 0000000..4dfb12b --- /dev/null +++ b/test/test_infeasibility_false_positive.py @@ -0,0 +1,117 @@ +# Copyright 2025 Haihao Lu +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Regression tests for false infeasibility claims. + +The PDLP ray criterion (constraint violation relative to the ray objective) +is not scale invariant: on a bounded, feasible LP whose objective vector or +matrix contains huge coefficients, a direction that violates the constraints +badly can still look like a certificate. These LPs must solve to OPTIMAL with +the default infeasibility tolerance, with and without presolve. +""" + +import numpy as np +import pytest +from cupdlpx import Model, PDLP + + +def _solve(c, A, l, u, lb, ub, presolve): + # Model defaults variable bounds to (-inf, inf); these LPs are stated with x >= 0. + n = len(c) + if lb is None: + lb = np.zeros(n) + if ub is None: + ub = np.full(n, np.inf) + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=presolve, TimeLimit=30.0) + model.optimize() + return model + + +@pytest.mark.parametrize("presolve", [False, True]) +def test_huge_objective_with_dummy_equality(presolve, atol): + """ + Minimize -1e10*x - 1.001e10*y + Subject to x + y <= 2 + 1e10*z == 0 + x, y, z >= 0 + Optimum (0, 2, 0) with objective -2.002e10. The direction (0, 1, 0) + improves the objective by 1e10 per unit while violating x + y <= 0 by 1; + the large unrelated coefficient on z must not hide that violation. + """ + c = np.array([-1e10, -1.001e10, 0.0]) + A = np.array([[1.0, 1.0, 0.0], [0.0, 0.0, 1e10]]) + l = np.array([-np.inf, 0.0]) + u = np.array([2.0, 0.0]) + model = _solve(c, A, l, u, None, None, presolve) + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" + assert abs(model.ObjVal / -2.002e10 - 1.0) < atol + + +@pytest.mark.parametrize("presolve", [False, True]) +def test_huge_objective_and_huge_row(presolve, atol): + """ + Same LP, but the large coefficients sit in a real constraint: + Minimize -1e10*x - 1.001e10*y + Subject to x + y <= 2 + -1e10*x + 1e10*y >= 0 + x, y >= 0 + Optimum (0, 2) with objective -2.002e10. + """ + c = np.array([-1e10, -1.001e10]) + A = np.array([[1.0, 1.0], [-1e10, 1e10]]) + l = np.array([-np.inf, 0.0]) + u = np.array([2.0, np.inf]) + model = _solve(c, A, l, u, None, None, presolve) + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" + assert abs(model.ObjVal / -2.002e10 - 1.0) < atol + + +@pytest.mark.parametrize("presolve", [False, True]) +@pytest.mark.parametrize("seed", range(6)) +def test_random_bounded_lp_with_dummy_equality(presolve, seed, atol): + """ + Random bounded LPs: positive constraint coefficients, x = 1 feasible, + objective scaled by 1e10, plus an unrelated equality with a 1e12 + coefficient. All are feasible and bounded. + """ + rng = np.random.default_rng(seed) + m, n = 6, 4 + A_rows = rng.uniform(0.1, 2.0, (m, n)) + b = A_rows @ np.ones(n) + c = np.concatenate([-rng.uniform(0.1, 2.0, n) * 1e10, [0.0]]) + A = np.zeros((m + 1, n + 1)) + A[:m, :n] = A_rows + A[m, n] = 1e12 + l = np.concatenate([np.full(m, -np.inf), [0.0]]) + u = np.concatenate([b, [0.0]]) + model = _solve(c, A, l, u, None, None, presolve) + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" + # a feasible point (x = 1, z = 0) bounds the optimum from above + assert model.ObjVal <= c[:n].sum() * (1.0 - atol) + atol * abs(c[:n].sum()) + + +@pytest.mark.parametrize("presolve", [False, True]) +@pytest.mark.parametrize("scale", [1e-6, 1.0, 1e8]) +def test_infeasible_lp_detection_is_scale_invariant(presolve, scale, base_lp_data): + """ + The infeasible LP from test_infeasible_unbounded.py with the objective + scaled must still be declared primal infeasible with the default tolerance. + """ + c, A, l, u, lb, ub = base_lp_data + l, u = l.copy(), u.copy() + l[0], u[0] = 10, 10 + model = _solve(c * scale, A, l, u, lb, ub, presolve) + assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" diff --git a/test/test_infeasible_unbounded.py b/test/test_infeasible_unbounded.py index 60eefe9..bfadcca 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -74,4 +74,30 @@ def test_unbounded_lp(base_lp_data, atol): assert hasattr(model, "Status"), "Model.Status not exposed." assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check primal ray - assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}" \ No newline at end of file + assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}" + + +def test_infeasible_tol_zero_disables_detection(base_lp_data): + """ + With InfeasibleTol=0 the solver never declares infeasibility and runs the + infeasible LP into its iteration limit instead. + """ + c, A, l, u, lb, ub = base_lp_data + l, u = l.copy(), u.copy() + l[0], u[0] = 10, 10 + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False, InfeasibleTol=0.0, IterationLimit=2000) + model.optimize() + assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" + +def test_infeasibility_detection_flag_disables_detection(base_lp_data): + """ + InfeasibilityDetection=False has the same effect as InfeasibleTol=0. + """ + c, A, l, u, lb, ub = base_lp_data + l, u = l.copy(), u.copy() + l[0], u[0] = 10, 10 + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False, InfeasibilityDetection=False, IterationLimit=2000) + model.optimize() + assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" diff --git a/test/test_interface.c b/test/test_interface.c index f62e773..a3935c0 100644 --- a/test/test_interface.c +++ b/test/test_interface.c @@ -227,7 +227,7 @@ int main() // Test 9: GPU solver path (presolve disabled) -- forces hipBLAS/hipSPARSE execution printf("\n=== Test 9: CSR Matrix (presolve disabled, GPU solver) ===\n"); { - lp_problem_t *prob9 = create_lp_problem(c, &A_csr, l, u, NULL, NULL, NULL); + lp_problem_t *prob9 = create_lp_problem(c, &A_csr, l, u, NULL, NULL, NULL, NULL); if (!prob9) { fprintf(stderr, "[test] create_lp_problem failed for Test 9.\n"); From 8a9c06a46e28c1891590e6dc99b1261efaeee563 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 04:05:36 +0000 Subject: [PATCH 02/15] Add infeasible_objective_margin parameter The relative margin that the objective term of a ray must clear to count as an infeasibility certificate (t3 >= delta) was a compile-time constant of 1e-8. Expose it as a solver parameter with the same default: --infeasible_obj_margin on the CLI, InfeasibleObjMargin / infeasible_objective_margin in Python. It must be non-negative; 0 reduces the condition to a plain sign test. Detection results are unchanged. --- README.md | 1 + include/cupdlpx_types.h | 1 + internal/infeasibility.h | 4 +++- python/README.md | 3 ++- python/cupdlpx/PDLP.py | 1 + python/cupdlpx/model.py | 6 +++++- python_bindings/_core_bindings.cpp | 4 ++++ src/cli.c | 9 +++++++++ src/infeasibility.cu | 30 ++++++++++++++---------------- src/solver.cu | 2 +- src/utils.cu | 4 ++++ test/test_api_surface.py | 4 ++++ test/test_infeasible_unbounded.py | 15 +++++++++++++++ 13 files changed, 64 insertions(+), 20 deletions(-) diff --git a/README.md b/README.md index 98cdc4d..a102d93 100644 --- a/README.md +++ b/README.md @@ -118,6 +118,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | | `--no_presolve` | `flag` | Disable presolve | `enabled` | | `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate | `1e-8` | +| `--infeasible_obj_margin` | `double` | Smallest relative objective margin of an accepted ray certificate (`0`: plain sign test) | `1e-8` | | `--no_infeasibility_detection` | `flag` | Disable infeasibility detection | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 25b9386..a4af175 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -88,6 +88,7 @@ extern "C" double eps_feasible_relative; double eps_feas_polish_relative; double eps_infeasible_relative; + double infeasible_objective_margin; double time_sec_limit; int iteration_limit; } termination_criteria_t; diff --git a/internal/infeasibility.h b/internal/infeasibility.h index acf2001..338e9ff 100644 --- a/internal/infeasibility.h +++ b/internal/infeasibility.h @@ -44,7 +44,9 @@ extern "C" void compute_infeasibility_device(pdhg_solver_state_t *state); - void compute_infeasibility_host(pdhg_solver_state_t *state); + // objective_margin is the relative margin the objective term of a ray must + // clear (termination_criteria_t::infeasible_objective_margin). + void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin); #ifdef __cplusplus } diff --git a/python/README.md b/python/README.md index c7debc5..a070044 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. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix), and its objective term has the right sign with a fixed relative margin. `0` disables infeasibility detection. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. `0` disables infeasibility detection. | +| `InfeasibleObjMargin` | `infeasible_objective_margin` | float | `1e-8` | Relative margin the objective term of a ray must clear to count as a certificate: `-c^T d >= margin * |c|^T |d|` for a primal ray `d`, and the same for the dual objective term of a dual ray. It only has to clear rounding noise; `0` reduces it to a sign test. | | `InfeasibilityDetection` | `infeasibility_detection` | bool | `True` | Whether to look for infeasibility certificates at all. Setting `InfeasibleTol` to `0` disables detection too. | | `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. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 412426f..dcc5a2b 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -42,6 +42,7 @@ "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", "InfeasibleTol": "eps_infeasible_relative", + "InfeasibleObjMargin": "infeasible_objective_margin", "InfeasibilityDetection": "infeasibility_detection", # scaling / step size "GeoMeanIters": "geometric_mean_iterations", diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index a28b857..da60e8c 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -55,6 +55,7 @@ "eps_optimal_relative", "eps_feasible_relative", "eps_infeasible_relative", + "infeasible_objective_margin", "time_sec_limit", "pock_chambolle_alpha", "artificial_restart_threshold", @@ -78,7 +79,10 @@ } ) # eps_infeasible_relative = 0 also disables detection, like infeasibility_detection=False -_NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol", "eps_infeasible_relative"}) +# infeasible_objective_margin = 0 reduces the objective condition of a certificate to a sign test +_NONNEGATIVE_FLOAT_PARAMS = frozenset( + {"time_sec_limit", "matrix_zero_tol", "eps_infeasible_relative", "infeasible_objective_margin"} +) _STRING_PARAMS = frozenset({"optimality_norm"}) # int params are stored as C int32 on the backend diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index 35f001f..df56b47 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -274,6 +274,7 @@ static py::dict get_default_params_py() d["eps_optimal_relative"] = p.termination_criteria.eps_optimal_relative; d["eps_feasible_relative"] = p.termination_criteria.eps_feasible_relative; d["eps_infeasible_relative"] = p.termination_criteria.eps_infeasible_relative; + d["infeasible_objective_margin"] = p.termination_criteria.infeasible_objective_margin; // limits d["time_sec_limit"] = p.termination_criteria.time_sec_limit; @@ -387,6 +388,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getf("eps_optimal_relative", p->termination_criteria.eps_optimal_relative); getf("eps_feasible_relative", p->termination_criteria.eps_feasible_relative); getf("eps_infeasible_relative", p->termination_criteria.eps_infeasible_relative); + getf("infeasible_objective_margin", p->termination_criteria.infeasible_objective_margin); // limits getf("time_sec_limit", p->termination_criteria.time_sec_limit); @@ -442,6 +444,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p if (p->termination_criteria.eps_infeasible_relative < 0.0) throw std::invalid_argument( "eps_infeasible_relative must be non-negative (0 disables infeasibility detection)."); + if (p->termination_criteria.infeasible_objective_margin < 0.0) + throw std::invalid_argument("infeasible_objective_margin must be non-negative (0 is a plain sign test)."); if (p->termination_criteria.eps_feas_polish_relative <= 0.0) throw std::invalid_argument("eps_feas_polish_relative must be positive."); if (p->sv_tol <= 0.0) diff --git a/src/cli.c b/src/cli.c index 2ae78b5..15f4960 100644 --- a/src/cli.c +++ b/src/cli.c @@ -215,6 +215,11 @@ void print_usage(const char *prog_name) "largest relative\n" " constraint residual of an accepted ray " "certificate (default: 1e-8; 0 also disables detection).\n"); + fprintf(stderr, + " --infeasible_obj_margin Relative margin the objective term of a ray " + "must clear\n" + " to count as a certificate (default: 1e-8; 0 is a " + "plain sign test).\n"); fprintf(stderr, " --opt_norm " "Norm for optimality criteria: l2 or linf (default: l2).\n"); @@ -246,6 +251,7 @@ int main(int argc, char *argv[]) {"eps_feas", required_argument, 0, 1004}, {"eps_feas_polish", required_argument, 0, 1006}, {"eps_infeasible", required_argument, 0, 1019}, + {"infeasible_obj_margin", required_argument, 0, 1021}, {"feasibility_polishing", no_argument, 0, 'f'}, {"geo_mean_iter", required_argument, 0, 1018}, {"l_inf_ruiz_iter", required_argument, 0, 1007}, @@ -294,6 +300,9 @@ int main(int argc, char *argv[]) case 1019: // --eps_infeasible params.termination_criteria.eps_infeasible_relative = atof(optarg); break; + case 1021: // --infeasible_obj_margin + params.termination_criteria.infeasible_objective_margin = atof(optarg); + break; case 'f': // --feasibility_polishing params.feasibility_polishing = true; break; diff --git a/src/infeasibility.cu b/src/infeasibility.cu index a3f1e71..eb57c4e 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -382,12 +382,6 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) free(ctx); } -// Threshold on the objective-side margin of a certificate (t3 below). It only -// has to clear rounding noise and is deliberately independent of -// eps_infeasible, so that loosening the tolerance can only make detection -// easier. -#define INFEAS_OBJECTIVE_MARGIN 1e-8 - // On an infeasible problem the PDHG iterates diverge along the infimal // displacement vector v, whose primal or dual part is an infeasibility // certificate (Applegate et al. 2021; Lu and Yang 2024 for Halpern PDHG). Three @@ -404,9 +398,12 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) // normalized by its infinity norm in the original space; its primal part d is // a certificate of dual infeasibility and its dual part y one of primal // infeasibility when -// the objective term has the right sign with a relative margin, -// -c^T d >= INFEAS_OBJECTIVE_MARGIN |c|^T |d| (t3, primal side) -// Omega(y) >= INFEAS_OBJECTIVE_MARGIN sum |terms| (t3, dual side) +// the objective term has the right sign with a relative margin delta (the +// parameter infeasible_objective_margin; it only has to clear rounding noise +// and is deliberately separate from eps_infeasible, so that loosening the +// tolerance can only make detection easier), +// -c^T d >= delta |c|^T |d| (t3, primal side) +// Omega(y) >= delta sum |terms| (t3, dual side) // and the violation of the homogeneous constraints, relative to the // activity the ray produces, is at most eps_infeasible: // viol(A d) / max_i (|A||d|)_i <= eps (t2, primal side) @@ -457,9 +454,10 @@ void compute_infeasibility_device(pdhg_solver_state_t *state) } // Statistics of every candidate in the original space and the objective-side -// margin test, from the fetched first-round scalars. Returns whether any -// candidate passed the test, in which case the second round is needed. -static bool evaluate_objective_margins(pdhg_solver_state_t *state) +// margin test with the relative margin objective_margin, from the fetched +// first-round scalars. Returns whether any candidate passed the test, in which +// case the second round is needed. +static bool evaluate_objective_margins(pdhg_solver_state_t *state, double objective_margin) { infeas_ctx_t *ctx = state->infeas_ctx; const double bc = state->constraint_bound_rescaling, oc = state->objective_vector_rescaling; @@ -480,9 +478,9 @@ static bool evaluate_objective_margins(pdhg_solver_state_t *state) : 0.0; s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max_orig / oc / dual_norm : 0.0; ctx->need_primal[c] = s->primal_ray_linear_objective < 0.0 && s->primal_objective_abs > 0.0 && - -s->primal_ray_linear_objective >= INFEAS_OBJECTIVE_MARGIN * s->primal_objective_abs; + -s->primal_ray_linear_objective >= objective_margin * s->primal_objective_abs; ctx->need_dual[c] = s->dual_ray_objective > 0.0 && s->dual_objective_abs > 0.0 && - s->dual_ray_objective >= INFEAS_OBJECTIVE_MARGIN * s->dual_objective_abs; + s->dual_ray_objective >= objective_margin * s->dual_objective_abs; any = any || ctx->need_primal[c] || ctx->need_dual[c]; } return any; @@ -578,9 +576,9 @@ static void select_best_candidates(pdhg_solver_state_t *state) // statistics of the state. Host arithmetic only, unless some candidate passes // its objective margin: then the second round runs here, with one more copy // and synchronization. -void compute_infeasibility_host(pdhg_solver_state_t *state) +void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin) { - if (evaluate_objective_margins(state)) + if (evaluate_objective_margins(state, objective_margin)) { launch_abs_products(state); fetch_device_scalars(state); diff --git a/src/solver.cu b/src/solver.cu index 48bfca1..2040e26 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -221,7 +221,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * compute_residual_host(state, params->optimality_norm); if (detect_infeasibility) { - compute_infeasibility_host(state); + compute_infeasibility_host(state, params->termination_criteria.infeasible_objective_margin); } else { diff --git a/src/utils.cu b/src/utils.cu index 059f8ba..af5ffac 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -277,6 +277,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->termination_criteria.iteration_limit = INT32_MAX; params->termination_criteria.eps_feas_polish_relative = 1e-6; params->termination_criteria.eps_infeasible_relative = 1e-8; + params->termination_criteria.infeasible_objective_margin = 1e-8; params->restart_params.artificial_restart_threshold = 0.36; params->restart_params.sufficient_reduction_for_restart = 0.2; @@ -588,6 +589,9 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro PRINT_DIFF_DBL("eps_infeasible_relative", params->termination_criteria.eps_infeasible_relative, default_params.termination_criteria.eps_infeasible_relative); + PRINT_DIFF_DBL("infeasible_objective_margin", + params->termination_criteria.infeasible_objective_margin, + default_params.termination_criteria.infeasible_objective_margin); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); PRINT_DIFF_BOOL("infeasibility_detection", params->infeasibility_detection, default_params.infeasibility_detection); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); diff --git a/test/test_api_surface.py b/test/test_api_surface.py index b418ba5..191ae1e 100644 --- a/test/test_api_surface.py +++ b/test/test_api_surface.py @@ -243,6 +243,8 @@ def test_param_value_validation(base_lp_data): model.setParam("OptimalityTol", 0.0) with pytest.raises(ValueError): model.setParam("TimeLimit", -1.0) + with pytest.raises(ValueError): + model.setParam("InfeasibleObjMargin", -1e-8) with pytest.raises(ValueError): model.setParam("OptimalityNorm", "l1") @@ -259,6 +261,8 @@ def test_direct_core_param_value_validation(base_lp_data): solve_once(A, c, None, lb, ub, l, u, params={"time_sec_limit": float("nan")}) with pytest.raises(ValueError): solve_once(A, c, None, lb, ub, l, u, params={"optimality_norm": "l1"}) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"infeasible_objective_margin": -1.0}) def test_direct_core_model_data_validation(base_lp_data): diff --git a/test/test_infeasible_unbounded.py b/test/test_infeasible_unbounded.py index bfadcca..aa4603c 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -101,3 +101,18 @@ def test_infeasibility_detection_flag_disables_detection(base_lp_data): model.setParams(OutputFlag=False, Presolve=False, InfeasibilityDetection=False, IterationLimit=2000) model.optimize() assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" + +def test_infeasible_objective_margin_blocks_detection(base_lp_data): + """ + InfeasibleObjMargin is the relative margin the objective term of a ray must + clear. A margin above 1 can never be met, because the objective term is at + most the sum of the absolute values of its terms, so the infeasible LP of + test_infeasible_lp runs into its iteration limit instead. + """ + c, A, l, u, lb, ub = base_lp_data + l, u = l.copy(), u.copy() + l[0], u[0] = 10, 10 + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False, InfeasibleObjMargin=2.0, IterationLimit=2000) + model.optimize() + assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" From cc317e05db8127958414735868010a7fab2cf1b7 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 04:09:34 +0000 Subject: [PATCH 03/15] Evaluate the certificate ratio t2 in the units of the original problem The infeasibility check accepts a ray when max viol / max |A||d| is at most eps_infeasible. Both maxima were taken in the equilibrated space of the preconditioner; take them in the units of the original problem instead: eps is then a relative perturbation of the matrix as given, independent of the rescaling settings, and judged in the same units as the optimality residuals. OSQP, SCS and PDLP test their certificates in the original space as well. On the 334 MIPLIB infeasible/unbounded instances the detected set is unchanged (318 with and without presolve) and detection takes on average 6% (presolve off) and 3% (on) fewer iterations. No false claim on 413 feasible instances, but the closest run moves from 3.2e-5 to 1.08e-8 (ns1688926): a row with very large coefficients that the projected ray satisfies by sign can dominate the activity and hide a violated row, so test_huge_objective_and_huge_row with presolve off is now an expected failure and the Python README asks for coefficient ranges within about 1/InfeasibleTol. Details in docs/infeasibility_detection. --- README.md | 2 +- internal/internal_types.h | 32 ++++++++++-------- python/README.md | 2 +- src/infeasibility.cu | 40 +++++++++++++++-------- test/test_infeasibility_false_positive.py | 16 ++++++++- 5 files changed, 62 insertions(+), 30 deletions(-) diff --git a/README.md b/README.md index a102d93..3c1739e 100644 --- a/README.md +++ b/README.md @@ -117,7 +117,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` | | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | | `--no_presolve` | `flag` | Disable presolve | `enabled` | -| `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate | `1e-8` | +| `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate, relative to the activity `|A||x|` of the ray on the problem as given | `1e-8` | | `--infeasible_obj_margin` | `double` | Smallest relative objective margin of an accepted ray certificate (`0`: plain sign test) | `1e-8` | | `--no_infeasibility_detection` | `flag` | Disable infeasibility detection | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | diff --git a/internal/internal_types.h b/internal/internal_types.h index fe73a3a..614642c 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -26,8 +26,10 @@ limitations under the License. #include // The scalars a termination check produces on the device, copied to the host -// in one go; see pdhg_solver_state_t::scalars_h. All of them are in the -// rescaled space; the host halves unscale. Name parts: res = residual, +// in one go; see pdhg_solver_state_t::scalars_h. The optimality scalars are in +// the rescaled space and the host halves unscale; the certificate scalars of +// the infeasibility check are formed in original units on the device, see +// primal_ray_scalars_t. Name parts: res = residual, // sq = sum of squares, max = infinity norm, dot = inner product. // // Each member is the destination of one reduction (reduce.h). @@ -40,20 +42,23 @@ limitations under the License. // infinity norm: the certificate criterion is defined with it, and the // reported statistics follow PDLP's normalization. // -// The violation exists twice because the two are different maxima over the -// rows. The scaled one is the numerator of the certificate error -// viol_max_scaled / product_abs_max, evaluated in the rescaled space because -// the ratio is not invariant to non-uniform row scaling and the rescaling has -// equilibrated the rows. The orig one is in the units of the original problem -// and is only reported (max_primal_ray_infeasibility and its dual twin). +// The violation of the sign conditions and the activity |A||d| are both in +// the units of the original problem (the row and column rescaling factors are +// multiplied back on the device), so the certificate error +// viol_max / product_abs_max is a relative perturbation of the matrix as +// given and does not depend on the preconditioner. The two scalar rescaling +// factors beta_b, beta_c cancel in the ratio and are left out; the reported +// statistic divides viol_max by them and by max_abs. Neither maximum is +// invariant to a non-uniform row scaling of the input, like the optimality +// residuals; OSQP, SCS and PDLP test their certificates in the original space +// as well. typedef struct { double max_abs; // ||d||_inf in the original space; normalizes the reported statistics double obj; // c^T d, the objective term double obj_abs; // |c|^T |d|, the same with absolute values - double viol_max_scaled; // max violation of the sign conditions on A d - double product_abs_max; // max_i (|A||d|)_i, scaled like the violation above (second round) - double viol_max_orig; // the violation in original units + double viol_max; // max violation of the sign conditions on A d, original units + double product_abs_max; // max_i (|A||d|)_i, original units (second round) } primal_ray_scalars_t; typedef struct @@ -63,9 +68,8 @@ typedef struct double obj_var_bounds; // l_v^T r^+ - u_v^T r^-, the variable-bound part double obj_con_bounds_abs; // |l_c|^T y^+ + |u_c|^T y^-, the constraint-bound part with absolute values double obj_var_bounds_abs; // |l_v|^T r^+ + |u_v|^T r^-, the variable-bound part with absolute values - double viol_max_scaled; // max violation of the sign conditions on r - double product_abs_max; // max_j (|A^T||y|)_j, scaled like the violation above (second round) - double viol_max_orig; // the violation in original units + double viol_max; // max violation of the sign conditions on r, original units + double product_abs_max; // max_j (|A^T||y|)_j, original units (second round) } dual_ray_scalars_t; typedef struct diff --git a/python/README.md b/python/README.md index a070044..71037c6 100644 --- a/python/README.md +++ b/python/README.md @@ -163,7 +163,7 @@ 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. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. `0` disables infeasibility detection. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix as given, in the units of the original problem and independent of the solver's internal rescaling), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. The maximum is over rows, so a model whose coefficient magnitudes span more than about `1/InfeasibleTol` can produce a false certificate; scale such models. `0` disables infeasibility detection. | | `InfeasibleObjMargin` | `infeasible_objective_margin` | float | `1e-8` | Relative margin the objective term of a ray must clear to count as a certificate: `-c^T d >= margin * |c|^T |d|` for a primal ray `d`, and the same for the dual objective term of a dual ray. It only has to clear rounding noise; `0` reduces it to a sign test. | | `InfeasibilityDetection` | `infeasibility_detection` | bool | `True` | Whether to look for infeasibility certificates at all. Setting `InfeasibleTol` to `0` disables detection too. | | `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. | diff --git a/src/infeasibility.cu b/src/infeasibility.cu index eb57c4e..a212c3c 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -183,8 +183,8 @@ __global__ void abs_columns_kernel(double *__restrict__ X, int rows) } // First-round scalars of candidate c from its columns of D, Y and the products -// A D and A^T Y. Violations come in the scaled space (they decide the -// certificate) and in the original space (they are reported). +// A D and A^T Y. The violations are in original units; they decide the +// certificate and are reported. static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -219,12 +219,11 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * d[j]; }, &out->primal_ray.obj); reduce_sum( state, n, [=] __host__ __device__(int j) { return fabs(objective[j] * d[j]); }, &out->primal_ray.obj_abs); - reduce_max(state, m, row_violation, &out->primal_ray.viol_max_scaled); reduce_max( state, m, [=] __host__ __device__(int i) { return row_violation(i) * constraint_rescaling[i]; }, - &out->primal_ray.viol_max_orig); + &out->primal_ray.viol_max); // the dual ray y and its reduced cost r = -A^T y: the norm of y in the // original space, the two parts of the dual objective term with and @@ -266,12 +265,11 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) n, [=] __host__ __device__(int j) { return fabs(r_lower_term(j)) + fabs(r_upper_term(j)); }, &out->dual_ray.obj_var_bounds_abs); - reduce_max(state, n, reduced_cost_violation, &out->dual_ray.viol_max_scaled); reduce_max( state, n, [=] __host__ __device__(int j) { return reduced_cost_violation(j) * variable_rescaling[j]; }, - &out->dual_ray.viol_max_orig); + &out->dual_ray.viol_max); } static void spmm_buffer( @@ -411,7 +409,21 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) // Both conditions are relative to the data: a ray with t2 <= eps is an exact // certificate for a problem whose matrix differs from A by a relative amount // eps, and its objective condition survives relative perturbations of c below -// the margin. The plain PDLP ratio viol / (-c^T d) is not used: it is not scale +// the margin. Both ratios are formed in the units of the original problem: the +// violation and the activity are multiplied back by the row (column) rescaling +// factors on the device, so eps refers to the matrix as given and the verdict +// does not depend on the preconditioner. OSQP, SCS and PDLP test their +// certificates in the original space as well. The maximum over rows is not +// invariant to a non-uniform row scaling of the input, like the optimality +// residuals, and a row with very large coefficients that the projected ray +// satisfies by its sign can dominate the activity and hide a violated row of +// ordinary size (test_huge_objective_and_huge_row, coefficient range 1e10, is +// the known case). On the MIPLIB infeasible set the original and the +// equilibrated space detect the same instances, the original one on average +// a few percent earlier; on 413 feasible instances neither makes a false +// claim, but the closest run is 1.08e-8 in the original space against 3.2e-5 +// in the equilibrated one (docs/infeasibility_detection). The plain PDLP +// ratio viol / (-c^T d) is not used: it is not scale // invariant, and a large objective vector lets junk directions pass. The // certificate error stored in the state is t2 when the margin holds and // INFINITY otherwise; separately for the primal and the dual ray, the candidate @@ -470,13 +482,13 @@ static bool evaluate_objective_margins(pdhg_solver_state_t *state, double object double primal_norm = v->primal_ray.max_abs, dual_norm = v->dual_ray.max_abs; s->primal_ray_linear_objective = primal_norm > 0.0 ? v->primal_ray.obj / (bc * oc) / primal_norm : 0.0; s->primal_objective_abs = primal_norm > 0.0 ? v->primal_ray.obj_abs / (bc * oc) / primal_norm : 0.0; - s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max_orig / bc / primal_norm : 0.0; + s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max / bc / primal_norm : 0.0; s->dual_ray_objective = dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds + v->dual_ray.obj_var_bounds) / (bc * oc) / dual_norm : 0.0; s->dual_objective_abs = dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds_abs + v->dual_ray.obj_var_bounds_abs) / (bc * oc) / dual_norm : 0.0; - s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max_orig / oc / dual_norm : 0.0; + s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max / oc / dual_norm : 0.0; ctx->need_primal[c] = s->primal_ray_linear_objective < 0.0 && s->primal_objective_abs > 0.0 && -s->primal_ray_linear_objective >= objective_margin * s->primal_objective_abs; ctx->need_dual[c] = s->dual_ray_objective > 0.0 && s->dual_objective_abs > 0.0 && @@ -508,10 +520,11 @@ static void launch_abs_products(pdhg_solver_state_t *state) for (int c = 0; c < K; c++) { const double *p = ctx->prod_m + (size_t)c * m; + const double *constraint_rescaling = state->constraint_rescaling; reduce_max( state, m, - [=] __host__ __device__(int i) { return p[i]; }, + [=] __host__ __device__(int i) { return p[i] * constraint_rescaling[i]; }, &state->scalars_d->infeas[c].primal_ray.product_abs_max); } } @@ -522,10 +535,11 @@ static void launch_abs_products(pdhg_solver_state_t *state) for (int c = 0; c < K; c++) { const double *p = ctx->prod_n + (size_t)c * n; + const double *variable_rescaling = state->variable_rescaling; reduce_max( state, n, - [=] __host__ __device__(int j) { return p[j]; }, + [=] __host__ __device__(int j) { return p[j] * variable_rescaling[j]; }, &state->scalars_d->infeas[c].dual_ray.product_abs_max); } } @@ -545,12 +559,12 @@ static void select_best_candidates(pdhg_solver_state_t *state) double perr = INFINITY, derr = INFINITY; if (ctx->need_primal[c]) { - double viol = v->primal_ray.viol_max_scaled, act = v->primal_ray.product_abs_max; + double viol = v->primal_ray.viol_max, act = v->primal_ray.product_abs_max; perr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); } if (ctx->need_dual[c]) { - double viol = v->dual_ray.viol_max_scaled, act = v->dual_ray.product_abs_max; + double viol = v->dual_ray.viol_max, act = v->dual_ray.product_abs_max; derr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); } if (perr < best_perr) diff --git a/test/test_infeasibility_false_positive.py b/test/test_infeasibility_false_positive.py index 4dfb12b..8c42e3a 100644 --- a/test/test_infeasibility_false_positive.py +++ b/test/test_infeasibility_false_positive.py @@ -60,7 +60,21 @@ def test_huge_objective_with_dummy_equality(presolve, atol): assert abs(model.ObjVal / -2.002e10 - 1.0) < atol -@pytest.mark.parametrize("presolve", [False, True]) +@pytest.mark.parametrize( + "presolve", + [ + pytest.param( + False, + marks=pytest.mark.xfail( + strict=True, + reason="known limitation of the certificate ratio in original units: the 1e10 row is " + "satisfied by the sign of the projected ray, dominates the activity |A||d| and hides the " + "violated row x + y <= 2, so the LP is declared DUAL_INFEASIBLE; presolve removes the row", + ), + ), + True, + ], +) def test_huge_objective_and_huge_row(presolve, atol): """ Same LP, but the large coefficients sit in a real constraint: From 75bded1a30e7825dd71a3c69668a1e35a6dd669a Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 18:22:24 +0000 Subject: [PATCH 04/15] change back scaled LP to t2 --- README.md | 2 +- internal/internal_types.h | 38 +++++++++-------- python/README.md | 2 +- src/infeasibility.cu | 51 +++++++++++------------ test/test_infeasibility_false_positive.py | 16 +------ 5 files changed, 48 insertions(+), 61 deletions(-) diff --git a/README.md b/README.md index 3c1739e..8d978ee 100644 --- a/README.md +++ b/README.md @@ -117,7 +117,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` | | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | | `--no_presolve` | `flag` | Disable presolve | `enabled` | -| `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate, relative to the activity `|A||x|` of the ray on the problem as given | `1e-8` | +| `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate, relative to the activity `|A||x|` of the ray on the rescaled problem | `1e-8` | | `--infeasible_obj_margin` | `double` | Smallest relative objective margin of an accepted ray certificate (`0`: plain sign test) | `1e-8` | | `--no_infeasibility_detection` | `flag` | Disable infeasibility detection | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | diff --git a/internal/internal_types.h b/internal/internal_types.h index 614642c..4b33c13 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -26,10 +26,9 @@ limitations under the License. #include // The scalars a termination check produces on the device, copied to the host -// in one go; see pdhg_solver_state_t::scalars_h. The optimality scalars are in -// the rescaled space and the host halves unscale; the certificate scalars of -// the infeasibility check are formed in original units on the device, see -// primal_ray_scalars_t. Name parts: res = residual, +// in one go; see pdhg_solver_state_t::scalars_h. They are in the rescaled +// space and the host halves unscale, except the reported ray violations +// (viol_max_orig below). Name parts: res = residual, // sq = sum of squares, max = infinity norm, dot = inner product. // // Each member is the destination of one reduction (reduce.h). @@ -42,23 +41,25 @@ limitations under the License. // infinity norm: the certificate criterion is defined with it, and the // reported statistics follow PDLP's normalization. // -// The violation of the sign conditions and the activity |A||d| are both in -// the units of the original problem (the row and column rescaling factors are -// multiplied back on the device), so the certificate error -// viol_max / product_abs_max is a relative perturbation of the matrix as -// given and does not depend on the preconditioner. The two scalar rescaling -// factors beta_b, beta_c cancel in the ratio and are left out; the reported -// statistic divides viol_max by them and by max_abs. Neither maximum is -// invariant to a non-uniform row scaling of the input, like the optimality -// residuals; OSQP, SCS and PDLP test their certificates in the original space -// as well. +// The violation exists twice because the two are different maxima over the +// rows. The scaled one is the numerator of the certificate error +// viol_max_scaled / product_abs_max, evaluated in the rescaled space: the +// ratio of the two maxima is not invariant to non-uniform row scaling, and in +// original units a row with very large coefficients that the ray satisfies by +// its sign dominates the activity and hides a violated row of ordinary size, +// which turns junk directions of feasible problems into certificates (see the +// comment above compute_infeasibility_device). The rescaling has equilibrated +// the rows, so the maxima compare like with like. The orig one is in the +// units of the original problem and is only reported +// (max_primal_ray_infeasibility and its dual twin). typedef struct { double max_abs; // ||d||_inf in the original space; normalizes the reported statistics double obj; // c^T d, the objective term double obj_abs; // |c|^T |d|, the same with absolute values - double viol_max; // max violation of the sign conditions on A d, original units - double product_abs_max; // max_i (|A||d|)_i, original units (second round) + double viol_max_scaled; // max violation of the sign conditions on A d, rescaled space + double product_abs_max; // max_i (|A||d|)_i, scaled like the violation above (second round) + double viol_max_orig; // the violation in original units, reported only } primal_ray_scalars_t; typedef struct @@ -68,8 +69,9 @@ typedef struct double obj_var_bounds; // l_v^T r^+ - u_v^T r^-, the variable-bound part double obj_con_bounds_abs; // |l_c|^T y^+ + |u_c|^T y^-, the constraint-bound part with absolute values double obj_var_bounds_abs; // |l_v|^T r^+ + |u_v|^T r^-, the variable-bound part with absolute values - double viol_max; // max violation of the sign conditions on r, original units - double product_abs_max; // max_j (|A^T||y|)_j, original units (second round) + double viol_max_scaled; // max violation of the sign conditions on r, rescaled space + double product_abs_max; // max_j (|A^T||y|)_j, scaled like the violation above (second round) + double viol_max_orig; // the violation in original units, reported only } dual_ray_scalars_t; typedef struct diff --git a/python/README.md b/python/README.md index 71037c6..5434aec 100644 --- a/python/README.md +++ b/python/README.md @@ -163,7 +163,7 @@ 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. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix as given, in the units of the original problem and independent of the solver's internal rescaling), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. The maximum is over rows, so a model whose coefficient magnitudes span more than about `1/InfeasibleTol` can produce a false certificate; scale such models. `0` disables infeasibility detection. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix, evaluated after the solver's row and column equilibration so that the scaling of the input does not decide the verdict), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. `0` disables infeasibility detection. | | `InfeasibleObjMargin` | `infeasible_objective_margin` | float | `1e-8` | Relative margin the objective term of a ray must clear to count as a certificate: `-c^T d >= margin * |c|^T |d|` for a primal ray `d`, and the same for the dual objective term of a dual ray. It only has to clear rounding noise; `0` reduces it to a sign test. | | `InfeasibilityDetection` | `infeasibility_detection` | bool | `True` | Whether to look for infeasibility certificates at all. Setting `InfeasibleTol` to `0` disables detection too. | | `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. | diff --git a/src/infeasibility.cu b/src/infeasibility.cu index a212c3c..d1983cc 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -183,8 +183,8 @@ __global__ void abs_columns_kernel(double *__restrict__ X, int rows) } // First-round scalars of candidate c from its columns of D, Y and the products -// A D and A^T Y. The violations are in original units; they decide the -// certificate and are reported. +// A D and A^T Y. Violations come in the rescaled space (they decide the +// certificate) and in original units (they are reported). static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -219,11 +219,12 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * d[j]; }, &out->primal_ray.obj); reduce_sum( state, n, [=] __host__ __device__(int j) { return fabs(objective[j] * d[j]); }, &out->primal_ray.obj_abs); + reduce_max(state, m, row_violation, &out->primal_ray.viol_max_scaled); reduce_max( state, m, [=] __host__ __device__(int i) { return row_violation(i) * constraint_rescaling[i]; }, - &out->primal_ray.viol_max); + &out->primal_ray.viol_max_orig); // the dual ray y and its reduced cost r = -A^T y: the norm of y in the // original space, the two parts of the dual objective term with and @@ -265,11 +266,12 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) n, [=] __host__ __device__(int j) { return fabs(r_lower_term(j)) + fabs(r_upper_term(j)); }, &out->dual_ray.obj_var_bounds_abs); + reduce_max(state, n, reduced_cost_violation, &out->dual_ray.viol_max_scaled); reduce_max( state, n, [=] __host__ __device__(int j) { return reduced_cost_violation(j) * variable_rescaling[j]; }, - &out->dual_ray.viol_max); + &out->dual_ray.viol_max_orig); } static void spmm_buffer( @@ -409,20 +411,19 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) // Both conditions are relative to the data: a ray with t2 <= eps is an exact // certificate for a problem whose matrix differs from A by a relative amount // eps, and its objective condition survives relative perturbations of c below -// the margin. Both ratios are formed in the units of the original problem: the -// violation and the activity are multiplied back by the row (column) rescaling -// factors on the device, so eps refers to the matrix as given and the verdict -// does not depend on the preconditioner. OSQP, SCS and PDLP test their -// certificates in the original space as well. The maximum over rows is not -// invariant to a non-uniform row scaling of the input, like the optimality -// residuals, and a row with very large coefficients that the projected ray -// satisfies by its sign can dominate the activity and hide a violated row of -// ordinary size (test_huge_objective_and_huge_row, coefficient range 1e10, is -// the known case). On the MIPLIB infeasible set the original and the -// equilibrated space detect the same instances, the original one on average -// a few percent earlier; on 413 feasible instances neither makes a false -// claim, but the closest run is 1.08e-8 in the original space against 3.2e-5 -// in the equilibrated one (docs/infeasibility_detection). The plain PDLP +// the margin. Both ratios are evaluated in the rescaled space of the +// preconditioner. Row by row viol_i / (|A||d|)_i does not depend on any scaling, +// but the ratio of the two maxima does: in original units a row with very +// large coefficients that the projected ray satisfies by its sign dominates the +// activity and hides a violated row of ordinary size, so that a junk direction +// of a feasible problem passes (test_huge_objective_and_huge_row is the +// two-variable case; on 413 feasible instances the closest run is 1.08e-8 in +// original units against 3.2e-5 after equilibration). The rescaling has +// equilibrated the rows, so there the maxima compare like with like. A genuine +// certificate passes in either space, since its violation goes to zero: on the +// MIPLIB infeasible set both detect the same instances, original units a few +// percent earlier (docs/infeasibility_detection). The statistics reported to +// the user are in original units. The plain PDLP // ratio viol / (-c^T d) is not used: it is not scale // invariant, and a large objective vector lets junk directions pass. The // certificate error stored in the state is t2 when the margin holds and @@ -482,13 +483,13 @@ static bool evaluate_objective_margins(pdhg_solver_state_t *state, double object double primal_norm = v->primal_ray.max_abs, dual_norm = v->dual_ray.max_abs; s->primal_ray_linear_objective = primal_norm > 0.0 ? v->primal_ray.obj / (bc * oc) / primal_norm : 0.0; s->primal_objective_abs = primal_norm > 0.0 ? v->primal_ray.obj_abs / (bc * oc) / primal_norm : 0.0; - s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max / bc / primal_norm : 0.0; + s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max_orig / bc / primal_norm : 0.0; s->dual_ray_objective = dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds + v->dual_ray.obj_var_bounds) / (bc * oc) / dual_norm : 0.0; s->dual_objective_abs = dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds_abs + v->dual_ray.obj_var_bounds_abs) / (bc * oc) / dual_norm : 0.0; - s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max / oc / dual_norm : 0.0; + s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max_orig / oc / dual_norm : 0.0; ctx->need_primal[c] = s->primal_ray_linear_objective < 0.0 && s->primal_objective_abs > 0.0 && -s->primal_ray_linear_objective >= objective_margin * s->primal_objective_abs; ctx->need_dual[c] = s->dual_ray_objective > 0.0 && s->dual_objective_abs > 0.0 && @@ -520,11 +521,10 @@ static void launch_abs_products(pdhg_solver_state_t *state) for (int c = 0; c < K; c++) { const double *p = ctx->prod_m + (size_t)c * m; - const double *constraint_rescaling = state->constraint_rescaling; reduce_max( state, m, - [=] __host__ __device__(int i) { return p[i] * constraint_rescaling[i]; }, + [=] __host__ __device__(int i) { return p[i]; }, &state->scalars_d->infeas[c].primal_ray.product_abs_max); } } @@ -535,11 +535,10 @@ static void launch_abs_products(pdhg_solver_state_t *state) for (int c = 0; c < K; c++) { const double *p = ctx->prod_n + (size_t)c * n; - const double *variable_rescaling = state->variable_rescaling; reduce_max( state, n, - [=] __host__ __device__(int j) { return p[j] * variable_rescaling[j]; }, + [=] __host__ __device__(int j) { return p[j]; }, &state->scalars_d->infeas[c].dual_ray.product_abs_max); } } @@ -559,12 +558,12 @@ static void select_best_candidates(pdhg_solver_state_t *state) double perr = INFINITY, derr = INFINITY; if (ctx->need_primal[c]) { - double viol = v->primal_ray.viol_max, act = v->primal_ray.product_abs_max; + double viol = v->primal_ray.viol_max_scaled, act = v->primal_ray.product_abs_max; perr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); } if (ctx->need_dual[c]) { - double viol = v->dual_ray.viol_max, act = v->dual_ray.product_abs_max; + double viol = v->dual_ray.viol_max_scaled, act = v->dual_ray.product_abs_max; derr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); } if (perr < best_perr) diff --git a/test/test_infeasibility_false_positive.py b/test/test_infeasibility_false_positive.py index 8c42e3a..4dfb12b 100644 --- a/test/test_infeasibility_false_positive.py +++ b/test/test_infeasibility_false_positive.py @@ -60,21 +60,7 @@ def test_huge_objective_with_dummy_equality(presolve, atol): assert abs(model.ObjVal / -2.002e10 - 1.0) < atol -@pytest.mark.parametrize( - "presolve", - [ - pytest.param( - False, - marks=pytest.mark.xfail( - strict=True, - reason="known limitation of the certificate ratio in original units: the 1e10 row is " - "satisfied by the sign of the projected ray, dominates the activity |A||d| and hides the " - "violated row x + y <= 2, so the LP is declared DUAL_INFEASIBLE; presolve removes the row", - ), - ), - True, - ], -) +@pytest.mark.parametrize("presolve", [False, True]) def test_huge_objective_and_huge_row(presolve, atol): """ Same LP, but the large coefficients sit in a real constraint: From d588582502fba1dc80a5a82d036674a58c34bf3b Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 18:29:06 +0000 Subject: [PATCH 05/15] cleanup comments --- internal/internal_types.h | 29 ++++------- internal/reduce.h | 17 +++--- src/infeasibility.cu | 106 +++++++++++++------------------------- 3 files changed, 53 insertions(+), 99 deletions(-) diff --git a/internal/internal_types.h b/internal/internal_types.h index 4b33c13..8214253 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -37,21 +37,13 @@ limitations under the License. // One candidate of the infeasibility check (infeasibility.cu): a primal ray d, // which certifies dual infeasibility, and a dual ray y with reduced cost -// r = -A^T y, which certifies primal infeasibility. Every norm here is the -// infinity norm: the certificate criterion is defined with it, and the -// reported statistics follow PDLP's normalization. +// r = -A^T y, which certifies primal infeasibility. All norms are infinity +// norms. // -// The violation exists twice because the two are different maxima over the -// rows. The scaled one is the numerator of the certificate error -// viol_max_scaled / product_abs_max, evaluated in the rescaled space: the -// ratio of the two maxima is not invariant to non-uniform row scaling, and in -// original units a row with very large coefficients that the ray satisfies by -// its sign dominates the activity and hides a violated row of ordinary size, -// which turns junk directions of feasible problems into certificates (see the -// comment above compute_infeasibility_device). The rescaling has equilibrated -// the rows, so the maxima compare like with like. The orig one is in the -// units of the original problem and is only reported -// (max_primal_ray_infeasibility and its dual twin). +// The violation is reduced twice. viol_max_scaled / product_abs_max is the +// certificate error, with both maxima in the rescaled space, where the rows +// have comparable scales. viol_max_orig is in original units and is only +// reported. typedef struct { double max_abs; // ||d||_inf in the original space; normalizes the reported statistics @@ -199,11 +191,10 @@ typedef struct infeas_ctx_t *infeas_ctx; // detection buffers, created before the main loop when detection is on // Scalars of a termination check (optimality, infeasibility) and of the - // restart; see device_scalars_t above. Every quantity is produced on the - // device by the reductions of reduce.h and the whole struct is copied to - // the host in one go (fetch_device_scalars in termination.cu), so a check - // costs one host round trip instead of one per scalar. Both norms of each residual are always computed, so that L2 - // or L-inf stays a host-side choice. + // restart; see device_scalars_t above. They are reduced on the device + // (reduce.h) and the whole struct is copied to the host at once + // (fetch_device_scalars in termination.cu). Both norms of each residual + // are always computed, so that L2 or L-inf stays a host-side choice. device_scalars_t *scalars_d; device_scalars_t *scalars_h; // pinned, so that the copy is asynchronous void *reduce_temp_d; // temporary storage of the reductions diff --git a/internal/reduce.h b/internal/reduce.h index e1b0a49..7c8fe38 100644 --- a/internal/reduce.h +++ b/internal/reduce.h @@ -16,17 +16,14 @@ // // Every scalar of a check is one cub::DeviceReduce call over a // thrust::transform_iterator: a device lambda computes the quantity from the -// element index on the fly, and the result is written to a member of the -// device-side device_scalars_t (internal_types.h), never to the host. The -// calls only queue kernels and never wait; the whole struct is copied to the -// host once per check (termination.h). -// DeviceReduce uses no atomics ("run-to-run" determinism), so a given problem -// size always produces the same result. +// element index, and the result is written to a member of the device-side +// device_scalars_t (internal_types.h). The calls only queue kernels; the whole +// struct is copied to the host once per check (termination.h). The results +// are run-to-run deterministic. // -// The lambdas capture device pointers by value from local variables (never -// the solver state itself) and are marked __host__ __device__ so that the -// iterator can deduce their result type on the host; nvcc needs -// --extended-lambda for that (set in CMakeLists.txt). +// The lambdas capture device pointers by value from local variables and are +// marked __host__ __device__ so that the iterator can deduce their result type +// on the host; nvcc needs --extended-lambda for that (set in CMakeLists.txt). #ifndef CUPDLPX_REDUCE_H #define CUPDLPX_REDUCE_H diff --git a/src/infeasibility.cu b/src/infeasibility.cu index d1983cc..4894a41 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -15,8 +15,7 @@ limitations under the License. */ // Detection of primal and dual infeasibility from ray certificates. See the -// comment above compute_infeasibility_device for the method, and -// docs/infeasibility_detection/report.tex for its derivation and evaluation. +// comment above compute_infeasibility_device for the method. #include "infeasibility.h" #include "reduce.h" @@ -35,9 +34,8 @@ typedef struct double primal_ray_linear_objective; double max_dual_ray_infeasibility; double dual_ray_objective; - // the objective terms with absolute values, |c|^T|d| and the absolute-value - // sum of the dual objective terms, in the same scale as the objective terms - // above; the objective term relative to this is the margin t3 + // the objective terms with absolute values, in the same scale as the + // objective terms above; the objective margin is measured against them double primal_objective_abs; double dual_objective_abs; } ray_candidate_stats_t; @@ -49,26 +47,16 @@ __global__ void abs_kernel(double *__restrict__ out, const double *__restrict__ out[i] = fabs(x[i]); } -// --------------------------------------------------------------------------- -// Batched infeasibility detection. The three candidate directions are stacked -// column by column into dense n x K and m x K matrices and multiplied by A and -// A^T with one SpMM each, so that the matrix is read once per product instead -// of once per candidate; the products |A||d| and |A^T||y| use the absolute -// values of the matrix entries stored on the same CSR indices. Every scalar of -// a check comes out of the fused passes below (one over the variables and one -// over the constraints per candidate, all quantities at once, deterministic -// two-stage reductions into the solver state's scalar array), so a check -// synchronizes twice at most: once before the objective margin is tested and -// once after the absolute-value products. +// The candidate directions are stacked column by column into dense n x K and +// m x K matrices, so that A and A^T are applied to all of them with one SpMM +// each. |A| and |A^T| share the CSR indices of A and A^T and only store the +// absolute values. // -// Memory: a check keeps four dense K-column matrices, the candidates D and Y -// and the products prod_m and prod_n. The first round fills prod_m with A D -// and prod_n with A^T Y and reduces everything that depends on them, including -// the violations that the second round's ratios need. The second round then no longer -// needs D, Y or the products, so it overwrites D with |D| in place and prod_m -// with |A||D| (likewise on the dual side) and reduces only their maxima. The -// next check rewrites all four matrices. -// --------------------------------------------------------------------------- +// A check has two rounds. The first fills prod_m with A D and prod_n with +// A^T Y and reduces every scalar that depends on them. The second runs only +// when some candidate passes its objective margin: it overwrites D with |D| +// and prod_m with |A||D| (likewise on the dual side) and reduces their maxima. +// The next check rewrites all four matrices. // The scalars of a check are the infeas[] members of device_scalars_t // (internal_types.h), one infeas_candidate_scalars_t per candidate. @@ -78,8 +66,8 @@ struct infeas_ctx double *abs_val_A; double *abs_val_AT; cusparseSpMatDescr_t matA, matAT, matAbsA, matAbsAT; - // column-major, one column per candidate; see "Memory" above for what - // each holds between the two rounds of a check + // column-major, one column per candidate; the comment above says what + // each holds in the two rounds of a check double *D; // n x K: primal candidates, then |D| double *Y; // m x K: dual candidates, then |Y| double *prod_m; // m x K: A D, then |A||D| @@ -382,53 +370,31 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) free(ctx); } -// On an infeasible problem the PDHG iterates diverge along the infimal -// displacement vector v, whose primal or dual part is an infeasibility -// certificate (Applegate et al. 2021; Lu and Yang 2024 for Halpern PDHG). Three -// candidate directions are evaluated at every termination check: -// 0. difference of iterates T(z^k) - z^k, which converges to v at rate -// sqrt(log k)/k within a restart epoch and linearly at the epoch starts; -// 1. the iterate T(z^k) itself, the anchor plus the displacement below; -// 2. displacement since the last restart T(z^k) - z^0, the normalized -// Halpern iterate of Lu and Yang, rate 1/k within an epoch. -// The three are complementary in when they are close to v, and no pair matches -// all three on the benchmarks. +// On an infeasible problem the PDHG iterates diverge along a direction whose +// primal or dual part is an infeasibility certificate (Applegate et al. 2021; +// Lu and Yang 2024 for Halpern PDHG). Three candidate directions are evaluated +// at every termination check: +// 0. the difference of iterates T(z^k) - z^k, +// 1. the iterate T(z^k), +// 2. the displacement since the last restart T(z^k) - z^0. // -// Each candidate is projected onto the sign conditions of a certificate and -// normalized by its infinity norm in the original space; its primal part d is -// a certificate of dual infeasibility and its dual part y one of primal +// Each candidate is projected onto the sign conditions of a certificate. Its +// primal part d certifies dual infeasibility and its dual part y primal // infeasibility when -// the objective term has the right sign with a relative margin delta (the -// parameter infeasible_objective_margin; it only has to clear rounding noise -// and is deliberately separate from eps_infeasible, so that loosening the -// tolerance can only make detection easier), -// -c^T d >= delta |c|^T |d| (t3, primal side) -// Omega(y) >= delta sum |terms| (t3, dual side) +// the objective term has the right sign with the relative margin +// delta = infeasible_objective_margin, +// -c^T d >= delta |c|^T |d| (primal ray) +// Omega(y) >= delta sum |terms of Omega(y)| (dual ray) // and the violation of the homogeneous constraints, relative to the -// activity the ray produces, is at most eps_infeasible: -// viol(A d) / max_i (|A||d|)_i <= eps (t2, primal side) -// viol(A^T y) / max_j (|A^T||y|)_j <= eps (t2, dual side). -// Both conditions are relative to the data: a ray with t2 <= eps is an exact -// certificate for a problem whose matrix differs from A by a relative amount -// eps, and its objective condition survives relative perturbations of c below -// the margin. Both ratios are evaluated in the rescaled space of the -// preconditioner. Row by row viol_i / (|A||d|)_i does not depend on any scaling, -// but the ratio of the two maxima does: in original units a row with very -// large coefficients that the projected ray satisfies by its sign dominates the -// activity and hides a violated row of ordinary size, so that a junk direction -// of a feasible problem passes (test_huge_objective_and_huge_row is the -// two-variable case; on 413 feasible instances the closest run is 1.08e-8 in -// original units against 3.2e-5 after equilibration). The rescaling has -// equilibrated the rows, so there the maxima compare like with like. A genuine -// certificate passes in either space, since its violation goes to zero: on the -// MIPLIB infeasible set both detect the same instances, original units a few -// percent earlier (docs/infeasibility_detection). The statistics reported to -// the user are in original units. The plain PDLP -// ratio viol / (-c^T d) is not used: it is not scale -// invariant, and a large objective vector lets junk directions pass. The -// certificate error stored in the state is t2 when the margin holds and -// INFINITY otherwise; separately for the primal and the dual ray, the candidate -// with the smallest error is kept and its statistics are reported. +// activity the ray produces, is at most eps = eps_infeasible_relative, +// viol(A d) / max_i (|A||d|)_i <= eps (primal ray) +// viol(A^T y) / max_j (|A^T||y|)_j <= eps (dual ray). +// The two maxima of a ratio are taken in the rescaled space, where the rows +// have comparable scales; in original units a row with very large coefficients +// can dominate the activity and hide a violated row. The certificate error +// stored in the state is this ratio when the margin holds and INFINITY +// otherwise. Per side the candidate with the smallest error is kept, and its +// statistics are reported in original units. // // First round, kernels only: the candidates, A D and A^T Y, and every scalar // that depends on them. infeasibility_context_create must have run. From 6a8e04b12170acae9622a103ca2adacc8a6d2f5a Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 19:30:31 +0000 Subject: [PATCH 06/15] clean up comments --- internal/infeasibility.h | 23 ++----- internal/internal_types.h | 29 +++------ internal/reduce.h | 30 +++------ internal/termination.h | 33 +++------- src/feasibility_polish.cu | 10 +-- src/infeasibility.cu | 129 +++++++++++--------------------------- src/reduce.cu | 4 +- src/solver.cu | 9 +-- src/termination.cu | 9 +-- 9 files changed, 78 insertions(+), 198 deletions(-) diff --git a/internal/infeasibility.h b/internal/infeasibility.h index 338e9ff..f17f601 100644 --- a/internal/infeasibility.h +++ b/internal/infeasibility.h @@ -23,29 +23,18 @@ extern "C" { #endif - // Detection of primal and dual infeasibility from candidate rays; see the - // comment above compute_infeasibility_device in infeasibility.cu for the - // method. Same protocol as termination.h: compute_infeasibility_device - // only queues kernels (the candidates, A d, A^T y and their reductions) and - // never waits; the caller fetches the scalars; compute_infeasibility_host - // turns them into the certificate errors of the state. - // - // A check has a second round, because |A||d| and |A^T||y| cost two more - // SpMMs and are only needed when some candidate passes its objective - // margin. compute_infeasibility_host runs it when needed, with one more - // copy and synchronization; otherwise it is host arithmetic only. - - // The context owns every buffer of a check. It is created once before the - // main loop, so that a check never allocates, and destroyed with the solver - // state. + // Same two halves as termination.h. |A||d| and |A^T||y| cost two more SpMMs + // and matter only when some candidate passes its objective margin, so the + // host half computes them itself, with one more synchronization, only then. + + // Owns every buffer of a check; created before the main loop so that a check + // never allocates. infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state); void infeasibility_context_destroy(infeas_ctx_t *ctx); void compute_infeasibility_device(pdhg_solver_state_t *state); - // objective_margin is the relative margin the objective term of a ray must - // clear (termination_criteria_t::infeasible_objective_margin). void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin); #ifdef __cplusplus diff --git a/internal/internal_types.h b/internal/internal_types.h index 8214253..ee5a48c 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -25,25 +25,16 @@ limitations under the License. #include #include -// The scalars a termination check produces on the device, copied to the host -// in one go; see pdhg_solver_state_t::scalars_h. They are in the rescaled -// space and the host halves unscale, except the reported ray violations -// (viol_max_orig below). Name parts: res = residual, -// sq = sum of squares, max = infinity norm, dot = inner product. -// -// Each member is the destination of one reduction (reduce.h). +// Scalars of a termination check, one reduction each, in the rescaled space +// except viol_max_orig. Name parts: res = residual, sq = sum of squares, +// max = infinity norm, dot = inner product. #define INFEAS_NUM_CANDIDATES 3 // candidate rays of one infeasibility check -// One candidate of the infeasibility check (infeasibility.cu): a primal ray d, -// which certifies dual infeasibility, and a dual ray y with reduced cost -// r = -A^T y, which certifies primal infeasibility. All norms are infinity -// norms. -// -// The violation is reduced twice. viol_max_scaled / product_abs_max is the -// certificate error, with both maxima in the rescaled space, where the rows -// have comparable scales. viol_max_orig is in original units and is only -// reported. +// One infeasibility candidate: a primal ray d, which certifies dual +// infeasibility, and a dual ray y with reduced cost r = -A^T y, which certifies +// primal infeasibility. viol_max_scaled / product_abs_max is the certificate +// error; viol_max_orig is only reported. typedef struct { double max_abs; // ||d||_inf in the original space; normalizes the reported statistics @@ -190,11 +181,7 @@ typedef struct void *spmv_ctx; infeas_ctx_t *infeas_ctx; // detection buffers, created before the main loop when detection is on - // Scalars of a termination check (optimality, infeasibility) and of the - // restart; see device_scalars_t above. They are reduced on the device - // (reduce.h) and the whole struct is copied to the host at once - // (fetch_device_scalars in termination.cu). Both norms of each residual - // are always computed, so that L2 or L-inf stays a host-side choice. + // Both norms of each residual are kept so that the norm type is a host-side choice. device_scalars_t *scalars_d; device_scalars_t *scalars_h; // pinned, so that the copy is asynchronous void *reduce_temp_d; // temporary storage of the reductions diff --git a/internal/reduce.h b/internal/reduce.h index 7c8fe38..a21cfe9 100644 --- a/internal/reduce.h +++ b/internal/reduce.h @@ -12,18 +12,11 @@ // See the License for the specific language governing permissions and // limitations under the License. -// Reductions of a termination check. C++ only. -// -// Every scalar of a check is one cub::DeviceReduce call over a -// thrust::transform_iterator: a device lambda computes the quantity from the -// element index, and the result is written to a member of the device-side -// device_scalars_t (internal_types.h). The calls only queue kernels; the whole -// struct is copied to the host once per check (termination.h). The results -// are run-to-run deterministic. -// -// The lambdas capture device pointers by value from local variables and are -// marked __host__ __device__ so that the iterator can deduce their result type -// on the host; nvcc needs --extended-lambda for that (set in CMakeLists.txt). +// Reductions of a termination check (C++ only). Each scalar is one +// cub::DeviceReduce over a thrust::transform_iterator whose lambda computes the +// quantity from the index; the result stays in device memory. The lambdas are +// __host__ __device__ so the iterator can deduce their type, which needs nvcc +// --extended-lambda. #ifndef CUPDLPX_REDUCE_H #define CUPDLPX_REDUCE_H @@ -33,15 +26,13 @@ #include #include -// Temporary storage of cub::DeviceReduce, sized once for the longest vector of -// the problem and shared by every call (reduce.cu). +// Temporary storage of cub::DeviceReduce, sized once and shared by every call. void reduce_temp_storage_create(pdhg_solver_state_t *state); -// Exits with a message when a call would need more temporary storage than -// reduce_temp_storage_create allocated. +// Exits when a call needs more temporary storage than was allocated. void reduce_check_temp_storage(const pdhg_solver_state_t *state, size_t bytes); -// out = sum over i < n of f(i), on the stream, into device memory. +// out = sum over i < n of f(i). template void reduce_sum(pdhg_solver_state_t *state, int n, F f, double *out) { if (n <= 0) @@ -57,9 +48,8 @@ template void reduce_sum(pdhg_solver_state_t *state, int n, F f, do CUDA_CHECK(cub::DeviceReduce::Sum(state->reduce_temp_d, bytes, it, out, n, state->stream)); } -// out = max over i < n of f(i), on the stream, into device memory. Every -// maximum the solver takes is a norm or a violation, hence nonnegative, so an -// empty vector yields 0. +// out = max over i < n of f(i); 0 when n is 0, since every maximum taken here is +// nonnegative. template void reduce_max(pdhg_solver_state_t *state, int n, F f, double *out) { if (n <= 0) diff --git a/internal/termination.h b/internal/termination.h index b72a894..a1dde82 100644 --- a/internal/termination.h +++ b/internal/termination.h @@ -23,22 +23,12 @@ extern "C" { #endif - // The termination check: fixed point error, residuals and objective values - // of the current iterate, and the criteria that decide whether the solver - // stops. The infeasibility check (infeasibility.h) follows the same - // protocol. - // - // Every scalar of a check is produced on the device by the reductions of - // reduce.h and lands in state->scalars_d. A check is computed in halves so - // that it synchronizes once: a _device half only queues kernels on the - // stream and never waits, so the host can queue a whole check while the - // GPU is still running the iterations before it; a _host half only does - // host arithmetic on state->scalars_h. In between, the caller calls - // fetch_device_scalars. The functions without a suffix run a whole check, - // synchronization included. - - // Copies the whole scalar struct to the host and waits for the stream: - // the one synchronization of a check. + // A check runs in two halves so that it synchronizes once: a _device half + // only queues kernels, a _host half only does arithmetic on + // state->scalars_h, and the caller calls fetch_device_scalars in between. + // The functions without a suffix run all three. + + // Copies the scalar struct to the host and waits. void fetch_device_scalars(pdhg_solver_state_t *state); // state->delta_primal_solution = to_primal - from_primal, and the same for @@ -49,14 +39,9 @@ extern "C" const double *from_dual, const double *to_dual); - // Residual vectors and their reductions from the products already in the - // state (A x in primal_product, A^T y in dual_product): the primal - // residual, the dual residual, the primal slack (the contribution of y to - // the dual objective against the finite constraint bounds given), c^T x - // and (dual slack)^T x. The optimality check passes the state's own - // vectors; feasibility polishing, whose states carry a modified problem, - // passes the original problem's objective vector, primal solution and - // finite bounds where its formulas need them. + // Residuals and objective values from A x and A^T y already in the state. + // Feasibility polishing passes vectors of the original problem where its own + // state holds a modified one. void launch_residual_check(pdhg_solver_state_t *state, const double *objective_vector, const double *primal_solution, diff --git a/src/feasibility_polish.cu b/src/feasibility_polish.cu index a023d87..7a05830 100644 --- a/src/feasibility_polish.cu +++ b/src/feasibility_polish.cu @@ -159,13 +159,9 @@ static void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state } } -// Residuals and objective values of a polishing iterate. The same reductions -// as the main termination check (termination.h) with one synchronization per -// check; the polish only needs its own side, so the other side's slots hold -// junk that is never read. The polishing states carry a modified problem, so -// the vectors the objective values need come from the original state: the -// primal polish zeroes its objective, the dual polish zeroes its primal -// iterate and its finite bounds. +// Residuals and objective values of a polishing iterate; only its own side is +// read. The objective terms use the original state's vectors, because the +// primal polish zeroes its objective and the dual polish its iterate and bounds. static void compute_polish_residual_device(pdhg_solver_state_t *state, const pdhg_solver_state_t *ori_state, bool is_primal_polish) { diff --git a/src/infeasibility.cu b/src/infeasibility.cu index 4894a41..088d283 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -14,8 +14,7 @@ See the License for the specific language governing permissions and limitations under the License. */ -// Detection of primal and dual infeasibility from ray certificates. See the -// comment above compute_infeasibility_device for the method. +// Detection of primal and dual infeasibility from ray certificates. #include "infeasibility.h" #include "reduce.h" @@ -34,8 +33,7 @@ typedef struct double primal_ray_linear_objective; double max_dual_ray_infeasibility; double dual_ray_objective; - // the objective terms with absolute values, in the same scale as the - // objective terms above; the objective margin is measured against them + // the objective terms with absolute values; the margin is measured against them double primal_objective_abs; double dual_objective_abs; } ray_candidate_stats_t; @@ -47,41 +45,28 @@ __global__ void abs_kernel(double *__restrict__ out, const double *__restrict__ out[i] = fabs(x[i]); } -// The candidate directions are stacked column by column into dense n x K and -// m x K matrices, so that A and A^T are applied to all of them with one SpMM -// each. |A| and |A^T| share the CSR indices of A and A^T and only store the -// absolute values. -// -// A check has two rounds. The first fills prod_m with A D and prod_n with -// A^T Y and reduces every scalar that depends on them. The second runs only -// when some candidate passes its objective margin: it overwrites D with |D| -// and prod_m with |A||D| (likewise on the dual side) and reduces their maxima. -// The next check rewrites all four matrices. - -// The scalars of a check are the infeas[] members of device_scalars_t -// (internal_types.h), one infeas_candidate_scalars_t per candidate. - +// The candidates are the columns of dense n x K and m x K matrices, so one SpMM +// applies A or A^T to all of them. |A| and |A^T| reuse the CSR indices of A and +// A^T. struct infeas_ctx { double *abs_val_A; double *abs_val_AT; cusparseSpMatDescr_t matA, matAT, matAbsA, matAbsAT; - // column-major, one column per candidate; the comment above says what - // each holds in the two rounds of a check + // column-major, one column per candidate; first round, then second round double *D; // n x K: primal candidates, then |D| double *Y; // m x K: dual candidates, then |Y| double *prod_m; // m x K: A D, then |A||D| double *prod_n; // n x K: A^T Y, then |A^T||Y| cusparseDnMatDescr_t dnD, dnY, dnProdM, dnProdN; void *buf_A, *buf_AT, *buf_absA, *buf_absAT; - // carried from the margin test to the finish of a check + // from the margin test, read when the best candidate is selected ray_candidate_stats_t stats[INFEAS_NUM_CANDIDATES]; bool need_primal[INFEAS_NUM_CANDIDATES]; bool need_dual[INFEAS_NUM_CANDIDATES]; }; -// The candidate directions, each x + beta * y of two device vectors (y may be -// NULL for x alone); passed to the kernels by value. +// Candidate c is x[c] + beta[c] * y[c], or x[c] alone when y[c] is NULL. typedef struct { const double *x[INFEAS_NUM_CANDIDATES]; @@ -94,25 +79,22 @@ __device__ __forceinline__ double candidate_entry(const candidate_sources_t &src return src.y[c] ? src.x[c][i] + src.beta[c] * src.y[c][i] : src.x[c][i]; } -// Violation of the sign conditions of a primal ray: (A d)_i must be >= 0 when -// the constraint has a finite lower bound and <= 0 when it has a finite upper -// bound. +// Sign conditions of a primal ray: (A d)_i >= 0 under a finite lower bound and +// <= 0 under a finite upper bound. __host__ __device__ __forceinline__ double primal_ray_row_violation(double ad, bool lb_finite, bool ub_finite) { return (lb_finite ? fmax(-ad, 0.0) : 0.0) + (ub_finite ? fmax(ad, 0.0) : 0.0); } -// Violation of the sign conditions of a dual ray on its reduced cost -// r_j = -(A^T y)_j: r_j > 0 needs a finite lower bound on x_j and r_j < 0 a -// finite upper bound. +// Sign conditions of a dual ray on r = -A^T y: r_j > 0 needs a finite lower bound +// on x_j and r_j < 0 a finite upper bound. __host__ __device__ __forceinline__ double dual_ray_reduced_cost_violation(double r, bool lb_finite, bool ub_finite) { return (lb_finite ? 0.0 : fmax(r, 0.0)) + (ub_finite ? 0.0 : fmax(-r, 0.0)); } -// The two summands of the contribution of a dual ray entry v (y_i, or the -// reduced cost r_j) to the dual objective: v^+ times the finite lower bound -// and v^- times the finite upper bound (0 when the bound is infinite). +// Dual objective terms of an entry v (y_i or r_j): v^+ times the finite lower +// bound and v^- times the finite upper bound. __host__ __device__ __forceinline__ double lower_bound_term(double v, double lower_bound_finite) { return fmax(v, 0.0) * lower_bound_finite; @@ -123,8 +105,7 @@ __host__ __device__ __forceinline__ double upper_bound_term(double v, double upp return fmin(v, 0.0) * upper_bound_finite; } -// Forms candidate blockIdx.y of a primal ray and projects it onto the sign -// conditions of a certificate. +// Forms candidate blockIdx.y and projects it onto the sign conditions. __global__ void form_primal_candidates_kernel(candidate_sources_t src, const double *__restrict__ var_lb, const double *__restrict__ var_ub, @@ -161,8 +142,6 @@ __global__ void form_dual_candidates_kernel(candidate_sources_t src, Y[(size_t)c * m + i] = v; } -// Replaces every entry of column blockIdx.y of a column-major matrix with its -// absolute value. __global__ void abs_columns_kernel(double *__restrict__ X, int rows) { int i = blockIdx.x * blockDim.x + threadIdx.x; @@ -170,9 +149,8 @@ __global__ void abs_columns_kernel(double *__restrict__ X, int rows) X[(size_t)blockIdx.y * rows + i] = fabs(X[(size_t)blockIdx.y * rows + i]); } -// First-round scalars of candidate c from its columns of D, Y and the products -// A D and A^T Y. Violations come in the rescaled space (they decide the -// certificate) and in original units (they are reported). +// First-round scalars of candidate c. Violations are reduced in the rescaled +// space, which decides the certificate, and in original units, which are reported. static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -194,9 +172,7 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) const double *constraint_lower_finite = state->constraint_lower_bound_finite_val; const double *constraint_upper_finite = state->constraint_upper_bound_finite_val; - // the primal ray d: its norm in the original space, the objective term - // c^T d with and without absolute values, and the violation of the sign - // conditions on A d + // primal ray d auto row_violation = [=] __host__ __device__(int i) { return primal_ray_row_violation(Ad[i], isfinite(constraint_lower[i]), isfinite(constraint_upper[i])); }; reduce_max( @@ -214,9 +190,7 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) [=] __host__ __device__(int i) { return row_violation(i) * constraint_rescaling[i]; }, &out->primal_ray.viol_max_orig); - // the dual ray y and its reduced cost r = -A^T y: the norm of y in the - // original space, the two parts of the dual objective term with and - // without absolute values, and the violation of the sign conditions on r + // dual ray y and its reduced cost r = -A^T y auto reduced_cost = [=] __host__ __device__(int j) { return -ATy[j]; }; auto y_lower_term = [=] __host__ __device__(int i) { return lower_bound_term(y[i], constraint_lower_finite[i]); }; auto y_upper_term = [=] __host__ __device__(int i) { return upper_bound_term(y[i], constraint_upper_finite[i]); }; @@ -371,39 +345,20 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) } // On an infeasible problem the PDHG iterates diverge along a direction whose -// primal or dual part is an infeasibility certificate (Applegate et al. 2021; -// Lu and Yang 2024 for Halpern PDHG). Three candidate directions are evaluated -// at every termination check: -// 0. the difference of iterates T(z^k) - z^k, -// 1. the iterate T(z^k), -// 2. the displacement since the last restart T(z^k) - z^0. -// -// Each candidate is projected onto the sign conditions of a certificate. Its -// primal part d certifies dual infeasibility and its dual part y primal -// infeasibility when -// the objective term has the right sign with the relative margin -// delta = infeasible_objective_margin, -// -c^T d >= delta |c|^T |d| (primal ray) -// Omega(y) >= delta sum |terms of Omega(y)| (dual ray) -// and the violation of the homogeneous constraints, relative to the -// activity the ray produces, is at most eps = eps_infeasible_relative, -// viol(A d) / max_i (|A||d|)_i <= eps (primal ray) -// viol(A^T y) / max_j (|A^T||y|)_j <= eps (dual ray). -// The two maxima of a ratio are taken in the rescaled space, where the rows -// have comparable scales; in original units a row with very large coefficients -// can dominate the activity and hide a violated row. The certificate error -// stored in the state is this ratio when the margin holds and INFINITY -// otherwise. Per side the candidate with the smallest error is kept, and its -// statistics are reported in original units. -// -// First round, kernels only: the candidates, A D and A^T Y, and every scalar -// that depends on them. infeasibility_context_create must have run. +// primal or dual part is an infeasibility certificate. Three candidates are +// tried: T(z) - z, T(z), and T(z) - z^0 (the displacement since the last +// restart). After projection onto the sign conditions, a primal ray d +// certifies dual infeasibility and a dual ray y primal infeasibility when +// -c^T d >= delta |c|^T |d| and viol(A d) / max_i (|A||d|)_i <= eps +// Omega(y) >= delta sum |Omega terms| and viol(A^T y) / max_j (|A^T||y|)_j <= eps +// with delta = infeasible_objective_margin and eps = eps_infeasible_relative. +// The ratio is taken in the rescaled space, where rows have comparable scales; +// in original units a large row can hide a violated one. void compute_infeasibility_device(pdhg_solver_state_t *state) { infeas_ctx_t *ctx = state->infeas_ctx; const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; - // 1. the candidates, one per column, projected candidate_sources_t ps, ds; // T(z) - z = reflected - pdhg, since reflected = 2 T(z) - z ps.x[0] = state->reflected_primal_solution, ps.y[0] = state->pdhg_primal_solution, ps.beta[0] = -1.0; @@ -419,23 +374,18 @@ void compute_infeasibility_device(pdhg_solver_state_t *state) form_dual_candidates_kernel<<num_blocks_dual, K), THREADS_PER_BLOCK, 0, state->stream>>>( ds, state->constraint_lower_bound, state->constraint_upper_bound, ctx->Y, m); - // 2. A D and A^T Y in one product each spmm(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, ctx->buf_A); spmm(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, ctx->buf_AT); - // 3. every scalar that depends on them, including the violations the - // certificate ratios need, so that the second round may overwrite D, Y - // and the products + // everything that needs D, Y or the products, before the second round overwrites them for (int c = 0; c < K; c++) { launch_candidate_reductions(state, c); } } -// Statistics of every candidate in the original space and the objective-side -// margin test with the relative margin objective_margin, from the fetched -// first-round scalars. Returns whether any candidate passed the test, in which -// case the second round is needed. +// Ray statistics in original units and the objective margin test. Returns +// whether any candidate passed, that is, whether the second round is needed. static bool evaluate_objective_margins(pdhg_solver_state_t *state, double objective_margin) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -465,10 +415,8 @@ static bool evaluate_objective_margins(pdhg_solver_state_t *state, double object return any; } -// Second round, kernels only: |A||d| or |A^T||y| and their maxima for each side on -// which some candidate passed the margin (all candidates of that side at -// once). D becomes |D| in place and prod_m receives |A||D|; likewise Y and -// prod_n. +// Second round, for each side on which some candidate passed the margin: D +// becomes |D| in place and prod_m becomes |A||D|; likewise on the dual side. static void launch_abs_products(pdhg_solver_state_t *state) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -510,9 +458,8 @@ static void launch_abs_products(pdhg_solver_state_t *state) } } -// Certificate errors and the best candidate per side. When no candidate of a -// side passed the margin, that side reports the statistics of candidate 0 -// with an infinite error, so the termination test never accepts it. +// Keeps the candidate with the smallest certificate error per side. A side with +// no candidate past the margin reports candidate 0 with an infinite error. static void select_best_candidates(pdhg_solver_state_t *state) { infeas_ctx_t *ctx = state->infeas_ctx; @@ -551,10 +498,8 @@ static void select_best_candidates(pdhg_solver_state_t *state) state->dual_ray_certificate_error = best_derr; } -// Turns the fetched first-round scalars into the certificate errors and ray -// statistics of the state. Host arithmetic only, unless some candidate passes -// its objective margin: then the second round runs here, with one more copy -// and synchronization. +// Runs the second round here, with one more synchronization, only when some +// candidate passes its objective margin. void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin) { if (evaluate_objective_margins(state, objective_margin)) diff --git a/src/reduce.cu b/src/reduce.cu index 7e472bd..4656788 100644 --- a/src/reduce.cu +++ b/src/reduce.cu @@ -31,9 +31,7 @@ void reduce_temp_storage_create(pdhg_solver_state_t *state) size_t sum_bytes = 0, max_bytes = 0; CUDA_CHECK(cub::DeviceReduce::Sum(NULL, sum_bytes, in, out, n_max, state->stream)); CUDA_CHECK(cub::DeviceReduce::Max(NULL, max_bytes, in, out, n_max, state->stream)); - // twice the queried size and a margin: the size depends only on the vector - // length and the value type, but reduce_check_temp_storage catches a - // transform iterator that needs more anyway + // headroom over the queried size; reduce_check_temp_storage catches the rest state->reduce_temp_bytes = 2 * (sum_bytes > max_bytes ? sum_bytes : max_bytes) + 256; CUDA_CHECK(cudaMalloc(&state->reduce_temp_d, state->reduce_temp_bytes)); } diff --git a/src/solver.cu b/src/solver.cu index 2040e26..6ffcac8 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -203,13 +203,8 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * } CUDA_CHECK(cudaGraphLaunch(graphExec, state->stream)); - // Termination check. The graph launch returns at once, so the host - // queues the kernels of the check (fixed point error, residuals and, if - // enabled, the infeasibility candidates) behind it while the GPU is - // still iterating; every scalar stays on the device until the one copy - // and synchronization of fetch_device_scalars. Only the infeasibility - // check may synchronize a second time, inside its host half, when some - // candidate passes its objective margin. + // Termination check, queued behind the graph while the GPU is still + // iterating, then one copy and synchronization. compute_fixed_point_error_device(state); compute_residual_device(state); if (detect_infeasibility) diff --git a/src/termination.cu b/src/termination.cu index 592179f..9a29c8a 100644 --- a/src/termination.cu +++ b/src/termination.cu @@ -14,9 +14,8 @@ See the License for the specific language governing permissions and limitations under the License. */ -// The termination check: fixed point error, residuals and objective values of -// the current iterate, and the criteria that decide whether the solver stops. -// See termination.h for the _device / _host protocol. +// Termination check: fixed point error, residuals, objective values, and the +// stopping criteria. #include "reduce.h" #include "termination.h" @@ -130,8 +129,6 @@ void launch_residual_check(pdhg_solver_state_t *state, state->num_constraints, state->num_variables); - // the two norms of each residual (so that the norm type is a host-side - // choice), the primal objective c^T x and the two dual objective terms device_scalars_t *out = state->scalars_d; const int m = state->num_constraints, n = state->num_variables; const double *primal_residual = state->primal_residual, *primal_slack = state->primal_slack; @@ -151,8 +148,6 @@ void launch_residual_check(pdhg_solver_state_t *state, reduce_sum(state, n, [=] __host__ __device__(int j) { return dual_slack[j] * x[j]; }, &out->dual_slack_dot); } -// Queues the residuals and objective values of the current iterate; the -// results land in state->scalars_d. void compute_residual_device(pdhg_solver_state_t *state) { cupdlpx_spmv_Ax(state->sparse_handle, state->spmv_ctx, state->pdhg_primal_solution, state->primal_product); From 64d5e5fc962a28d4a362ec2b8b5094df4b29d909 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 19:41:21 +0000 Subject: [PATCH 07/15] rename infeasible_obj_margin to infeasible_ray_objective_margin --- README.md | 2 +- include/cupdlpx_types.h | 2 +- python/README.md | 4 ++-- python/cupdlpx/PDLP.py | 2 +- python/cupdlpx/model.py | 6 +++--- python_bindings/_core_bindings.cpp | 8 ++++---- src/cli.c | 10 +++++----- src/infeasibility.cu | 2 +- src/solver.cu | 2 +- src/utils.cu | 8 ++++---- test/test_api_surface.py | 4 ++-- test/test_infeasible_unbounded.py | 6 +++--- 12 files changed, 28 insertions(+), 28 deletions(-) diff --git a/README.md b/README.md index 8d978ee..e8da529 100644 --- a/README.md +++ b/README.md @@ -118,7 +118,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | | `--no_presolve` | `flag` | Disable presolve | `enabled` | | `--eps_infeasible` | `double` | Largest relative constraint residual of an accepted ray certificate, relative to the activity `|A||x|` of the ray on the rescaled problem | `1e-8` | -| `--infeasible_obj_margin` | `double` | Smallest relative objective margin of an accepted ray certificate (`0`: plain sign test) | `1e-8` | +| `--infeasible_ray_objective_margin` | `double` | Smallest relative objective margin of an accepted ray certificate (`0`: plain sign test) | `1e-8` | | `--no_infeasibility_detection` | `flag` | Disable infeasibility detection | `enabled` | | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index a4af175..5641c1d 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -88,7 +88,7 @@ extern "C" double eps_feasible_relative; double eps_feas_polish_relative; double eps_infeasible_relative; - double infeasible_objective_margin; + double infeasible_ray_objective_margin; double time_sec_limit; int iteration_limit; } termination_criteria_t; diff --git a/python/README.md b/python/README.md index 5434aec..30a2c15 100644 --- a/python/README.md +++ b/python/README.md @@ -163,8 +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. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix, evaluated after the solver's row and column equilibration so that the scaling of the input does not decide the verdict), and its objective term has the right sign with the relative margin `InfeasibleObjMargin`. `0` disables infeasibility detection. | -| `InfeasibleObjMargin` | `infeasible_objective_margin` | float | `1e-8` | Relative margin the objective term of a ray must clear to count as a certificate: `-c^T d >= margin * |c|^T |d|` for a primal ray `d`, and the same for the dual objective term of a dual ray. It only has to clear rounding noise; `0` reduces it to a sign test. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-8` | Tolerance for declaring primal/dual infeasibility: a ray is accepted as a certificate when its constraint violation is at most this fraction of the activity `|A||x|` it produces (a relative backward error on the matrix, evaluated after the solver's row and column equilibration so that the scaling of the input does not decide the verdict), and its objective term has the right sign with the relative margin `InfeasibleRayObjectiveMargin`. `0` disables infeasibility detection. | +| `InfeasibleRayObjectiveMargin` | `infeasible_ray_objective_margin` | float | `1e-8` | Relative margin the objective term of a ray must clear to count as a certificate: `-c^T d >= margin * |c|^T |d|` for a primal ray `d`, and the same for the dual objective term of a dual ray. It only has to clear rounding noise; `0` reduces it to a sign test. | | `InfeasibilityDetection` | `infeasibility_detection` | bool | `True` | Whether to look for infeasibility certificates at all. Setting `InfeasibleTol` to `0` disables detection too. | | `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. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index dcc5a2b..27b8a63 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -42,7 +42,7 @@ "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", "InfeasibleTol": "eps_infeasible_relative", - "InfeasibleObjMargin": "infeasible_objective_margin", + "InfeasibleRayObjectiveMargin": "infeasible_ray_objective_margin", "InfeasibilityDetection": "infeasibility_detection", # scaling / step size "GeoMeanIters": "geometric_mean_iterations", diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index da60e8c..ea40bde 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -55,7 +55,7 @@ "eps_optimal_relative", "eps_feasible_relative", "eps_infeasible_relative", - "infeasible_objective_margin", + "infeasible_ray_objective_margin", "time_sec_limit", "pock_chambolle_alpha", "artificial_restart_threshold", @@ -79,9 +79,9 @@ } ) # eps_infeasible_relative = 0 also disables detection, like infeasibility_detection=False -# infeasible_objective_margin = 0 reduces the objective condition of a certificate to a sign test +# infeasible_ray_objective_margin = 0 reduces the objective condition of a certificate to a sign test _NONNEGATIVE_FLOAT_PARAMS = frozenset( - {"time_sec_limit", "matrix_zero_tol", "eps_infeasible_relative", "infeasible_objective_margin"} + {"time_sec_limit", "matrix_zero_tol", "eps_infeasible_relative", "infeasible_ray_objective_margin"} ) _STRING_PARAMS = frozenset({"optimality_norm"}) diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index df56b47..cefcdf5 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -274,7 +274,7 @@ static py::dict get_default_params_py() d["eps_optimal_relative"] = p.termination_criteria.eps_optimal_relative; d["eps_feasible_relative"] = p.termination_criteria.eps_feasible_relative; d["eps_infeasible_relative"] = p.termination_criteria.eps_infeasible_relative; - d["infeasible_objective_margin"] = p.termination_criteria.infeasible_objective_margin; + d["infeasible_ray_objective_margin"] = p.termination_criteria.infeasible_ray_objective_margin; // limits d["time_sec_limit"] = p.termination_criteria.time_sec_limit; @@ -388,7 +388,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getf("eps_optimal_relative", p->termination_criteria.eps_optimal_relative); getf("eps_feasible_relative", p->termination_criteria.eps_feasible_relative); getf("eps_infeasible_relative", p->termination_criteria.eps_infeasible_relative); - getf("infeasible_objective_margin", p->termination_criteria.infeasible_objective_margin); + getf("infeasible_ray_objective_margin", p->termination_criteria.infeasible_ray_objective_margin); // limits getf("time_sec_limit", p->termination_criteria.time_sec_limit); @@ -444,8 +444,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p if (p->termination_criteria.eps_infeasible_relative < 0.0) throw std::invalid_argument( "eps_infeasible_relative must be non-negative (0 disables infeasibility detection)."); - if (p->termination_criteria.infeasible_objective_margin < 0.0) - throw std::invalid_argument("infeasible_objective_margin must be non-negative (0 is a plain sign test)."); + if (p->termination_criteria.infeasible_ray_objective_margin < 0.0) + throw std::invalid_argument("infeasible_ray_objective_margin must be non-negative (0 is a plain sign test)."); if (p->termination_criteria.eps_feas_polish_relative <= 0.0) throw std::invalid_argument("eps_feas_polish_relative must be positive."); if (p->sv_tol <= 0.0) diff --git a/src/cli.c b/src/cli.c index 15f4960..6742c1f 100644 --- a/src/cli.c +++ b/src/cli.c @@ -216,8 +216,8 @@ void print_usage(const char *prog_name) " constraint residual of an accepted ray " "certificate (default: 1e-8; 0 also disables detection).\n"); fprintf(stderr, - " --infeasible_obj_margin Relative margin the objective term of a ray " - "must clear\n" + " --infeasible_ray_objective_margin \n" + " Relative margin the objective term of a ray must clear\n" " to count as a certificate (default: 1e-8; 0 is a " "plain sign test).\n"); fprintf(stderr, @@ -251,7 +251,7 @@ int main(int argc, char *argv[]) {"eps_feas", required_argument, 0, 1004}, {"eps_feas_polish", required_argument, 0, 1006}, {"eps_infeasible", required_argument, 0, 1019}, - {"infeasible_obj_margin", required_argument, 0, 1021}, + {"infeasible_ray_objective_margin", required_argument, 0, 1021}, {"feasibility_polishing", no_argument, 0, 'f'}, {"geo_mean_iter", required_argument, 0, 1018}, {"l_inf_ruiz_iter", required_argument, 0, 1007}, @@ -300,8 +300,8 @@ int main(int argc, char *argv[]) case 1019: // --eps_infeasible params.termination_criteria.eps_infeasible_relative = atof(optarg); break; - case 1021: // --infeasible_obj_margin - params.termination_criteria.infeasible_objective_margin = atof(optarg); + case 1021: // --infeasible_ray_objective_margin + params.termination_criteria.infeasible_ray_objective_margin = atof(optarg); break; case 'f': // --feasibility_polishing params.feasibility_polishing = true; diff --git a/src/infeasibility.cu b/src/infeasibility.cu index 088d283..c268378 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -351,7 +351,7 @@ void infeasibility_context_destroy(infeas_ctx_t *ctx) // certifies dual infeasibility and a dual ray y primal infeasibility when // -c^T d >= delta |c|^T |d| and viol(A d) / max_i (|A||d|)_i <= eps // Omega(y) >= delta sum |Omega terms| and viol(A^T y) / max_j (|A^T||y|)_j <= eps -// with delta = infeasible_objective_margin and eps = eps_infeasible_relative. +// with delta = infeasible_ray_objective_margin and eps = eps_infeasible_relative. // The ratio is taken in the rescaled space, where rows have comparable scales; // in original units a large row can hide a violated one. void compute_infeasibility_device(pdhg_solver_state_t *state) diff --git a/src/solver.cu b/src/solver.cu index 6ffcac8..3bcfbd1 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -216,7 +216,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * compute_residual_host(state, params->optimality_norm); if (detect_infeasibility) { - compute_infeasibility_host(state, params->termination_criteria.infeasible_objective_margin); + compute_infeasibility_host(state, params->termination_criteria.infeasible_ray_objective_margin); } else { diff --git a/src/utils.cu b/src/utils.cu index af5ffac..f62381b 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -277,7 +277,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->termination_criteria.iteration_limit = INT32_MAX; params->termination_criteria.eps_feas_polish_relative = 1e-6; params->termination_criteria.eps_infeasible_relative = 1e-8; - params->termination_criteria.infeasible_objective_margin = 1e-8; + params->termination_criteria.infeasible_ray_objective_margin = 1e-8; params->restart_params.artificial_restart_threshold = 0.36; params->restart_params.sufficient_reduction_for_restart = 0.2; @@ -589,9 +589,9 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro PRINT_DIFF_DBL("eps_infeasible_relative", params->termination_criteria.eps_infeasible_relative, default_params.termination_criteria.eps_infeasible_relative); - PRINT_DIFF_DBL("infeasible_objective_margin", - params->termination_criteria.infeasible_objective_margin, - default_params.termination_criteria.infeasible_objective_margin); + PRINT_DIFF_DBL("infeasible_ray_objective_margin", + params->termination_criteria.infeasible_ray_objective_margin, + default_params.termination_criteria.infeasible_ray_objective_margin); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); PRINT_DIFF_BOOL("infeasibility_detection", params->infeasibility_detection, default_params.infeasibility_detection); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); diff --git a/test/test_api_surface.py b/test/test_api_surface.py index 191ae1e..7bff31d 100644 --- a/test/test_api_surface.py +++ b/test/test_api_surface.py @@ -244,7 +244,7 @@ def test_param_value_validation(base_lp_data): with pytest.raises(ValueError): model.setParam("TimeLimit", -1.0) with pytest.raises(ValueError): - model.setParam("InfeasibleObjMargin", -1e-8) + model.setParam("InfeasibleRayObjectiveMargin", -1e-8) with pytest.raises(ValueError): model.setParam("OptimalityNorm", "l1") @@ -262,7 +262,7 @@ def test_direct_core_param_value_validation(base_lp_data): with pytest.raises(ValueError): solve_once(A, c, None, lb, ub, l, u, params={"optimality_norm": "l1"}) with pytest.raises(ValueError): - solve_once(A, c, None, lb, ub, l, u, params={"infeasible_objective_margin": -1.0}) + solve_once(A, c, None, lb, ub, l, u, params={"infeasible_ray_objective_margin": -1.0}) def test_direct_core_model_data_validation(base_lp_data): diff --git a/test/test_infeasible_unbounded.py b/test/test_infeasible_unbounded.py index aa4603c..4e2f613 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -102,9 +102,9 @@ def test_infeasibility_detection_flag_disables_detection(base_lp_data): model.optimize() assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" -def test_infeasible_objective_margin_blocks_detection(base_lp_data): +def test_infeasible_ray_objective_margin_blocks_detection(base_lp_data): """ - InfeasibleObjMargin is the relative margin the objective term of a ray must + InfeasibleRayObjectiveMargin is the relative margin the objective term of a ray must clear. A margin above 1 can never be met, because the objective term is at most the sum of the absolute values of its terms, so the infeasible LP of test_infeasible_lp runs into its iteration limit instead. @@ -113,6 +113,6 @@ def test_infeasible_objective_margin_blocks_detection(base_lp_data): l, u = l.copy(), u.copy() l[0], u[0] = 10, 10 model = Model(c, A, l, u, lb, ub) - model.setParams(OutputFlag=False, Presolve=False, InfeasibleObjMargin=2.0, IterationLimit=2000) + model.setParams(OutputFlag=False, Presolve=False, InfeasibleRayObjectiveMargin=2.0, IterationLimit=2000) model.optimize() assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" From 868c3036cb484a44434768f41293f08233dadee7 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 20:24:04 +0000 Subject: [PATCH 08/15] remove reduce_check_temp_storage --- internal/reduce.h | 16 ++++------------ src/reduce.cu | 18 +----------------- 2 files changed, 5 insertions(+), 29 deletions(-) diff --git a/internal/reduce.h b/internal/reduce.h index a21cfe9..13fe3c8 100644 --- a/internal/reduce.h +++ b/internal/reduce.h @@ -26,12 +26,10 @@ #include #include -// Temporary storage of cub::DeviceReduce, sized once and shared by every call. +// Temporary storage of cub::DeviceReduce, sized once for the longest vector and +// shared by every call. cub rejects a call that needs more. void reduce_temp_storage_create(pdhg_solver_state_t *state); -// Exits when a call needs more temporary storage than was allocated. -void reduce_check_temp_storage(const pdhg_solver_state_t *state, size_t bytes); - // out = sum over i < n of f(i). template void reduce_sum(pdhg_solver_state_t *state, int n, F f, double *out) { @@ -41,10 +39,7 @@ template void reduce_sum(pdhg_solver_state_t *state, int n, F f, do return; } auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); - size_t bytes = 0; - CUDA_CHECK(cub::DeviceReduce::Sum(NULL, bytes, it, out, n, state->stream)); - reduce_check_temp_storage(state, bytes); - bytes = state->reduce_temp_bytes; + size_t bytes = state->reduce_temp_bytes; CUDA_CHECK(cub::DeviceReduce::Sum(state->reduce_temp_d, bytes, it, out, n, state->stream)); } @@ -58,10 +53,7 @@ template void reduce_max(pdhg_solver_state_t *state, int n, F f, do return; } auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); - size_t bytes = 0; - CUDA_CHECK(cub::DeviceReduce::Max(NULL, bytes, it, out, n, state->stream)); - reduce_check_temp_storage(state, bytes); - bytes = state->reduce_temp_bytes; + size_t bytes = state->reduce_temp_bytes; CUDA_CHECK(cub::DeviceReduce::Max(state->reduce_temp_d, bytes, it, out, n, state->stream)); } diff --git a/src/reduce.cu b/src/reduce.cu index 4656788..e865a06 100644 --- a/src/reduce.cu +++ b/src/reduce.cu @@ -18,9 +18,6 @@ limitations under the License. #include "reduce.h" -#include -#include - void reduce_temp_storage_create(pdhg_solver_state_t *state) { int n_max = state->num_variables > state->num_constraints ? state->num_variables : state->num_constraints; @@ -31,19 +28,6 @@ void reduce_temp_storage_create(pdhg_solver_state_t *state) size_t sum_bytes = 0, max_bytes = 0; CUDA_CHECK(cub::DeviceReduce::Sum(NULL, sum_bytes, in, out, n_max, state->stream)); CUDA_CHECK(cub::DeviceReduce::Max(NULL, max_bytes, in, out, n_max, state->stream)); - // headroom over the queried size; reduce_check_temp_storage catches the rest - state->reduce_temp_bytes = 2 * (sum_bytes > max_bytes ? sum_bytes : max_bytes) + 256; + state->reduce_temp_bytes = sum_bytes > max_bytes ? sum_bytes : max_bytes; CUDA_CHECK(cudaMalloc(&state->reduce_temp_d, state->reduce_temp_bytes)); } - -void reduce_check_temp_storage(const pdhg_solver_state_t *state, size_t bytes) -{ - if (bytes > state->reduce_temp_bytes) - { - fprintf(stderr, - "Error: a reduction needs %zu bytes of temporary storage, %zu allocated.\n", - bytes, - state->reduce_temp_bytes); - exit(EXIT_FAILURE); - } -} From 3089d5e3b6fae0bc3f032b9cf959bcb75bd4ad89 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 21 Sep 2026 22:19:53 +0000 Subject: [PATCH 09/15] improve code readability --- internal/reduce.h | 3 + src/feasibility_polish.cu | 2 +- src/infeasibility.cu | 392 ++++++++++++++++++-------------------- src/solver.cu | 4 +- src/termination.cu | 34 +--- 5 files changed, 205 insertions(+), 230 deletions(-) diff --git a/internal/reduce.h b/internal/reduce.h index 13fe3c8..1dcf130 100644 --- a/internal/reduce.h +++ b/internal/reduce.h @@ -26,6 +26,9 @@ #include #include +// The element function of a reduction: DEVICE_FN(i) { return ...; } +#define DEVICE_FN(i) [=] __host__ __device__(int i) + // Temporary storage of cub::DeviceReduce, sized once for the longest vector and // shared by every call. cub rejects a call that needs more. void reduce_temp_storage_create(pdhg_solver_state_t *state); diff --git a/src/feasibility_polish.cu b/src/feasibility_polish.cu index 7a05830..85ef991 100644 --- a/src/feasibility_polish.cu +++ b/src/feasibility_polish.cu @@ -218,7 +218,7 @@ static void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, const double *delta = is_primal_polish ? state->delta_primal_solution : state->delta_dual_solution; int len = is_primal_polish ? state->num_variables : state->num_constraints; double *out = is_primal_polish ? &state->scalars_d->fixed_point_primal_sq : &state->scalars_d->fixed_point_dual_sq; - reduce_sum(state, len, [=] __host__ __device__(int i) { return delta[i] * delta[i]; }, out); + reduce_sum(state, len, DEVICE_FN(i) { return delta[i] * delta[i]; }, out); } static void compute_polish_fixed_point_error_host(pdhg_solver_state_t *state, bool is_primal_polish) diff --git a/src/infeasibility.cu b/src/infeasibility.cu index c268378..0c9f13f 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -26,7 +26,8 @@ limitations under the License. static const double HOST_ONE = 1.0; static const double HOST_ZERO = 0.0; -// Statistics of one candidate ray, in the original problem space. +// Statistics of one candidate ray in original units, and whether each side +// passed its objective margin. Lives for one compute_infeasibility_host call. typedef struct { double max_primal_ray_infeasibility; @@ -36,9 +37,12 @@ typedef struct // the objective terms with absolute values; the margin is measured against them double primal_objective_abs; double dual_objective_abs; + bool primal_passes_margin; + bool dual_passes_margin; } ray_candidate_stats_t; -__global__ void abs_kernel(double *__restrict__ out, const double *__restrict__ x, int n) +// out[i] = |x[i]|; out may be x. +__global__ void abs_kernel(double *out, const double *x, int n) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i < n) @@ -60,23 +64,24 @@ struct infeas_ctx double *prod_n; // n x K: A^T Y, then |A^T||Y| cusparseDnMatDescr_t dnD, dnY, dnProdM, dnProdN; void *buf_A, *buf_AT, *buf_absA, *buf_absAT; - // from the margin test, read when the best candidate is selected - ray_candidate_stats_t stats[INFEAS_NUM_CANDIDATES]; - bool need_primal[INFEAS_NUM_CANDIDATES]; - bool need_dual[INFEAS_NUM_CANDIDATES]; }; -// Candidate c is x[c] + beta[c] * y[c], or x[c] alone when y[c] is NULL. +// The direction to - from; from == NULL means to itself. +typedef struct +{ + const double *to; + const double *from; +} direction_t; + +// Wrapped in a struct so that the array is passed to a kernel by value. typedef struct { - const double *x[INFEAS_NUM_CANDIDATES]; - const double *y[INFEAS_NUM_CANDIDATES]; - double beta[INFEAS_NUM_CANDIDATES]; -} candidate_sources_t; + direction_t candidate[INFEAS_NUM_CANDIDATES]; +} candidate_directions_t; -__device__ __forceinline__ double candidate_entry(const candidate_sources_t &src, int c, int i) +__device__ __forceinline__ double direction_entry(const direction_t &d, int i) { - return src.y[c] ? src.x[c][i] + src.beta[c] * src.y[c][i] : src.x[c][i]; + return d.from ? d.to[i] - d.from[i] : d.to[i]; } // Sign conditions of a primal ray: (A d)_i >= 0 under a finite lower bound and @@ -93,20 +98,21 @@ __host__ __device__ __forceinline__ double dual_ray_reduced_cost_violation(doubl return (lb_finite ? 0.0 : fmax(r, 0.0)) + (ub_finite ? 0.0 : fmax(-r, 0.0)); } -// Dual objective terms of an entry v (y_i or r_j): v^+ times the finite lower -// bound and v^- times the finite upper bound. -__host__ __device__ __forceinline__ double lower_bound_term(double v, double lower_bound_finite) +// Dual objective term of an entry v (y_i or r_j): v^+ times the finite lower +// bound plus v^- times the finite upper bound. +__host__ __device__ __forceinline__ double bound_term(double v, double lower_finite, double upper_finite) { - return fmax(v, 0.0) * lower_bound_finite; + return fmax(v, 0.0) * lower_finite + fmin(v, 0.0) * upper_finite; } -__host__ __device__ __forceinline__ double upper_bound_term(double v, double upper_bound_finite) +// The same with each summand in absolute value. +__host__ __device__ __forceinline__ double bound_term_abs(double v, double lower_finite, double upper_finite) { - return fmin(v, 0.0) * upper_bound_finite; + return fabs(fmax(v, 0.0) * lower_finite) + fabs(fmin(v, 0.0) * upper_finite); } // Forms candidate blockIdx.y and projects it onto the sign conditions. -__global__ void form_primal_candidates_kernel(candidate_sources_t src, +__global__ void form_primal_candidates_kernel(candidate_directions_t directions, const double *__restrict__ var_lb, const double *__restrict__ var_ub, double *__restrict__ D, @@ -116,7 +122,7 @@ __global__ void form_primal_candidates_kernel(candidate_sources_t src, int c = blockIdx.y; if (j >= n) return; - double v = candidate_entry(src, c, j); + double v = direction_entry(directions.candidate[c], j); if (isfinite(var_lb[j])) v = fmax(v, 0.0); if (isfinite(var_ub[j])) @@ -124,7 +130,7 @@ __global__ void form_primal_candidates_kernel(candidate_sources_t src, D[(size_t)c * n + j] = v; } -__global__ void form_dual_candidates_kernel(candidate_sources_t src, +__global__ void form_dual_candidates_kernel(candidate_directions_t directions, const double *__restrict__ con_lb, const double *__restrict__ con_ub, double *__restrict__ Y, @@ -134,7 +140,7 @@ __global__ void form_dual_candidates_kernel(candidate_sources_t src, int c = blockIdx.y; if (i >= m) return; - double v = candidate_entry(src, c, i); + double v = direction_entry(directions.candidate[c], i); if (!isfinite(con_lb[i])) v = fmin(v, 0.0); if (!isfinite(con_ub[i])) @@ -142,13 +148,6 @@ __global__ void form_dual_candidates_kernel(candidate_sources_t src, Y[(size_t)c * m + i] = v; } -__global__ void abs_columns_kernel(double *__restrict__ X, int rows) -{ - int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < rows) - X[(size_t)blockIdx.y * rows + i] = fabs(X[(size_t)blockIdx.y * rows + i]); -} - // First-round scalars of candidate c. Violations are reduced in the rescaled // space, which decides the certificate, and in original units, which are reported. static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) @@ -173,67 +172,40 @@ static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) const double *constraint_upper_finite = state->constraint_upper_bound_finite_val; // primal ray d - auto row_violation = [=] __host__ __device__(int i) + auto d_abs_orig = DEVICE_FN(j) { return fabs(d[j]) / (variable_rescaling[j] * bound_rescaling); }; + auto objective_term = DEVICE_FN(j) { return objective[j] * d[j]; }; + auto objective_term_abs = DEVICE_FN(j) { return fabs(objective[j] * d[j]); }; + auto row_violation = DEVICE_FN(i) { return primal_ray_row_violation(Ad[i], isfinite(constraint_lower[i]), isfinite(constraint_upper[i])); }; - reduce_max( - state, - n, - [=] __host__ __device__(int j) { return fabs(d[j]) / (variable_rescaling[j] * bound_rescaling); }, - &out->primal_ray.max_abs); - reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * d[j]; }, &out->primal_ray.obj); - reduce_sum( - state, n, [=] __host__ __device__(int j) { return fabs(objective[j] * d[j]); }, &out->primal_ray.obj_abs); + auto row_violation_orig = DEVICE_FN(i) { return row_violation(i) * constraint_rescaling[i]; }; + + reduce_max(state, n, d_abs_orig, &out->primal_ray.max_abs); + reduce_sum(state, n, objective_term, &out->primal_ray.obj); + reduce_sum(state, n, objective_term_abs, &out->primal_ray.obj_abs); reduce_max(state, m, row_violation, &out->primal_ray.viol_max_scaled); - reduce_max( - state, - m, - [=] __host__ __device__(int i) { return row_violation(i) * constraint_rescaling[i]; }, - &out->primal_ray.viol_max_orig); + reduce_max(state, m, row_violation_orig, &out->primal_ray.viol_max_orig); // dual ray y and its reduced cost r = -A^T y - auto reduced_cost = [=] __host__ __device__(int j) { return -ATy[j]; }; - auto y_lower_term = [=] __host__ __device__(int i) { return lower_bound_term(y[i], constraint_lower_finite[i]); }; - auto y_upper_term = [=] __host__ __device__(int i) { return upper_bound_term(y[i], constraint_upper_finite[i]); }; - auto r_lower_term = [=] __host__ __device__(int j) - { return lower_bound_term(reduced_cost(j), variable_lower_finite[j]); }; - auto r_upper_term = [=] __host__ __device__(int j) - { return upper_bound_term(reduced_cost(j), variable_upper_finite[j]); }; - auto reduced_cost_violation = [=] __host__ __device__(int j) - { - return dual_ray_reduced_cost_violation( - reduced_cost(j), isfinite(variable_lower[j]), isfinite(variable_upper[j])); - }; - reduce_max( - state, - m, - [=] __host__ __device__(int i) { return fabs(y[i]) / (constraint_rescaling[i] * objective_rescaling); }, - &out->dual_ray.max_abs); - reduce_sum( - state, - m, - [=] __host__ __device__(int i) { return y_lower_term(i) + y_upper_term(i); }, - &out->dual_ray.obj_con_bounds); - reduce_sum( - state, - n, - [=] __host__ __device__(int j) { return r_lower_term(j) + r_upper_term(j); }, - &out->dual_ray.obj_var_bounds); - reduce_sum( - state, - m, - [=] __host__ __device__(int i) { return fabs(y_lower_term(i)) + fabs(y_upper_term(i)); }, - &out->dual_ray.obj_con_bounds_abs); - reduce_sum( - state, - n, - [=] __host__ __device__(int j) { return fabs(r_lower_term(j)) + fabs(r_upper_term(j)); }, - &out->dual_ray.obj_var_bounds_abs); + auto y_abs_orig = DEVICE_FN(i) { return fabs(y[i]) / (constraint_rescaling[i] * objective_rescaling); }; + auto con_bounds_term = DEVICE_FN(i) + { return bound_term(y[i], constraint_lower_finite[i], constraint_upper_finite[i]); }; + auto con_bounds_term_abs = DEVICE_FN(i) + { return bound_term_abs(y[i], constraint_lower_finite[i], constraint_upper_finite[i]); }; + auto var_bounds_term = + DEVICE_FN(j) { return bound_term(-ATy[j], variable_lower_finite[j], variable_upper_finite[j]); }; + auto var_bounds_term_abs = DEVICE_FN(j) + { return bound_term_abs(-ATy[j], variable_lower_finite[j], variable_upper_finite[j]); }; + auto reduced_cost_violation = DEVICE_FN(j) + { return dual_ray_reduced_cost_violation(-ATy[j], isfinite(variable_lower[j]), isfinite(variable_upper[j])); }; + auto reduced_cost_violation_orig = DEVICE_FN(j) { return reduced_cost_violation(j) * variable_rescaling[j]; }; + + reduce_max(state, m, y_abs_orig, &out->dual_ray.max_abs); + reduce_sum(state, m, con_bounds_term, &out->dual_ray.obj_con_bounds); + reduce_sum(state, n, var_bounds_term, &out->dual_ray.obj_var_bounds); + reduce_sum(state, m, con_bounds_term_abs, &out->dual_ray.obj_con_bounds_abs); + reduce_sum(state, n, var_bounds_term_abs, &out->dual_ray.obj_var_bounds_abs); reduce_max(state, n, reduced_cost_violation, &out->dual_ray.viol_max_scaled); - reduce_max( - state, - n, - [=] __host__ __device__(int j) { return reduced_cost_violation(j) * variable_rescaling[j]; }, - &out->dual_ray.viol_max_orig); + reduce_max(state, n, reduced_cost_violation_orig, &out->dual_ray.viol_max_orig); } static void spmm_buffer( @@ -270,10 +242,38 @@ spmm(cusparseHandle_t handle, cusparseSpMatDescr_t mat, cusparseDnMatDescr_t in, buffer)); } +// A CSR descriptor on the indices of M with the given values. +static cusparseSpMatDescr_t create_csr_descr(const cu_sparse_matrix_csr_t *M, double *values) +{ + cusparseSpMatDescr_t descr; + CUSPARSE_CHECK(cusparseCreateCsr(&descr, + M->num_rows, + M->num_cols, + M->num_nonzeros, + M->row_ptr, + M->col_ind, + values, + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_BASE_ZERO, + CUDA_R_64F)); + return descr; +} + +// Allocates a column-major rows x K matrix and creates its descriptor. +static double *create_dense(int rows, cusparseDnMatDescr_t *descr) +{ + const int K = INFEAS_NUM_CANDIDATES; + double *values = NULL; + CUDA_CHECK(cudaMalloc(&values, (size_t)rows * K * sizeof(double))); + CUSPARSE_CHECK(cusparseCreateDnMat(descr, rows, K, rows, values, CUDA_R_64F, CUSPARSE_ORDER_COL)); + return values; +} + infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state) { infeas_ctx_t *ctx = (infeas_ctx_t *)safe_calloc(1, sizeof(infeas_ctx_t)); - const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; + const int n = state->num_variables, m = state->num_constraints; const cu_sparse_matrix_csr_t *A = state->constraint_matrix, *AT = state->constraint_matrix_t; const int nnz = A->num_nonzeros; @@ -282,32 +282,15 @@ infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state) abs_kernel<<num_blocks_nnz, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->abs_val_A, A->val, nnz); abs_kernel<<num_blocks_nnz, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->abs_val_AT, AT->val, nnz); -#define CREATE_CSR(descr, M, vals) \ - CUSPARSE_CHECK(cusparseCreateCsr(&(descr), \ - (M)->num_rows, \ - (M)->num_cols, \ - (M)->num_nonzeros, \ - (M)->row_ptr, \ - (M)->col_ind, \ - (vals), \ - CUSPARSE_INDEX_32I, \ - CUSPARSE_INDEX_32I, \ - CUSPARSE_INDEX_BASE_ZERO, \ - CUDA_R_64F)) - CREATE_CSR(ctx->matA, A, A->val); - CREATE_CSR(ctx->matAT, AT, AT->val); - CREATE_CSR(ctx->matAbsA, A, ctx->abs_val_A); - CREATE_CSR(ctx->matAbsAT, AT, ctx->abs_val_AT); -#undef CREATE_CSR - -#define CREATE_DN(ptr, descr, rows) \ - CUDA_CHECK(cudaMalloc(&(ptr), (size_t)(rows) * K * sizeof(double))); \ - CUSPARSE_CHECK(cusparseCreateDnMat(&(descr), (rows), K, (rows), (ptr), CUDA_R_64F, CUSPARSE_ORDER_COL)) - CREATE_DN(ctx->D, ctx->dnD, n); - CREATE_DN(ctx->prod_n, ctx->dnProdN, n); - CREATE_DN(ctx->Y, ctx->dnY, m); - CREATE_DN(ctx->prod_m, ctx->dnProdM, m); -#undef CREATE_DN + ctx->matA = create_csr_descr(A, A->val); + ctx->matAT = create_csr_descr(AT, AT->val); + ctx->matAbsA = create_csr_descr(A, ctx->abs_val_A); + ctx->matAbsAT = create_csr_descr(AT, ctx->abs_val_AT); + + ctx->D = create_dense(n, &ctx->dnD); + ctx->prod_n = create_dense(n, &ctx->dnProdN); + ctx->Y = create_dense(m, &ctx->dnY); + ctx->prod_m = create_dense(m, &ctx->dnProdM); spmm_buffer(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, &ctx->buf_A); spmm_buffer(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, &ctx->buf_AT); @@ -359,20 +342,21 @@ void compute_infeasibility_device(pdhg_solver_state_t *state) infeas_ctx_t *ctx = state->infeas_ctx; const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; - candidate_sources_t ps, ds; - // T(z) - z = reflected - pdhg, since reflected = 2 T(z) - z - ps.x[0] = state->reflected_primal_solution, ps.y[0] = state->pdhg_primal_solution, ps.beta[0] = -1.0; - ds.x[0] = state->reflected_dual_solution, ds.y[0] = state->pdhg_dual_solution, ds.beta[0] = -1.0; - // T(z) - ps.x[1] = state->pdhg_primal_solution, ps.y[1] = NULL, ps.beta[1] = 0.0; - ds.x[1] = state->pdhg_dual_solution, ds.y[1] = NULL, ds.beta[1] = 0.0; - // T(z) - z^0 - ps.x[2] = state->pdhg_primal_solution, ps.y[2] = state->initial_primal_solution, ps.beta[2] = -1.0; - ds.x[2] = state->pdhg_dual_solution, ds.y[2] = state->initial_dual_solution, ds.beta[2] = -1.0; + // reflected = 2 T(z) - z, so reflected - pdhg = T(z) - z + const candidate_directions_t primal = {{ + {state->reflected_primal_solution, state->pdhg_primal_solution}, // T(z) - z + {state->pdhg_primal_solution, NULL}, // T(z) + {state->pdhg_primal_solution, state->initial_primal_solution}, // T(z) - z^0 + }}; + const candidate_directions_t dual = {{ + {state->reflected_dual_solution, state->pdhg_dual_solution}, + {state->pdhg_dual_solution, NULL}, + {state->pdhg_dual_solution, state->initial_dual_solution}, + }}; form_primal_candidates_kernel<<num_blocks_primal, K), THREADS_PER_BLOCK, 0, state->stream>>>( - ps, state->variable_lower_bound, state->variable_upper_bound, ctx->D, n); + primal, state->variable_lower_bound, state->variable_upper_bound, ctx->D, n); form_dual_candidates_kernel<<num_blocks_dual, K), THREADS_PER_BLOCK, 0, state->stream>>>( - ds, state->constraint_lower_bound, state->constraint_upper_bound, ctx->Y, m); + dual, state->constraint_lower_bound, state->constraint_upper_bound, ctx->Y, m); spmm(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, ctx->buf_A); spmm(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, ctx->buf_AT); @@ -384,128 +368,130 @@ void compute_infeasibility_device(pdhg_solver_state_t *state) } } -// Ray statistics in original units and the objective margin test. Returns -// whether any candidate passed, that is, whether the second round is needed. -static bool evaluate_objective_margins(pdhg_solver_state_t *state, double objective_margin) +// a / b, or 0 when b is 0: a ray of zero norm has no statistics. +static double divide_or_zero(double a, double b) { - infeas_ctx_t *ctx = state->infeas_ctx; - const double bc = state->constraint_bound_rescaling, oc = state->objective_vector_rescaling; + return b > 0.0 ? a / b : 0.0; +} + +// violation / product_abs_max. A zero product certifies only a zero violation. +static double certificate_error(double violation, double product_abs_max) +{ + if (product_abs_max > 0.0) + return violation / product_abs_max; + return violation > 0.0 ? INFINITY : 0.0; +} + +// Ray statistics in original units and the objective margin test of every candidate. +static void +evaluate_objective_margins(const pdhg_solver_state_t *state, double objective_margin, ray_candidate_stats_t *stats) +{ + const double bound_rescaling = state->constraint_bound_rescaling; + const double objective_rescaling = state->objective_vector_rescaling; + const double objective_scale = bound_rescaling * objective_rescaling; - bool any = false; for (int c = 0; c < INFEAS_NUM_CANDIDATES; c++) { - const infeas_candidate_scalars_t *v = &state->scalars_h->infeas[c]; - ray_candidate_stats_t *s = &ctx->stats[c]; - double primal_norm = v->primal_ray.max_abs, dual_norm = v->dual_ray.max_abs; - s->primal_ray_linear_objective = primal_norm > 0.0 ? v->primal_ray.obj / (bc * oc) / primal_norm : 0.0; - s->primal_objective_abs = primal_norm > 0.0 ? v->primal_ray.obj_abs / (bc * oc) / primal_norm : 0.0; - s->max_primal_ray_infeasibility = primal_norm > 0.0 ? v->primal_ray.viol_max_orig / bc / primal_norm : 0.0; - s->dual_ray_objective = - dual_norm > 0.0 ? (v->dual_ray.obj_con_bounds + v->dual_ray.obj_var_bounds) / (bc * oc) / dual_norm : 0.0; - s->dual_objective_abs = dual_norm > 0.0 - ? (v->dual_ray.obj_con_bounds_abs + v->dual_ray.obj_var_bounds_abs) / (bc * oc) / dual_norm - : 0.0; - s->max_dual_ray_infeasibility = dual_norm > 0.0 ? v->dual_ray.viol_max_orig / oc / dual_norm : 0.0; - ctx->need_primal[c] = s->primal_ray_linear_objective < 0.0 && s->primal_objective_abs > 0.0 && - -s->primal_ray_linear_objective >= objective_margin * s->primal_objective_abs; - ctx->need_dual[c] = s->dual_ray_objective > 0.0 && s->dual_objective_abs > 0.0 && - s->dual_ray_objective >= objective_margin * s->dual_objective_abs; - any = any || ctx->need_primal[c] || ctx->need_dual[c]; + const primal_ray_scalars_t *primal = &state->scalars_h->infeas[c].primal_ray; + const dual_ray_scalars_t *dual = &state->scalars_h->infeas[c].dual_ray; + ray_candidate_stats_t *out = &stats[c]; + + out->primal_ray_linear_objective = divide_or_zero(primal->obj / objective_scale, primal->max_abs); + out->primal_objective_abs = divide_or_zero(primal->obj_abs / objective_scale, primal->max_abs); + out->max_primal_ray_infeasibility = divide_or_zero(primal->viol_max_orig / bound_rescaling, primal->max_abs); + out->dual_ray_objective = + divide_or_zero((dual->obj_con_bounds + dual->obj_var_bounds) / objective_scale, dual->max_abs); + out->dual_objective_abs = + divide_or_zero((dual->obj_con_bounds_abs + dual->obj_var_bounds_abs) / objective_scale, dual->max_abs); + out->max_dual_ray_infeasibility = divide_or_zero(dual->viol_max_orig / objective_rescaling, dual->max_abs); + + out->primal_passes_margin = out->primal_ray_linear_objective < 0.0 && out->primal_objective_abs > 0.0 && + -out->primal_ray_linear_objective >= objective_margin * out->primal_objective_abs; + out->dual_passes_margin = out->dual_ray_objective > 0.0 && out->dual_objective_abs > 0.0 && + out->dual_ray_objective >= objective_margin * out->dual_objective_abs; } - return any; } -// Second round, for each side on which some candidate passed the margin: D -// becomes |D| in place and prod_m becomes |A||D|; likewise on the dual side. -static void launch_abs_products(pdhg_solver_state_t *state) +// Second round for the sides that need it: D becomes |D| in place and prod_m +// becomes |A||D|; likewise on the dual side. +static void launch_abs_products(pdhg_solver_state_t *state, bool primal_side, bool dual_side) { infeas_ctx_t *ctx = state->infeas_ctx; const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; - bool any_primal = false, any_dual = false; - for (int c = 0; c < K; c++) - { - any_primal = any_primal || ctx->need_primal[c]; - any_dual = any_dual || ctx->need_dual[c]; - } - if (any_primal) + if (primal_side) { - abs_columns_kernel<<num_blocks_primal, K), THREADS_PER_BLOCK, 0, state->stream>>>(ctx->D, n); + abs_kernel<<num_blocks_primal * K, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->D, ctx->D, n * K); spmm(state->sparse_handle, ctx->matAbsA, ctx->dnD, ctx->dnProdM, ctx->buf_absA); for (int c = 0; c < K; c++) { const double *p = ctx->prod_m + (size_t)c * m; - reduce_max( - state, - m, - [=] __host__ __device__(int i) { return p[i]; }, - &state->scalars_d->infeas[c].primal_ray.product_abs_max); + reduce_max(state, m, DEVICE_FN(i) { return p[i]; }, &state->scalars_d->infeas[c].primal_ray.product_abs_max); } } - if (any_dual) + if (dual_side) { - abs_columns_kernel<<num_blocks_dual, K), THREADS_PER_BLOCK, 0, state->stream>>>(ctx->Y, m); + abs_kernel<<num_blocks_dual * K, THREADS_PER_BLOCK, 0, state->stream>>>(ctx->Y, ctx->Y, m * K); spmm(state->sparse_handle, ctx->matAbsAT, ctx->dnY, ctx->dnProdN, ctx->buf_absAT); for (int c = 0; c < K; c++) { const double *p = ctx->prod_n + (size_t)c * n; - reduce_max( - state, - n, - [=] __host__ __device__(int j) { return p[j]; }, - &state->scalars_d->infeas[c].dual_ray.product_abs_max); + reduce_max(state, n, DEVICE_FN(j) { return p[j]; }, &state->scalars_d->infeas[c].dual_ray.product_abs_max); } } } // Keeps the candidate with the smallest certificate error per side. A side with // no candidate past the margin reports candidate 0 with an infinite error. -static void select_best_candidates(pdhg_solver_state_t *state) +static void select_best_candidates(pdhg_solver_state_t *state, const ray_candidate_stats_t *stats) { - infeas_ctx_t *ctx = state->infeas_ctx; - int best_p = 0, best_d = 0; - double best_perr = INFINITY, best_derr = INFINITY; + int best_primal = 0, best_dual = 0; + double best_primal_error = INFINITY, best_dual_error = INFINITY; for (int c = 0; c < INFEAS_NUM_CANDIDATES; c++) { - const infeas_candidate_scalars_t *v = &state->scalars_h->infeas[c]; - double perr = INFINITY, derr = INFINITY; - if (ctx->need_primal[c]) - { - double viol = v->primal_ray.viol_max_scaled, act = v->primal_ray.product_abs_max; - perr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); - } - if (ctx->need_dual[c]) - { - double viol = v->dual_ray.viol_max_scaled, act = v->dual_ray.product_abs_max; - derr = act > 0.0 ? viol / act : (viol > 0.0 ? INFINITY : 0.0); - } - if (perr < best_perr) + const primal_ray_scalars_t *primal = &state->scalars_h->infeas[c].primal_ray; + const dual_ray_scalars_t *dual = &state->scalars_h->infeas[c].dual_ray; + double primal_error = stats[c].primal_passes_margin + ? certificate_error(primal->viol_max_scaled, primal->product_abs_max) + : INFINITY; + double dual_error = + stats[c].dual_passes_margin ? certificate_error(dual->viol_max_scaled, dual->product_abs_max) : INFINITY; + if (primal_error < best_primal_error) { - best_perr = perr; - best_p = c; + best_primal_error = primal_error; + best_primal = c; } - if (derr < best_derr) + if (dual_error < best_dual_error) { - best_derr = derr; - best_d = c; + best_dual_error = dual_error; + best_dual = c; } } - state->max_primal_ray_infeasibility = ctx->stats[best_p].max_primal_ray_infeasibility; - state->primal_ray_linear_objective = ctx->stats[best_p].primal_ray_linear_objective; - state->primal_ray_certificate_error = best_perr; - state->max_dual_ray_infeasibility = ctx->stats[best_d].max_dual_ray_infeasibility; - state->dual_ray_objective = ctx->stats[best_d].dual_ray_objective; - state->dual_ray_certificate_error = best_derr; + state->max_primal_ray_infeasibility = stats[best_primal].max_primal_ray_infeasibility; + state->primal_ray_linear_objective = stats[best_primal].primal_ray_linear_objective; + state->primal_ray_certificate_error = best_primal_error; + state->max_dual_ray_infeasibility = stats[best_dual].max_dual_ray_infeasibility; + state->dual_ray_objective = stats[best_dual].dual_ray_objective; + state->dual_ray_certificate_error = best_dual_error; } // Runs the second round here, with one more synchronization, only when some // candidate passes its objective margin. void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin) { - if (evaluate_objective_margins(state, objective_margin)) + ray_candidate_stats_t stats[INFEAS_NUM_CANDIDATES]; + evaluate_objective_margins(state, objective_margin, stats); + + bool primal_side = false, dual_side = false; + for (int c = 0; c < INFEAS_NUM_CANDIDATES; c++) + { + primal_side = primal_side || stats[c].primal_passes_margin; + dual_side = dual_side || stats[c].dual_passes_margin; + } + if (primal_side || dual_side) { - launch_abs_products(state); + launch_abs_products(state, primal_side, dual_side); fetch_device_scalars(state); } - select_best_candidates(state); + select_best_candidates(state, stats); } diff --git a/src/solver.cu b/src/solver.cu index 3bcfbd1..11950e8 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -962,12 +962,12 @@ static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t reduce_sum( state, state->num_variables, - [=] __host__ __device__(int j) { return delta_primal[j] * delta_primal[j]; }, + DEVICE_FN(j) { return delta_primal[j] * delta_primal[j]; }, &state->scalars_d->restart_primal_sq); reduce_sum( state, state->num_constraints, - [=] __host__ __device__(int i) { return delta_dual[i] * delta_dual[i]; }, + DEVICE_FN(i) { return delta_dual[i] * delta_dual[i]; }, &state->scalars_d->restart_dual_sq); fetch_device_scalars(state); double primal_dist = sqrt(state->scalars_h->restart_primal_sq); diff --git a/src/termination.cu b/src/termination.cu index 9a29c8a..e6adadb 100644 --- a/src/termination.cu +++ b/src/termination.cu @@ -134,18 +134,13 @@ void launch_residual_check(pdhg_solver_state_t *state, const double *primal_residual = state->primal_residual, *primal_slack = state->primal_slack; const double *dual_residual = state->dual_residual, *dual_slack = state->dual_slack; const double *objective = objective_vector, *x = primal_solution; - reduce_sum( - state, - m, - [=] __host__ __device__(int i) { return primal_residual[i] * primal_residual[i]; }, - &out->primal_res_sq); - reduce_max(state, m, [=] __host__ __device__(int i) { return fabs(primal_residual[i]); }, &out->primal_res_max); - reduce_sum(state, m, [=] __host__ __device__(int i) { return primal_slack[i]; }, &out->primal_slack_sum); - reduce_sum( - state, n, [=] __host__ __device__(int j) { return dual_residual[j] * dual_residual[j]; }, &out->dual_res_sq); - reduce_max(state, n, [=] __host__ __device__(int j) { return fabs(dual_residual[j]); }, &out->dual_res_max); - reduce_sum(state, n, [=] __host__ __device__(int j) { return objective[j] * x[j]; }, &out->objective_dot); - reduce_sum(state, n, [=] __host__ __device__(int j) { return dual_slack[j] * x[j]; }, &out->dual_slack_dot); + reduce_sum(state, m, DEVICE_FN(i) { return primal_residual[i] * primal_residual[i]; }, &out->primal_res_sq); + reduce_max(state, m, DEVICE_FN(i) { return fabs(primal_residual[i]); }, &out->primal_res_max); + reduce_sum(state, m, DEVICE_FN(i) { return primal_slack[i]; }, &out->primal_slack_sum); + reduce_sum(state, n, DEVICE_FN(j) { return dual_residual[j] * dual_residual[j]; }, &out->dual_res_sq); + reduce_max(state, n, DEVICE_FN(j) { return fabs(dual_residual[j]); }, &out->dual_res_max); + reduce_sum(state, n, DEVICE_FN(j) { return objective[j] * x[j]; }, &out->objective_dot); + reduce_sum(state, n, DEVICE_FN(j) { return dual_slack[j] * x[j]; }, &out->dual_slack_dot); } void compute_residual_device(pdhg_solver_state_t *state) @@ -204,18 +199,9 @@ void compute_fixed_point_error_device(pdhg_solver_state_t *state) const int m = state->num_constraints, n = state->num_variables; const double *delta_primal = state->delta_primal_solution, *delta_dual = state->delta_dual_solution; const double *delta_dual_product = state->delta_dual_product; - reduce_sum( - state, m, [=] __host__ __device__(int i) { return delta_dual[i] * delta_dual[i]; }, &out->fixed_point_dual_sq); - reduce_sum( - state, - n, - [=] __host__ __device__(int j) { return delta_primal[j] * delta_primal[j]; }, - &out->fixed_point_primal_sq); - reduce_sum( - state, - n, - [=] __host__ __device__(int j) { return delta_dual_product[j] * delta_primal[j]; }, - &out->fixed_point_cross); + reduce_sum(state, m, DEVICE_FN(i) { return delta_dual[i] * delta_dual[i]; }, &out->fixed_point_dual_sq); + reduce_sum(state, n, DEVICE_FN(j) { return delta_primal[j] * delta_primal[j]; }, &out->fixed_point_primal_sq); + reduce_sum(state, n, DEVICE_FN(j) { return delta_dual_product[j] * delta_primal[j]; }, &out->fixed_point_cross); } void compute_fixed_point_error_host(pdhg_solver_state_t *state) From dfc663877c59208a022dad334469640540f05f9c Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Tue, 22 Sep 2026 18:11:21 +0000 Subject: [PATCH 10/15] update clang-format config --- .clang-format | 5 ++++- _typos.toml | 3 +++ 2 files changed, 7 insertions(+), 1 deletion(-) create mode 100644 _typos.toml diff --git a/.clang-format b/.clang-format index 188fb31..8c771ef 100644 --- a/.clang-format +++ b/.clang-format @@ -43,4 +43,7 @@ AttributeMacros: - __shared__ - __constant__ - __managed__ - - __launch_bounds__ \ No newline at end of file + - __launch_bounds__ + +Macros: + - "DEVICE_FN(i)=[=](int i)" diff --git a/_typos.toml b/_typos.toml new file mode 100644 index 0000000..5e0a861 --- /dev/null +++ b/_typos.toml @@ -0,0 +1,3 @@ +[default] +# benchmark instance names (netlib, MIPLIB) quoted in docs/ +extend-ignore-re = ["\\b(shs|ns|neos)[0-9-]+\\b"] From 5436a159c0a46a2abdbecf95febe9df7335273c1 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Wed, 23 Sep 2026 01:19:36 +0000 Subject: [PATCH 11/15] record the termination check into the CUDA graph --- internal/active_set_boost.h | 4 +++- internal/internal_types.h | 2 +- src/active_set_boost.cu | 32 ++++++++++--------------- src/solver.cu | 48 +++++++++++++++++++++++-------------- 4 files changed, 47 insertions(+), 39 deletions(-) diff --git a/internal/active_set_boost.h b/internal/active_set_boost.h index 2a9a938..39aa8d7 100644 --- a/internal/active_set_boost.h +++ b/internal/active_set_boost.h @@ -40,7 +40,9 @@ extern "C" void active_set_boost_free(pdhg_solver_state_t *state); - void active_set_boost_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params); + void active_set_boost_window_device(pdhg_solver_state_t *state, const pdhg_parameters_t *params); + + void active_set_boost_window_host(pdhg_solver_state_t *state); asb_action_t active_set_boost_check(pdhg_solver_state_t *state, const pdhg_parameters_t *params); diff --git a/internal/internal_types.h b/internal/internal_types.h index a56b47c..b591609 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -78,6 +78,7 @@ typedef struct double restart_primal_sq; // ||x - x_anchor||^2, the primal weight's PID input double restart_dual_sq; // ||y - y_anchor||^2, the same on the dual side infeas_candidate_scalars_t infeas[INFEAS_NUM_CANDIDATES]; + int asb_window_delta[4]; // active-set window: variables added, removed, constraints added, removed } device_scalars_t; // Private state of the infeasibility detection, defined in infeasibility.cu. @@ -239,7 +240,6 @@ typedef struct double *d_asb_dual_anchor; double *d_asb_dual_slack_anchor; double *d_asb_dual_projection_input; - int *d_asb_count; cusparseHandle_t sparse_handle; cublasHandle_t blas_handle; diff --git a/src/active_set_boost.cu b/src/active_set_boost.cu index bf5802f..311717e 100644 --- a/src/active_set_boost.cu +++ b/src/active_set_boost.cu @@ -113,7 +113,6 @@ void active_set_boost_init(pdhg_solver_state_t *state) CUDA_CHECK(cudaMalloc(&state->d_asb_dual_anchor, m * sizeof(double))); CUDA_CHECK(cudaMalloc(&state->d_asb_dual_slack_anchor, n * sizeof(double))); CUDA_CHECK(cudaMalloc(&state->d_asb_dual_projection_input, m * sizeof(double))); - CUDA_CHECK(cudaMalloc(&state->d_asb_count, 4 * sizeof(int))); CUDA_CHECK(cudaMemset(state->d_asb_var_last_free, 0xff, n * sizeof(int))); /* -1 */ CUDA_CHECK(cudaMemset(state->d_asb_row_last_binding, 0xff, m * sizeof(int))); /* -1 */ CUDA_CHECK(cudaMemset(state->d_asb_col_mask, 0, n * sizeof(bool))); @@ -143,8 +142,6 @@ void active_set_boost_free(pdhg_solver_state_t *state) cudaFree(state->d_asb_dual_slack_anchor); if (state->d_asb_dual_projection_input) cudaFree(state->d_asb_dual_projection_input); - if (state->d_asb_count) - cudaFree(state->d_asb_count); } static void asb_restore_anchor(pdhg_solver_state_t *state) @@ -172,12 +169,15 @@ static void asb_restore_anchor(pdhg_solver_state_t *state) CUDA_CHECK(cudaStreamSynchronize(state->stream)); } -static void -asb_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params, int *delta_add, int *delta_remove) +// Queues the window update; the counts arrive with the next fetch_device_scalars. +void active_set_boost_window_device(pdhg_solver_state_t *state, const pdhg_parameters_t *params) { + if (state->asb_phase == ASB_PHASE_OFF) + return; int now = state->total_count; int window_start = now - params->asb_window_iter; - CUDA_CHECK(cudaMemsetAsync(state->d_asb_count, 0, 4 * sizeof(int), state->stream)); + int *delta_count = state->scalars_d->asb_window_delta; + CUDA_CHECK(cudaMemsetAsync(delta_count, 0, 4 * sizeof(int), state->stream)); asb_var_window_kernel<<num_blocks_primal, THREADS_PER_BLOCK, 0, state->stream>>>( state->variable_lower_bound, state->variable_upper_bound, @@ -188,7 +188,7 @@ asb_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params, i state->d_asb_col_mask, now, window_start, - state->d_asb_count, + delta_count, state->num_variables); asb_row_window_kernel<<num_blocks_dual, THREADS_PER_BLOCK, 0, state->stream>>>( state->d_asb_dual_projection_input, @@ -199,15 +199,8 @@ asb_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params, i state->d_asb_row_mask, now, window_start, - state->d_asb_count, + delta_count, state->num_constraints); - int counts[4] = {0, 0, 0, 0}; - CUDA_CHECK(cudaMemcpyAsync(counts, state->d_asb_count, 4 * sizeof(int), cudaMemcpyDeviceToHost, state->stream)); - CUDA_CHECK(cudaStreamSynchronize(state->stream)); - state->asb_free_variables += counts[0] - counts[1]; - state->asb_binding_constraints += counts[2] - counts[3]; - *delta_add = counts[0] + counts[2]; - *delta_remove = counts[1] + counts[3]; } typedef enum @@ -364,13 +357,14 @@ static double asb_refresh_target(pdhg_solver_state_t *state, const pdhg_paramete return state->step_size; } -void active_set_boost_update_window(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +void active_set_boost_window_host(pdhg_solver_state_t *state) { if (state->asb_phase == ASB_PHASE_OFF) return; - int delta_add = 0, delta_remove = 0; - asb_update_window(state, params, &delta_add, &delta_remove); - state->asb_changes_since_estimate += delta_add + delta_remove; + const int *counts = state->scalars_h->asb_window_delta; + state->asb_free_variables += counts[0] - counts[1]; + state->asb_binding_constraints += counts[2] - counts[3]; + state->asb_changes_since_estimate += counts[0] + counts[1] + counts[2] + counts[3]; } asb_action_t active_set_boost_check(pdhg_solver_state_t *state, const pdhg_parameters_t *params) diff --git a/src/solver.cu b/src/solver.cu index 44eb1ad..d533074 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -154,8 +154,23 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * active_set_boost_init(state); } + bool detect_infeasibility = + params->infeasibility_detection && params->termination_criteria.eps_infeasible_relative > 0.0; + if (detect_infeasibility) + { + state->infeas_ctx = infeasibility_context_create(state); + } + + // residuals of the starting point; also runs every check kernel once before the capture display_iteration_header(params); - compute_residual(state, params->optimality_norm); + compute_fixed_point_error_device(state); + compute_residual_device(state); + if (detect_infeasibility) + { + compute_infeasibility_device(state); + } + fetch_device_scalars(state); + compute_residual_host(state, params->optimality_norm); display_iteration_stats(state, params); state->start_time = clock(); @@ -164,13 +179,6 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * cudaGraphExec_t graphExec = NULL; bool graph_created = false; - bool detect_infeasibility = - params->infeasibility_detection && params->termination_criteria.eps_infeasible_relative > 0.0; - if (detect_infeasibility) - { - state->infeas_ctx = infeasibility_context_create(state); - } - while (state->total_count < params->termination_criteria.iteration_limit) { sync_inner_count_to_gpu(state); @@ -200,6 +208,14 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * compute_next_dual_solution( state, params->termination_evaluation_frequency, params->reflection_coefficient, true); + // termination check; the copy to the host stays outside the graph + compute_fixed_point_error_device(state); + compute_residual_device(state); + if (detect_infeasibility) + { + compute_infeasibility_device(state); + } + // end CUDA graph capture cudaGraph_t graph; @@ -209,14 +225,13 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * graph_created = true; } CUDA_CHECK(cudaGraphLaunch(graphExec, state->stream)); + state->inner_count += params->termination_evaluation_frequency; + state->total_count += params->termination_evaluation_frequency; - // Termination check, queued behind the graph while the GPU is still - // iterating, then one copy and synchronization. - compute_fixed_point_error_device(state); - compute_residual_device(state); - if (detect_infeasibility) + // iteration counts are kernel arguments, so the window update stays outside the graph + if (params->active_set_boost) { - compute_infeasibility_device(state); + active_set_boost_window_device(state, params); } fetch_device_scalars(state); compute_fixed_point_error_host(state); @@ -231,12 +246,9 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * state->primal_ray_certificate_error = INFINITY; state->dual_ray_certificate_error = INFINITY; } - state->inner_count += params->termination_evaluation_frequency; - state->total_count += params->termination_evaluation_frequency; - if (params->active_set_boost) { - active_set_boost_update_window(state, params); + active_set_boost_window_host(state); } // Check Termination From fa015e9bdddffd5f9802a3316d640dd7ece884d6 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Wed, 23 Sep 2026 01:21:58 +0000 Subject: [PATCH 12/15] update build workflow --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index be4c87a..9854a55 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -90,7 +90,7 @@ jobs: run: | apt-get install -y --no-install-recommends \ cmake ninja-build zlib1g-dev \ - hipblas-dev hipsparse-dev hipcub-dev rocprim-dev + hipblas-dev hipsparse-dev hipcub-dev rocprim-dev rocthrust-dev # The ROCm apt packages don't add /opt/rocm/bin to PATH for non-interactive shells. - name: Add ROCm to PATH From dfc8e170d5eca4bb4158275053a2ab02a3b917e1 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Wed, 23 Sep 2026 03:49:43 +0000 Subject: [PATCH 13/15] create the CUDA context in a background thread --- docs/C_API.md | 2 ++ include/cupdlpx.h | 4 ++++ internal/cuda_to_hip.h | 1 + python_bindings/_core_bindings.cpp | 1 + src/cli.c | 1 + src/cupdlpx.c | 1 + src/infeasibility.cu | 21 +++++++++++++++++++++ src/utils.cu | 11 +++++++++++ 8 files changed, 42 insertions(+) diff --git a/docs/C_API.md b/docs/C_API.md index 4ae6950..18ca014 100644 --- a/docs/C_API.md +++ b/docs/C_API.md @@ -40,6 +40,8 @@ cupdlpx_result_t* solve_lp_problem( - `prob`: An LP problem built with `create_LP_problem`. The solver may clean up the matrix (e.g., drop near-zero entries), so the struct must be mutable. - `params`: Solver parameters. If `NULL`, the solver will use default parameters. +Optional: `cupdlpx_gpu_init_async()` creates the CUDA context in a background thread. `solve_lp_problem` calls it itself; calling it before building or reading the problem overlaps that work with the cold start. + #### Example: Solving a Small LP ```c #include "cupdlpx.h" diff --git a/include/cupdlpx.h b/include/cupdlpx.h index f556fc5..a5eb7d9 100644 --- a/include/cupdlpx.h +++ b/include/cupdlpx.h @@ -42,6 +42,10 @@ extern "C" // parameter void set_default_parameters(pdhg_parameters_t *params); + /* Create the CUDA context in a background thread. solve_lp_problem calls + it; calling it before reading the problem overlaps the read as well. */ + void cupdlpx_gpu_init_async(void); + /* Return 0 if params are valid. Otherwise, return nonzero and write the first error to error_message when provided. */ int cupdlpx_validate_parameters(const pdhg_parameters_t *params, char *error_message, size_t error_message_size); diff --git a/internal/cuda_to_hip.h b/internal/cuda_to_hip.h index fea1b77..4c17546 100644 --- a/internal/cuda_to_hip.h +++ b/internal/cuda_to_hip.h @@ -162,6 +162,7 @@ static inline const char *cublasGetStatusName(hipblasStatus_t status) #define cusparseDestroyDnMat hipsparseDestroyDnMat #define cusparseSpMM hipsparseSpMM #define cusparseSpMM_bufferSize hipsparseSpMM_bufferSize +#define cusparseSpMM_preprocess hipsparseSpMM_preprocess #define CUSPARSE_ORDER_COL HIPSPARSE_ORDER_COL #define CUSPARSE_SPMM_ALG_DEFAULT HIPSPARSE_SPMM_ALG_DEFAULT #define cusparseCsr2cscEx2 hipsparseCsr2cscEx2 diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index edbfd2a..158b09a 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -865,6 +865,7 @@ static py::dict read_mps_py(const std::string &filename) { // release GIL during file parsing py::gil_scoped_release release; + cupdlpx_gpu_init_async(); prob = read_mps_file(filename.c_str()); } if (!prob) diff --git a/src/cli.c b/src/cli.c index 5b2d472..ce30d08 100644 --- a/src/cli.c +++ b/src/cli.c @@ -468,6 +468,7 @@ int main(int argc, char *argv[]) return 1; } + cupdlpx_gpu_init_async(); lp_problem_t *problem = read_mps_file(filename); if (problem == NULL) diff --git a/src/cupdlpx.c b/src/cupdlpx.c index 2806b37..dff6af2 100644 --- a/src/cupdlpx.c +++ b/src/cupdlpx.c @@ -269,6 +269,7 @@ cupdlpx_result_t *solve_lp_problem(lp_problem_t *prob, const pdhg_parameters_t * fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n"); return NULL; } + cupdlpx_gpu_init_async(); // prepare parameters: use defaults if not provided pdhg_parameters_t local_params; diff --git a/src/infeasibility.cu b/src/infeasibility.cu index 0c9f13f..fd6c100 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -226,6 +226,23 @@ static void spmm_buffer( CUDA_CHECK(cudaMalloc(buffer, size > 0 ? size : 16)); } +// One-time analysis of the matrix for the repeated products. +static void spmm_preprocess( + cusparseHandle_t handle, cusparseSpMatDescr_t mat, cusparseDnMatDescr_t in, cusparseDnMatDescr_t out, void *buffer) +{ + CUSPARSE_CHECK(cusparseSpMM_preprocess(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + CUSPARSE_OPERATION_NON_TRANSPOSE, + &HOST_ONE, + mat, + in, + &HOST_ZERO, + out, + CUDA_R_64F, + CUSPARSE_SPMM_ALG_DEFAULT, + buffer)); +} + static void spmm(cusparseHandle_t handle, cusparseSpMatDescr_t mat, cusparseDnMatDescr_t in, cusparseDnMatDescr_t out, void *buffer) { @@ -296,6 +313,10 @@ infeas_ctx_t *infeasibility_context_create(pdhg_solver_state_t *state) spmm_buffer(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, &ctx->buf_AT); spmm_buffer(state->sparse_handle, ctx->matAbsA, ctx->dnD, ctx->dnProdM, &ctx->buf_absA); spmm_buffer(state->sparse_handle, ctx->matAbsAT, ctx->dnY, ctx->dnProdN, &ctx->buf_absAT); + spmm_preprocess(state->sparse_handle, ctx->matA, ctx->dnD, ctx->dnProdM, ctx->buf_A); + spmm_preprocess(state->sparse_handle, ctx->matAT, ctx->dnY, ctx->dnProdN, ctx->buf_AT); + spmm_preprocess(state->sparse_handle, ctx->matAbsA, ctx->dnD, ctx->dnProdM, ctx->buf_absA); + spmm_preprocess(state->sparse_handle, ctx->matAbsAT, ctx->dnY, ctx->dnProdN, ctx->buf_absAT); return ctx; } diff --git a/src/utils.cu b/src/utils.cu index 05100c1..ef157db 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -14,10 +14,13 @@ See the License for the specific language governing permissions and limitations under the License. */ +#include "cupdlpx.h" #include "utils.h" #include +#include #include #include +#include #ifndef CUPDLPX_VERSION #define CUPDLPX_VERSION "unknown" @@ -383,6 +386,14 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state, return reason != NULL; } +// The context is created once per process; a later CUDA call on any thread +// waits for it inside the driver. +void cupdlpx_gpu_init_async(void) +{ + static std::once_flag once; + std::call_once(once, [] { std::thread([] { cudaFree(0); }).detach(); }); +} + void set_default_parameters(pdhg_parameters_t *params) { params->geometric_mean_iterations = 12; From 375ace2c1d16837e05ed1179073bbcf1733c9436 Mon Sep 17 00:00:00 2001 From: ZedongPeng Date: Mon, 26 Oct 2026 13:08:26 -0400 Subject: [PATCH 14/15] fix the Windows build: extended lambdas need an enclosing function with external linkage --- internal/feasibility_polish.h | 2 ++ internal/infeasibility.h | 3 +++ src/feasibility_polish.cu | 3 +-- src/infeasibility.cu | 4 ++-- src/solver.cu | 4 ++-- 5 files changed, 10 insertions(+), 6 deletions(-) diff --git a/internal/feasibility_polish.h b/internal/feasibility_polish.h index 88f77c2..2a3d1ee 100644 --- a/internal/feasibility_polish.h +++ b/internal/feasibility_polish.h @@ -9,6 +9,8 @@ extern "C" void feasibility_polish(const pdhg_parameters_t *params, pdhg_solver_state_t *state); + void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish); + #ifdef __cplusplus } #endif diff --git a/internal/infeasibility.h b/internal/infeasibility.h index f17f601..5a3a4fc 100644 --- a/internal/infeasibility.h +++ b/internal/infeasibility.h @@ -37,6 +37,9 @@ extern "C" void compute_infeasibility_host(pdhg_solver_state_t *state, double objective_margin); + void launch_candidate_reductions(pdhg_solver_state_t *state, int c); + void launch_abs_products(pdhg_solver_state_t *state, bool primal_side, bool dual_side); + #ifdef __cplusplus } #endif diff --git a/src/feasibility_polish.cu b/src/feasibility_polish.cu index 7bb48cf..c59af21 100644 --- a/src/feasibility_polish.cu +++ b/src/feasibility_polish.cu @@ -34,7 +34,6 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, int k_offset, double static void perform_primal_restart(pdhg_solver_state_t *state); static void perform_dual_restart(pdhg_solver_state_t *state); -static void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish); static void compute_polish_fixed_point_error_host(pdhg_solver_state_t *state, bool is_primal_polish); static void compute_polish_fixed_point_error(pdhg_solver_state_t *state, bool is_primal_polish); static void @@ -208,7 +207,7 @@ static void compute_polish_residual_host(pdhg_solver_state_t *state, norm_type_t // The fixed point error of a polishing iterate: only its own side of the // movement term, without the square root (the restart test compares ratios). -static void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish) +void compute_polish_fixed_point_error_device(pdhg_solver_state_t *state, bool is_primal_polish) { compute_solution_difference(state, state->pdhg_primal_solution, diff --git a/src/infeasibility.cu b/src/infeasibility.cu index fd6c100..68fc80a 100644 --- a/src/infeasibility.cu +++ b/src/infeasibility.cu @@ -150,7 +150,7 @@ __global__ void form_dual_candidates_kernel(candidate_directions_t directions, // First-round scalars of candidate c. Violations are reduced in the rescaled // space, which decides the certificate, and in original units, which are reported. -static void launch_candidate_reductions(pdhg_solver_state_t *state, int c) +void launch_candidate_reductions(pdhg_solver_state_t *state, int c) { infeas_ctx_t *ctx = state->infeas_ctx; const int n = state->num_variables, m = state->num_constraints; @@ -435,7 +435,7 @@ evaluate_objective_margins(const pdhg_solver_state_t *state, double objective_ma // Second round for the sides that need it: D becomes |D| in place and prod_m // becomes |A||D|; likewise on the dual side. -static void launch_abs_products(pdhg_solver_state_t *state, bool primal_side, bool dual_side) +void launch_abs_products(pdhg_solver_state_t *state, bool primal_side, bool dual_side) { infeas_ctx_t *ctx = state->infeas_ctx; const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; diff --git a/src/solver.cu b/src/solver.cu index d533074..2bdcb59 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -53,7 +53,7 @@ __global__ void rescale_solution_kernel(double *__restrict__ primal_solution, int n_vars, int n_cons); static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, const lp_problem_t *original_problem); -static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params); +void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params); static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, const pdhg_parameters_t *params); static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_problem, const pdhg_parameters_t *params, @@ -955,7 +955,7 @@ void compute_next_dual_solution(pdhg_solver_state_t *state, } } -static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params) +void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params) { if (params->debug) { From 23215aaa790031b1ed84be0ef51d0c225618a20c Mon Sep 17 00:00:00 2001 From: ZedongPeng Date: Mon, 26 Oct 2026 13:09:42 -0400 Subject: [PATCH 15/15] map the pinned host memory calls to HIP --- internal/cuda_to_hip.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/internal/cuda_to_hip.h b/internal/cuda_to_hip.h index 4c17546..04814b7 100644 --- a/internal/cuda_to_hip.h +++ b/internal/cuda_to_hip.h @@ -47,6 +47,9 @@ limitations under the License. // Memory management #define cudaMalloc hipMalloc #define cudaFree hipFree +#define cudaHostAlloc hipHostMalloc +#define cudaHostAllocDefault hipHostMallocDefault +#define cudaFreeHost hipHostFree #define cudaMemcpy hipMemcpy #define cudaMemcpyAsync hipMemcpyAsync #define cudaMemset hipMemset