Primal infeasibility result with C but not with Python interface

Hi,

I recently encountered a case where the OSQP solver gives an incorrect primal infeasible result for a small QP problem. While trying to determine the cause of the issue, I recreated the problem setup with the python interface. To my surprise I found that despite identical input data, settings, and OSQP version, I did not have the same issue with the python interface.

I created some minimal standalone examples in C and Python. Links to code and output below.

C code: main.c · GitHub
C output:

-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 155, constraints m = 262
          nnz(P) + nnz(A) = 599
settings: linear system solver = qdldl,
          eps_abs = 1.0e-02, eps_rel = 1.0e-04,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 5),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, 
iter   objective    pri res    dua res    rho
   1   1.3099e+06   2.80e+01   1.68e+06   1.00e-01
  55   1.0000e+30   2.50e-01   2.59e-04   8.68e-03

status:               primal infeasible
number of iterations: 55
optimal rho estimate: 3.28e-01

Python code: main.py · GitHub
Python output:

-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 155, constraints m = 262
          nnz(P) + nnz(A) = 599
settings: linear system solver = qdldl,
          eps_abs = 1.0e-02, eps_rel = 1.0e-04,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 5),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, time_limit: off

iter   objective    pri res    dua res    rho        time
   1   1.3099e+06   2.80e+01   1.68e+06   1.00e-01   3.57e-04s
 115   1.1925e+06   4.69e-02   3.01e-02   1.74e-03   1.45e-03s

status:               solved
number of iterations: 115
optimal objective:    1192485.8137
run time:             1.56e-03s
optimal rho estimate: 2.34e-03

The examples above use qdldl, but the result is the same for me with mkl pardiso.
Any help with diagnosing the issue would be appreciated.

Thanks,
Paul

Some discussion on github: Incorrect primal infeasibility result on a small QP · Issue #346 · osqp/osqp · GitHub