Add a pure C interface

CMakeLists.txt

@@ -18,8 +18,10 @@ else()
 endif()
 option(ELASTIQP_BUILD_TESTS "Build the elastiqp test suite" ${_elastiqp_top_level})
 option(ELASTIQP_BUILD_PYTHON "Build the Python bindings and JAX FFI" ${_elastiqp_top_level})
+option(ELASTIQP_BUILD_C "Build the C interface library (elastiqp_c)" ${_elastiqp_top_level})
 if(SKBUILD)  # pip wheel build: bindings only, no install/export
   set(ELASTIQP_BUILD_TESTS OFF)
+  set(ELASTIQP_BUILD_C OFF)
   set(ELASTIQP_BUILD_PYTHON ON)
 endif()
 
@@ -34,6 +36,27 @@ target_include_directories(elastiqp INTERFACE
 target_link_libraries(elastiqp INTERFACE Eigen3::Eigen)
 target_compile_features(elastiqp INTERFACE cxx_std_17)
 
+# C interface library
+if(ELASTIQP_BUILD_C)
+  enable_language(C)
+  add_library(elastiqp_c src/elastiqp_c.cc)
+  add_library(elastiqp::elastiqp_c ALIAS elastiqp_c)
+  target_link_libraries(elastiqp_c PRIVATE elastiqp::elastiqp)
+  target_include_directories(elastiqp_c PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:include>)
+  target_compile_definitions(elastiqp_c PRIVATE ELASTIQP_C_BUILDING)
+  if(BUILD_SHARED_LIBS)
+    target_compile_definitions(elastiqp_c PUBLIC ELASTIQP_C_SHARED)
+  endif()
+  set_target_properties(elastiqp_c PROPERTIES
+    CXX_VISIBILITY_PRESET hidden
+    VISIBILITY_INLINES_HIDDEN ON
+    POSITION_INDEPENDENT_CODE ON
+    VERSION ${PROJECT_VERSION}
+    SOVERSION ${PROJECT_VERSION_MAJOR})
+endif()
+
 # Headers for tests (never installed)
 add_library(elastiqp_testing INTERFACE)
 add_library(elastiqp::testing ALIAS elastiqp_testing)
@@ -75,6 +98,19 @@ if(ELASTIQP_BUILD_TESTS)
 
   elastiqp_test(test_lp)
   add_test(NAME elastiqp.lp COMMAND test_lp)
+
+  if(ELASTIQP_BUILD_C)
+    add_executable(test_c_api tests/test_c_api.c)
+    target_link_libraries(test_c_api PRIVATE elastiqp::elastiqp_c)
+    if(NOT MSVC)
+      target_link_libraries(test_c_api PRIVATE m)
+    endif()
+    add_test(NAME elastiqp.c_api COMMAND test_c_api)
+
+    elastiqp_test(test_c_parity)
+    target_link_libraries(test_c_parity PRIVATE elastiqp::elastiqp_c)
+    add_test(NAME elastiqp.c_parity COMMAND test_c_parity)
+  endif()
 endif()
 
 # Python bindings and JAX FFI

README.md

@@ -94,7 +94,7 @@ Note: if using UV, you can directly replace the above `pip` commands with `uv pi
 elastiqp::Solution sol = elastiqp::Solve(Q, q, A, b, G, h, penalty);
 // or without equalities: elastiqp::Solve(Q, q, G, h, penalty)
 
-// If you are solving multiple times in a control loop:
+// Repeated solves
 elastiqp::Solver solver;
 solver.setup(Q, q, A, b, G, h, penalty);
 while (running) {
@@ -104,6 +104,32 @@ while (running) {
 }
 ```
 
+### C
+
+```c
+#include "elastiqp/elastiqp_c.h"
+
+// One-shot solve; penalty >= ELASTIQP_INF makes a row hard.
+ElastiQPProblem qp = {.n = n, .m = m, .p = p, .Q = Q, .q = q, .A = A, .b = b,
+                      .G = G, .h = h, .penalty = penalty};
+ElastiQPSolution res = {.x = x, .z = z};  // caller-owned outputs; NULL ones are skipped
+elastiqp_solve_once(&res, &qp, NULL);     // NULL: default settings
+
+// Repeated solves
+ElastiQPWorkspace* work;
+elastiqp_setup(&work, &qp, NULL);
+while (running) {
+  elastiqp_set_q(work, q_k); elastiqp_set_h(work, h_k); elastiqp_set_b(work, b_k);
+  elastiqp_set_G(work, G_k);  // only rows that differ are refreshed
+  if (elastiqp_solve(&res, work) != ELASTIQP_SOLVED) { /* ... */ }
+}
+elastiqp_free(work);
+```
+
+Pass settings as an `ElastiQPSettings` struct or set them by name: `elastiqp_set_option(work, "eps_abs", 1e-8)`. Use `elastiqp_num_options` and `elastiqp_option_info` to list names, types, defaults, ranges, and descriptions.
+
+With CMake: `find_package(elastiqp)` and link `elastiqp::elastiqp_c`.
+
 ### Python
 
 ```python

cmake/elastiqpConfig.cmake.in

@@ -5,4 +5,14 @@ find_dependency(Eigen3)
 
 include("${CMAKE_CURRENT_LIST_DIR}/elastiqpTargets.cmake")
 
+# The static C interface requires the C++ linker and runtime.
+if(TARGET elastiqp::elastiqp_c)
+  get_target_property(_elastiqp_c_type elastiqp::elastiqp_c TYPE)
+  get_property(_elastiqp_languages GLOBAL PROPERTY ENABLED_LANGUAGES)
+  if(_elastiqp_c_type STREQUAL "STATIC_LIBRARY"
+     AND NOT "CXX" IN_LIST _elastiqp_languages)
+    enable_language(CXX)
+  endif()
+endif()
+
 check_required_components(elastiqp)

cmake/elastiqp_install.cmake

@@ -3,7 +3,14 @@ include(CMakePackageConfigHelpers)
 
 install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
 
-install(TARGETS elastiqp EXPORT elastiqpTargets)
+set(_elastiqp_install_targets elastiqp)
+if(ELASTIQP_BUILD_C)
+  list(APPEND _elastiqp_install_targets elastiqp_c)
+endif()
+install(TARGETS ${_elastiqp_install_targets} EXPORT elastiqpTargets
+  ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+  LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+  RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
 install(EXPORT elastiqpTargets
   NAMESPACE elastiqp::
   DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/elastiqp)

include/elastiqp/common.hpp

@@ -12,6 +12,9 @@ namespace elastiqp {
 using Eigen::MatrixXd;
 using Eigen::VectorXd;
 
+using MatRef = Eigen::Ref<const MatrixXd>;
+using VecRef = Eigen::Ref<const VectorXd>;
+
 enum class Status {
   kUnsolved = 0,
   kSolved = 1,
@@ -103,7 +106,7 @@ inline double ruiz_cost_gamma(const MatrixXd& Q) {
 // the reachable residual; nonzero only when A is rank-deficient.
 class EqCertificate {
  public:
-  void set_A(const MatrixXd& A) {
+  void set_A(const MatRef& A) {
     m_ = A.rows();
     n_ = A.cols();
     rank_deficient_ = false;
@@ -127,7 +130,7 @@ class EqCertificate {
   }
   bool rank_deficient() const { return rank_deficient_; }
 
-  double bound(const VectorXd& b) const {
+  double bound(const VecRef& b) const {
     if (!rank_deficient_ || m_ == 0) return 0.0;
     const VectorXd xls = qr_.solve(b);
     return (A_ * xls - b).norm() / std::sqrt(static_cast<double>(m_));

include/elastiqp/das.hpp

@@ -76,9 +76,9 @@ class Solver {
     kDropped
   };
 
-  void setup(const MatrixXd& Q, const VectorXd& q, const MatrixXd& A,
-             const VectorXd& b, const MatrixXd& G, const VectorXd& h,
-             const VectorXd& penalty) {
+  void setup(const MatRef& Q, const VecRef& q, const MatRef& A,
+             const VecRef& b, const MatRef& G, const VecRef& h,
+             const VecRef& penalty) {
     n_ = static_cast<int>(q.size());
     m_ = static_cast<int>(b.size());
     p_ = static_cast<int>(h.size());
@@ -158,43 +158,43 @@ class Solver {
     sol_.status = Status::kUnsolved;
     rw_.ready = false;
   }
-  void setup(const MatrixXd& Q, const VectorXd& q, const MatrixXd& A,
-             const VectorXd& b, const MatrixXd& G, const VectorXd& h,
+  void setup(const MatRef& Q, const VecRef& q, const MatRef& A,
+             const VecRef& b, const MatRef& G, const VecRef& h,
              double penalty) {
     setup(Q, q, A, b, G, h, VectorXd::Constant(h.size(), penalty));
   }
-  void setup(const MatrixXd& Q, const VectorXd& q, const MatrixXd& G,
-             const VectorXd& h, const VectorXd& penalty) {
+  void setup(const MatRef& Q, const VecRef& q, const MatRef& G,
+             const VecRef& h, const VecRef& penalty) {
     setup(Q, q, MatrixXd(0, q.size()), VectorXd(0), G, h, penalty);
   }
-  void setup(const MatrixXd& Q, const VectorXd& q, const MatrixXd& G,
-             const VectorXd& h, double penalty) {
+  void setup(const MatRef& Q, const VecRef& q, const MatRef& G,
+             const VecRef& h, double penalty) {
     setup(Q, q, G, h, VectorXd::Constant(h.size(), penalty));
   }
 
   int n() const { return n_; }
   int m() const { return m_; }
   int p() const { return p_; }
 
-  void set_Q(const MatrixXd& Q) {
+  void set_Q(const MatRef& Q) {
     Q_ = 0.5 * (Q + Q.transpose());
     Q_dirty_ = true;
   }
-  void set_q(const VectorXd& q) {
+  void set_q(const VecRef& q) {
     q_ = q;
     rhs_dirty_ = true;
   }
-  void set_A(const MatrixXd& A) { set_rows(A, 0); }
-  void set_b(const VectorXd& b) {
+  void set_A(const MatRef& A) { set_rows(A, 0); }
+  void set_b(const VecRef& b) {
     rhs_.head(m_) = b;
     rhs_dirty_ = true;
   }
-  void set_G(const MatrixXd& G) { set_rows(G, m_); }
-  void set_h(const VectorXd& h) {
+  void set_G(const MatRef& G) { set_rows(G, m_); }
+  void set_h(const VecRef& h) {
     rhs_.tail(p_) = h;
     rhs_dirty_ = true;
   }
-  void set_penalty(const VectorXd& penalty) {
+  void set_penalty(const VecRef& penalty) {
     penalty_ = penalty;
     penalty_dirty_ = true;
   }
@@ -203,13 +203,24 @@ class Solver {
 
   // Seeds states and multipliers from (x, y, z) on the next solve(), replacing
   // the internal warm start.
-  void set_warm_start(const VectorXd& x, const VectorXd& y, const VectorXd& z) {
+  void set_warm_start(const VecRef& x, const VecRef& y, const VecRef& z) {
     warm_x_ = x;
     warm_y_ = y;
     warm_z_ = z;
     explicit_warm_ = true;
   }
 
+  // Can relax/differentiate if the last solve succeeded and no data changed
+  bool can_relax() const {
+    if (!have_solution_ || Q_dirty_ || penalty_dirty_ || rhs_dirty_) {
+      return false;
+    }
+    for (char c : col_dirty_) {
+      if (c != 0) return false;
+    }
+    return true;
+  }
+
   RowState row_state(int i) const {
     return state_[static_cast<size_t>(m_ + i)];
   }
@@ -346,13 +357,13 @@ class Solver {
   // start untouched. The first call after setup() allocates its workspace.
   const Solution& relax(double kappa, double tol = 1e-6, int max_iter = 50) {
     if (p_ == 0 || kappa <= 0.0) return sol_;
-    bool stale = Q_dirty_ || penalty_dirty_ || rhs_dirty_;
-    for (char c : col_dirty_) stale |= c != 0;
-    if (stale && have_solution_) {
-      sol_.status = Status::kUnsolved;
-      sol_.converged = 0;
+    if (!can_relax()) {
+      if (have_solution_) {  // data changed since the solve
+        sol_.status = Status::kUnsolved;
+        sol_.converged = 0;
+      }
+      return sol_;
     }
-    if (stale || !have_solution_) return sol_;
     if (!rw_.ready) relax_alloc();
     RelaxWork& w = rw_;
     const double kappa_s = c_ * kappa;
@@ -670,7 +681,7 @@ class Solver {
   }
 
   // Marks only rows that actually changed.
-  void set_rows(const MatrixXd& R, int offset) {
+  void set_rows(const MatRef& R, int offset) {
     for (int i = 0; i < R.rows(); ++i) {
       bool same = true;
       for (int k = 0; k < n_ && same; ++k) same = Ct_(k, offset + i) == R(i, k);