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/.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 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/_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"] diff --git a/docs/guides/c-api.md b/docs/guides/c-api.md index 583b89a..8bd9cb2 100644 --- a/docs/guides/c-api.md +++ b/docs/guides/c-api.md @@ -120,6 +120,10 @@ int main(void) { } ``` +`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. + ## Warm starts To warm-start the solver, disable presolve and provide a primal vector, diff --git a/docs/guides/command-line.md b/docs/guides/command-line.md index 3c49e27..36e3a75 100644 --- a/docs/guides/command-line.md +++ b/docs/guides/command-line.md @@ -64,8 +64,8 @@ Common architecture values include `gfx90a` for MI200, `gfx1100` for RDNA3, and `gfx1201` for RDNA4. Replace the value with the architecture supported by the installed ROCm toolchain. -On HIP builds, cuBLAS, cuSPARSE, and CUB calls are mapped to hipBLAS, -hipSPARSE, and hipCUB. +On HIP builds, cuBLAS, cuSPARSE, CUB, and Thrust calls are mapped to +hipBLAS, hipSPARSE, hipCUB, and rocThrust. ### Build outputs diff --git a/docs/implementation/gpu-backends.md b/docs/implementation/gpu-backends.md index 72d457f..5e3d907 100644 --- a/docs/implementation/gpu-backends.md +++ b/docs/implementation/gpu-backends.md @@ -12,8 +12,8 @@ The compatibility layer selects the GPU libraries. | Target | Compiler | Dense algebra | Sparse algebra | Parallel primitives | | --- | --- | --- | --- | --- | -| NVIDIA CUDA | NVCC | cuBLAS | cuSPARSE | CUB | -| AMD ROCm | `hipcc` | hipBLAS | hipSPARSE | hipCUB / rocPRIM | +| NVIDIA CUDA | NVCC | cuBLAS | cuSPARSE | CUB / Thrust | +| AMD ROCm | `hipcc` | hipBLAS | hipSPARSE | hipCUB / rocPRIM / rocThrust | With `USE_HIP=ON`, compatibility headers map CUDA calls to HIP so that the same `.cu` sources compile for AMD GPUs. @@ -62,12 +62,14 @@ cmake -B build \ - NVIDIA: [CUDA Toolkit](https://docs.nvidia.com/cuda/), [cuBLAS](https://docs.nvidia.com/cuda/cublas/), - [cuSPARSE](https://docs.nvidia.com/cuda/cusparse/), and - [CUB](https://nvidia.github.io/cccl/unstable/cub/index.html). + [cuSPARSE](https://docs.nvidia.com/cuda/cusparse/), + [CUB](https://nvidia.github.io/cccl/unstable/cub/index.html), and + [Thrust](https://nvidia.github.io/cccl/unstable/thrust/index.html). - AMD: [ROCm](https://rocm.docs.amd.com/), [hipBLAS](https://rocm.docs.amd.com/projects/hipBLAS/en/latest/), - [hipSPARSE](https://rocm.docs.amd.com/projects/hipSPARSE/en/latest/), and - [rocPRIM](https://rocm.docs.amd.com/projects/rocPRIM/en/latest/). + [hipSPARSE](https://rocm.docs.amd.com/projects/hipSPARSE/en/latest/), + [rocPRIM](https://rocm.docs.amd.com/projects/rocPRIM/en/latest/), and + [rocThrust](https://rocm.docs.amd.com/projects/rocThrust/en/latest/). - The backend dispatch described here is implemented in [`internal/cusparse_compat.h`](https://github.com/MIT-Lu-Lab/cuPDLPx/blob/main/internal/cusparse_compat.h) and diff --git a/docs/reference/parameters.md b/docs/reference/parameters.md index b20358f..153b4d9 100644 --- a/docs/reference/parameters.md +++ b/docs/reference/parameters.md @@ -155,6 +155,35 @@ name, default, type, and allowed values for the selected interface. Residuals are evaluated in the unscaled model; when [presolve](../algorithm/presolve.md) is enabled, this is the presolved model. +## Infeasibility detection + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `InfeasibilityDetection` | `true` | boolean
`false`, `true` | Look for primal and dual infeasibility certificates. | + | `InfeasibleTol` | `1e-8` | float
$[0,\infty)$ | Largest relative constraint residual of an accepted ray certificate, relative to the activity $\lvert A\rvert\,\lvert x\rvert$ of the ray; `0` disables infeasibility detection. | + | `InfeasibleRayObjectiveMargin` | `1e-8` | float
$[0,\infty)$ | Smallest relative objective margin of an accepted ray certificate: $-c^\top d \ge \delta\,\lvert c\rvert^\top\lvert d\rvert$ for a primal ray $d$ with margin $\delta$, and likewise for the dual objective term of a dual ray; `0` reduces it to a sign test. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `infeasibility_detection` | `true` | boolean
`false`, `true` | Look for primal and dual infeasibility certificates. | + | `eps_infeasible_relative` | `1e-8` | float
$[0,\infty)$ | Largest relative constraint residual of an accepted ray certificate, relative to the activity $\lvert A\rvert\,\lvert x\rvert$ of the ray; `0` disables infeasibility detection. | + | `infeasible_ray_objective_margin` | `1e-8` | float
$[0,\infty)$ | Smallest relative objective margin of an accepted ray certificate: $-c^\top d \ge \delta\,\lvert c\rvert^\top\lvert d\rvert$ for a primal ray $d$ with margin $\delta$, and likewise for the dual objective term of a dual ray; `0` reduces it to a sign test. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--no_infeasibility_detection` | Not passed | flag | Disable infeasibility detection (enabled by default). | + | `--eps_infeasible` | `1e-8` | float
$[0,\infty)$ | Largest relative constraint residual of an accepted ray certificate, relative to the activity $\lvert A\rvert\,\lvert x\rvert$ of the ray; `0` disables infeasibility detection. | + | `--infeasible_ray_objective_margin` | `1e-8` | float
$[0,\infty)$ | Smallest relative objective margin of an accepted ray certificate: $-c^\top d \ge \delta\,\lvert c\rvert^\top\lvert d\rvert$ for a primal ray $d$ with margin $\delta$, and likewise for the dual objective term of a dual ray; `0` reduces it to a sign test. | + +The ray certificate is evaluated on the rescaled problem, after row and column +equilibration, so the scaling of the input does not decide the verdict. + ## Scaling and preprocessing === "Python" 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/include/cupdlpx_types.h b/include/cupdlpx_types.h index daf03a4..044bcbf 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_ray_objective_margin; double time_sec_limit; int iteration_limit; } termination_criteria_t; @@ -109,6 +110,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/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/cuda_to_hip.h b/internal/cuda_to_hip.h index ced2bf8..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 @@ -157,6 +160,14 @@ 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 cusparseSpMM_preprocess hipsparseSpMM_preprocess +#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/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 new file mode 100644 index 0000000..5a3a4fc --- /dev/null +++ b/internal/infeasibility.h @@ -0,0 +1,45 @@ +/* +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 + + // 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); + + 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/internal/internal_types.h b/internal/internal_types.h index 04e8b44..b591609 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -25,6 +25,65 @@ limitations under the License. #include #include +// 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 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 + 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, 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 +{ + 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, 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 +{ + 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]; + 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. +typedef struct infeas_ctx infeas_ctx_t; + typedef struct { int num_rows; @@ -120,7 +179,6 @@ typedef struct double objective_vector_rescaling; double *primal_slack; double *dual_slack; - double *infeasibility_dual_scratch; /* n-vector scratch for the ray certificate; dual_slack must survive */ double rescaling_time_sec; clock_t start_time; double cumulative_time_sec; @@ -139,6 +197,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; @@ -180,14 +240,19 @@ 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; void *spmv_ctx; + infeas_ctx_t *infeas_ctx; // detection buffers, created before the main loop when detection is on + + // 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 + size_t reduce_temp_bytes; - double *ones_primal_d; - double *ones_dual_d; + 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..1dcf130 --- /dev/null +++ b/internal/reduce.h @@ -0,0 +1,63 @@ +// 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). 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 + +#include "utils.h" + +#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); + +// 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) + { + CUDA_CHECK(cudaMemsetAsync(out, 0, sizeof(double), state->stream)); + return; + } + auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); + size_t 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); 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) + { + CUDA_CHECK(cudaMemsetAsync(out, 0, sizeof(double), state->stream)); + return; + } + auto it = thrust::make_transform_iterator(thrust::counting_iterator(0), f); + size_t 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..a1dde82 --- /dev/null +++ b/internal/termination.h @@ -0,0 +1,67 @@ +/* +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 + + // 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 + // 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); + + // 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, + 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 303e4db..94a4dad 100644 --- a/internal/utils.h +++ b/internal/utils.h @@ -125,10 +125,6 @@ extern "C" const restart_parameters_t *restart_params, int termination_evaluation_frequency); - bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol); - - 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); @@ -149,12 +145,8 @@ extern "C" int get_print_frequency(int iter); - void compute_residual(pdhg_solver_state_t *state, norm_type_t optimality_norm); - void sync_step_sizes_to_gpu(pdhg_solver_state_t *state); - 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); @@ -163,27 +155,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/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 93a30b8..967f8a3 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -43,6 +43,8 @@ "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", "InfeasibleTol": "eps_infeasible_relative", + "InfeasibleRayObjectiveMargin": "infeasible_ray_objective_margin", + "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 1439b9f..27c9072 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -38,6 +38,7 @@ "bound_objective_rescaling", "feasibility_polishing", "presolve", + "infeasibility_detection", "active_set_boost", } ) @@ -57,6 +58,7 @@ "eps_optimal_relative", "eps_feasible_relative", "eps_infeasible_relative", + "infeasible_ray_objective_margin", "time_sec_limit", "pock_chambolle_alpha", "artificial_restart_threshold", diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index b4d7708..158b09a 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_ray_objective_margin"] = p.termination_criteria.infeasible_ray_objective_margin; // limits d["time_sec_limit"] = p.termination_criteria.time_sec_limit; @@ -297,6 +298,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 @@ -421,6 +423,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_ray_objective_margin", p->termination_criteria.infeasible_ray_objective_margin); // limits getf("time_sec_limit", p->termination_criteria.time_sec_limit); @@ -444,6 +447,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 @@ -861,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/active_set_boost.cu b/src/active_set_boost.cu index c5c810e..311717e 100644 --- a/src/active_set_boost.cu +++ b/src/active_set_boost.cu @@ -15,6 +15,7 @@ limitations under the License. */ #include "active_set_boost.h" +#include "termination.h" #include "utils.h" #include #include @@ -112,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))); @@ -142,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) @@ -171,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, @@ -187,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, @@ -198,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 @@ -274,6 +268,8 @@ static asb_action_t asb_revert(pdhg_solver_state_t *state, const pdhg_parameters state->max_dual_ray_infeasibility = 0.0; state->primal_ray_linear_objective = 0.0; state->dual_ray_objective = 0.0; + state->primal_ray_certificate_error = INFINITY; + state->dual_ray_certificate_error = INFINITY; state->step_size = state->base_step_size; state->primal_weight = state->asb_anchor_primal_weight; state->primal_weight_error_sum = state->asb_anchor_pw_error_sum; @@ -361,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/cli.c b/src/cli.c index cc85c93..ce30d08 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, "ASB Step Raise Count: %d\n", result->asb_raise_count); @@ -212,9 +216,22 @@ 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, + " --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, " --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"); @@ -272,6 +289,8 @@ 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, 1031}, + {"infeasible_ray_objective_margin", required_argument, 0, 1033}, {"feasibility_polishing", no_argument, 0, 'f'}, {"geo_mean_iter", required_argument, 0, 1018}, {"l_inf_ruiz_iter", required_argument, 0, 1007}, @@ -282,6 +301,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, 1032}, {"no_presolve", no_argument, 0, 1015}, {"matrix_zero_tol", required_argument, 0, 1016}, {"infinite_bound", required_argument, 0, 1017}, @@ -328,6 +348,12 @@ int main(int argc, char *argv[]) case 1006: // --eps_feas_polish_relative params.termination_criteria.eps_feas_polish_relative = atof(optarg); break; + case 1031: // --eps_infeasible + params.termination_criteria.eps_infeasible_relative = atof(optarg); + break; + case 1033: // --infeasible_ray_objective_margin + params.termination_criteria.infeasible_ray_objective_margin = atof(optarg); + break; case 'f': // --feasibility_polishing params.feasibility_polishing = true; break; @@ -373,6 +399,9 @@ int main(int argc, char *argv[]) } } break; + case 1032: // --no_infeasibility_detection + params.infeasibility_detection = false; + break; case 1015: // --no_presolve params.presolve = false; break; @@ -439,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 10622e9..dff6af2 100644 --- a/src/cupdlpx.c +++ b/src/cupdlpx.c @@ -215,9 +215,12 @@ int cupdlpx_validate_parameters(const pdhg_parameters_t *p, char *error_message, CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_feasible_relative > 0.0, "eps_feasible_relative must be positive (got %g)", p->termination_criteria.eps_feasible_relative); - CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_infeasible_relative > 0.0, - "eps_infeasible_relative must be positive (got %g)", + CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_infeasible_relative >= 0.0, + "eps_infeasible_relative must be nonnegative, 0 disables infeasibility detection (got %g)", p->termination_criteria.eps_infeasible_relative); + CUPDLPX_CHECK_PARAM(p->termination_criteria.infeasible_ray_objective_margin >= 0.0, + "infeasible_ray_objective_margin must be nonnegative, 0 is a plain sign test (got %g)", + p->termination_criteria.infeasible_ray_objective_margin); CUPDLPX_CHECK_PARAM(p->termination_criteria.eps_feas_polish_relative > 0.0, "eps_feas_polish_relative must be positive (got %g)", p->termination_criteria.eps_feas_polish_relative); @@ -266,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/feasibility_polish.cu b/src/feasibility_polish.cu index 544da5c..c59af21 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,12 @@ 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_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 +51,189 @@ 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; 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) +{ + 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). +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, 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) +{ + 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 +320,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 +350,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; @@ -211,7 +395,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; } @@ -241,8 +425,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; @@ -557,50 +744,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..68fc80a --- /dev/null +++ b/src/infeasibility.cu @@ -0,0 +1,518 @@ +/* +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. + +#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 original units, and whether each side +// passed its objective margin. Lives for one compute_infeasibility_host call. +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; 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; + +// 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) + out[i] = fabs(x[i]); +} + +// 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; 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; +}; + +// 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 +{ + direction_t candidate[INFEAS_NUM_CANDIDATES]; +} candidate_directions_t; + +__device__ __forceinline__ double direction_entry(const direction_t &d, int 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 +// <= 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); +} + +// 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)); +} + +// 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_finite + fmin(v, 0.0) * upper_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 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_directions_t directions, + 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 = direction_entry(directions.candidate[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_directions_t directions, + 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 = direction_entry(directions.candidate[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; +} + +// First-round scalars of candidate c. Violations are reduced in the rescaled +// space, which decides the certificate, and in original units, which are reported. +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; + + // primal ray d + 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])); }; + 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, row_violation_orig, &out->primal_ray.viol_max_orig); + + // dual ray y and its reduced cost r = -A^T y + 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, reduced_cost_violation_orig, &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)); +} + +// 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) +{ + 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)); +} + +// 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; + 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); + + 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); + 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; +} + +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); +} + +// On an infeasible problem the PDHG iterates diverge along a direction whose +// 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_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) +{ + infeas_ctx_t *ctx = state->infeas_ctx; + const int n = state->num_variables, m = state->num_constraints, K = INFEAS_NUM_CANDIDATES; + + // 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>>>( + 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>>>( + 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); + + // 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); + } +} + +// 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) +{ + 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; + + for (int c = 0; c < INFEAS_NUM_CANDIDATES; 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; + } +} + +// Second round for the sides that need it: D becomes |D| in place and prod_m +// becomes |A||D|; likewise on the 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; + + if (primal_side) + { + 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, DEVICE_FN(i) { return p[i]; }, &state->scalars_d->infeas[c].primal_ray.product_abs_max); + } + } + if (dual_side) + { + 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, 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, const ray_candidate_stats_t *stats) +{ + 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 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_primal_error = primal_error; + best_primal = c; + } + if (dual_error < best_dual_error) + { + best_dual_error = dual_error; + best_dual = c; + } + } + 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) +{ + 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, primal_side, dual_side); + fetch_device_scalars(state); + } + select_best_candidates(state, stats); +} diff --git a/src/reduce.cu b/src/reduce.cu new file mode 100644 index 0000000..e865a06 --- /dev/null +++ b/src/reduce.cu @@ -0,0 +1,33 @@ +/* +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" + +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)); + state->reduce_temp_bytes = sum_bytes > max_bytes ? sum_bytes : max_bytes; + CUDA_CHECK(cudaMalloc(&state->reduce_temp_d, state->reduce_temp_bytes)); +} diff --git a/src/solver.cu b/src/solver.cu index 8317eb7..2bdcb59 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -17,10 +17,13 @@ limitations under the License. #include "active_set_boost.h" #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 @@ -49,21 +52,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); +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); @@ -160,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(); @@ -198,6 +207,15 @@ 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); + + // 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; @@ -207,16 +225,30 @@ 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); state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; + // iteration counts are kernel arguments, so the window update stays outside the graph if (params->active_set_boost) { - active_set_boost_update_window(state, params); + active_set_boost_window_device(state, params); + } + fetch_device_scalars(state); + compute_fixed_point_error_host(state); + compute_residual_host(state, params->optimality_norm); + if (detect_infeasibility) + { + compute_infeasibility_host(state, params->termination_criteria.infeasible_ray_objective_margin); + } + else + { + // detection off: no candidate evaluation, no |A| copies + state->primal_ray_certificate_error = INFINITY; + state->dual_ray_certificate_error = INFINITY; + } + if (params->active_set_boost) + { + active_set_boost_window_host(state); } // Check Termination @@ -465,9 +497,9 @@ static pdhg_solver_state_t *initialize_solver_state(const lp_problem_t *working_ ALLOC_ZERO(state->reflected_primal_solution, var_bytes); ALLOC_ZERO(state->dual_product, var_bytes); ALLOC_ZERO(state->dual_slack, var_bytes); - ALLOC_ZERO(state->infeasibility_dual_scratch, 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); @@ -617,26 +649,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, @@ -646,7 +666,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); return state; } @@ -849,27 +868,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, @@ -957,28 +955,33 @@ 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) { printf("[restart] iter %d: %s\n", state->total_count, state->last_restart_reason); } state->restart_count++; - 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)); + DEVICE_FN(j) { return delta_primal[j] * delta_primal[j]; }, + &state->scalars_d->restart_primal_sq); + reduce_sum( + state, + state->num_constraints, + 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); + double dual_dist = sqrt(state->scalars_h->restart_dual_sq); double ratio_infeas = state->relative_dual_residual / state->relative_primal_residual; @@ -1008,22 +1011,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; @@ -1076,47 +1083,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)); - // measure in the M-norm of the smaller (global) step: M is indefinite for a - // controller-boosted step - double norm_step = fmin(state->step_size, state->base_step_size); - interaction = 2 * norm_step * cross_term; - - state->fixed_point_error = sqrt(movement + interaction); -} - void pdhg_solver_state_free(pdhg_solver_state_t *state) { if (state == NULL) @@ -1130,6 +1096,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) @@ -1203,8 +1171,6 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state) CUDA_CHECK(cudaFree(state->primal_slack)); if (state->dual_slack) CUDA_CHECK(cudaFree(state->dual_slack)); - if (state->infeasibility_dual_scratch) - CUDA_CHECK(cudaFree(state->infeasibility_dual_scratch)); if (state->primal_residual) CUDA_CHECK(cudaFree(state->primal_residual)); if (state->dual_residual) @@ -1213,10 +1179,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..7125fb1 --- /dev/null +++ b/src/termination.cu @@ -0,0 +1,269 @@ +/* +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. +*/ + +// Termination check: fixed point error, residuals, objective values, and the +// stopping criteria. + +#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); + + 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, 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) +{ + 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, 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) +{ + 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; + // measure in the M-norm of the smaller (global) step: M is indefinite for a + // controller-boosted step + double norm_step = fmin(state->step_size, state->base_step_size); + double interaction = 2.0 * norm_step * 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 413797f..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" @@ -347,60 +350,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) @@ -437,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; @@ -457,7 +414,8 @@ 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->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; @@ -469,6 +427,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; params->active_set_boost = true; @@ -796,7 +755,11 @@ 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_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_BOOL("debug", params->debug, default_params.debug); 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); @@ -951,341 +914,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->infeasibility_dual_scratch, - 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->infeasibility_dual_scratch); - 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->infeasibility_dual_scratch, - state->variable_rescaling); - - state->max_primal_ray_infeasibility = - get_vector_inf_norm(state->blas_handle, state->num_constraints, state->primal_slack); - double dual_scratch_norm = - get_vector_inf_norm(state->blas_handle, state->num_variables, state->infeasibility_dual_scratch); - state->max_dual_ray_infeasibility = dual_scratch_norm; - - double scaling_factor = fmax(dual_ray_inf_norm, dual_scratch_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) { @@ -1448,241 +1076,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_api_surface.py b/test/test_api_surface.py index 29ec2ce..5381ef5 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("InfeasibleRayObjectiveMargin", -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_ray_objective_margin": -1.0}) with pytest.raises(ValueError, match="Unknown parameter"): solve_once(A, c, None, lb, ub, l, u, params={"eps_optimal_relatve": 1e-6}) with pytest.raises(ValueError, match="must be a number"): 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..4e2f613 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -74,4 +74,45 @@ 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}" + +def test_infeasible_ray_objective_margin_blocks_detection(base_lp_data): + """ + 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. + """ + 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, InfeasibleRayObjectiveMargin=2.0, IterationLimit=2000) + model.optimize() + assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}"