Improve infeasibility detection

.clang-format

@@ -43,4 +43,7 @@ AttributeMacros:
   - __shared__
   - __constant__
   - __managed__
-  - __launch_bounds__
+  - __launch_bounds__
+
+Macros:
+  - "DEVICE_FN(i)=[=](int i)"

.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

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 $<$<COMPILE_LANGUAGE:CUDA>:--extended-lambda>)
+endif()
+
 if(MSVC)
     target_compile_options(cupdlpx_compile_flags INTERFACE
         $<$<OR:$<COMPILE_LANGUAGE:C>,$<COMPILE_LANGUAGE:CXX>>:/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
   )

_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"]

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,

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
 

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

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<br>`false`, `true` | Look for primal and dual infeasibility certificates. |
+    | `InfeasibleTol` | `1e-8` | float<br>$[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<br>$[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<br>`false`, `true` | Look for primal and dual infeasibility certificates. |
+    | `eps_infeasible_relative` | `1e-8` | float<br>$[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<br>$[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<br>$[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<br>$[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"

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);

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;

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);
 

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
 

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

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

internal/internal_types.h

@@ -25,6 +25,65 @@ limitations under the License.
 #include <stdbool.h>
 #include <time.h>
 
+// Scalars of a termination check, one reduction each, in the rescaled space
+// except viol_max_orig. Name parts: res = residual, sq = sum of squares,
+// max = infinity norm, dot = inner product.
+
+#define INFEAS_NUM_CANDIDATES 3 // candidate rays of one infeasibility check
+
+// One 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;