From 47c5a792769927d7cb86b3f4313880f67ac37d7f Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 8 Oct 2026 13:11:09 +0200 Subject: [PATCH 1/8] Add a pure C interface --- CMakeLists.txt | 36 ++ README.md | 30 +- cmake/elastiqpConfig.cmake.in | 10 + cmake/elastiqp_install.cmake | 9 +- include/elastiqp/das.hpp | 12 +- include/elastiqp/elastiqp_c.h | 238 ++++++++++++++ src/elastiqp_c.cc | 599 ++++++++++++++++++++++++++++++++++ tests/test_c_api.c | 417 +++++++++++++++++++++++ tests/test_c_parity.cc | 245 ++++++++++++++ 9 files changed, 1593 insertions(+), 3 deletions(-) create mode 100644 include/elastiqp/elastiqp_c.h create mode 100644 src/elastiqp_c.cc create mode 100644 tests/test_c_api.c create mode 100644 tests/test_c_parity.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index f0b2ec7..296f726 100644 --- a/CMakeLists.txt +++ b/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 + $ + $) + 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 diff --git a/README.md b/README.md index b017788..5527130 100644 --- a/README.md +++ b/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,34 @@ while (running) { } ``` +### C + +The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Exit flags: `> 0` success, `0` unsolved, `< 0` failure. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. + +```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}; +ElastiQPResult res = {.x = x, .z = z}; // caller-owned outputs; NULL ones are skipped +elastiqp_quadprog(&res, &qp, NULL); // NULL: default settings + +// Repeated solves +ElastiQPWorkspace* work; +elastiqp_setup(&work, &qp, NULL); +while (running) { + elastiqp_update_q(work, q_k); elastiqp_update_h(work, h_k); elastiqp_update_b(work, b_k); + elastiqp_update_G_rows(work, first, count, rows_k); // changed rows + if (elastiqp_solve(&res, work) < 0) { /* ... */ } +} +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 diff --git a/cmake/elastiqpConfig.cmake.in b/cmake/elastiqpConfig.cmake.in index 0050ccb..1351299 100644 --- a/cmake/elastiqpConfig.cmake.in +++ b/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) diff --git a/cmake/elastiqp_install.cmake b/cmake/elastiqp_install.cmake index 216bc59..6a85f7f 100644 --- a/cmake/elastiqp_install.cmake +++ b/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) diff --git a/include/elastiqp/das.hpp b/include/elastiqp/das.hpp index 3dba2f9..81159ae 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -190,6 +190,15 @@ class Solver { rhs_dirty_ = true; } void set_G(const MatrixXd& G) { set_rows(G, m_); } + // Replace rows [first, first + R.rows()); read R in place. + template + void set_A_rows(int first, const Eigen::MatrixBase& R) { + set_rows(R, first); + } + template + void set_G_rows(int first, const Eigen::MatrixBase& R) { + set_rows(R, m_ + first); + } void set_h(const VectorXd& h) { rhs_.tail(p_) = h; rhs_dirty_ = true; @@ -670,7 +679,8 @@ class Solver { } // Marks only rows that actually changed. - void set_rows(const MatrixXd& R, int offset) { + template + void set_rows(const Eigen::MatrixBase& 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); diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h new file mode 100644 index 0000000..c16b9d7 --- /dev/null +++ b/include/elastiqp/elastiqp_c.h @@ -0,0 +1,238 @@ +/* + * ElastiQP C interface (dual active-set backend, see das.hpp). + * + * Solves + * + * minimize 0.5 x'Q x + q'x + penalty't + * subject to A x = b + * G x - t <= h, t >= 0 + * + * Dense row-major matrices: M(i, j) = M[i * cols + j]. + * Q is symmetric; only its upper triangle is read. + * Two-sided constraints and variable bounds use two rows of G (g and -g). + * Penalties must be > 0; values >= ELASTIQP_INF make the row hard (t = 0). + * Exit flags: > 0 success, 0 unsolved, < 0 failure. + * + * Use elastiqp_quadprog() for one-shot solves, or setup/update/solve/free + * for repeated solves. Input data is copied; caller buffers can be reused. + * Allocation-free update/solve coverage: tests/test_das_alloc.cc. + * + * Windows DLL users outside CMake: define ELASTIQP_C_SHARED before including. + */ + +#ifndef ELASTIQP_C_H +#define ELASTIQP_C_H + +#if defined(_WIN32) && defined(ELASTIQP_C_SHARED) +#ifdef ELASTIQP_C_BUILDING +#define ELASTIQP_C_API __declspec(dllexport) +#else +#define ELASTIQP_C_API __declspec(dllimport) +#endif +#elif defined(__GNUC__) +#define ELASTIQP_C_API __attribute__((visibility("default"))) +#else +#define ELASTIQP_C_API +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +/* Hard-row penalty threshold. */ +#define ELASTIQP_INF 1e30 + +/* Exit flags. */ +#define ELASTIQP_SOLVED 1 +#define ELASTIQP_UNSOLVED 0 /* no solve yet */ +#define ELASTIQP_MAX_ITER (-1) /* iteration limit reached */ +#define ELASTIQP_NUMERICS (-2) /* numerical failure */ +#define ELASTIQP_INFEASIBLE (-3) /* inconsistent equalities or hard rows */ +#define ELASTIQP_ERR_INVALID_ARG (-10) /* NULL, bad dimension, value, name */ +#define ELASTIQP_ERR_NOT_SOLVED (-11) /* relax() without a current solve */ +#define ELASTIQP_ERR_ALLOC (-12) /* out of memory */ +#define ELASTIQP_ERR_INTERNAL (-13) /* unexpected C++ exception */ + +typedef struct { + int n; /* number of variables, >= 1 */ + int m; /* number of equality constraints (rows of A) */ + int p; /* number of elastic inequality constraints (rows of G) */ + + const double* Q; /* n x n, upper triangle read; NULL: zero (LP) */ + const double* q; /* n; NULL: zero */ + const double* A; /* m x n, row-major */ + const double* b; /* m */ + const double* G; /* p x n, row-major */ + const double* h; /* p */ + const double* penalty; /* p, > 0; >= ELASTIQP_INF: hard row */ +} ElastiQPProblem; + +/* + * Booleans use 0 / nonzero. Fields are also accessible as named options; + * elastiqp_option_info() gives their ranges. Out-of-range values are rejected. + * + * Tolerances and iteration limits apply on the next solve. Scaling and + * factorization options (ruiz*, zero_tol, eps_prox, reuse_factorization) + * apply on refactorization; updating Q forces it. Disabling ruiz retains + * existing scaling. Explicit warm starts apply even with warm_start = 0. + */ +typedef struct { + double eps_abs; + double eps_rel; + double sing_tol; + double zero_tol; + double eps_prox; + double eta_prox; + double prox_relaxation; + int prox_escalations; + int max_iter; + int max_outer; + int warm_start; + int reuse_factorization; + + int ruiz; + int ruiz_max_iter; + double ruiz_tol; + double ruiz_refresh_ratio; + + int check_eq_consistency; + + double progress_tol; + int cycle_tol; + double refactor_tol; + + double relax_reg; + int relax_factor_retries; +} ElastiQPSettings; + +/* + * Caller-owned output arrays; NULL arrays are skipped. + * KKT: Q x + q + A'y + G'z = 0, 0 <= z <= penalty, t = max(G x - h, 0). + * z_t = penalty - z after solve; relax() computes it independently. + * On API errors, only exitflag is valid; statistics are zeroed. + */ +typedef struct { + double* x; /* n */ + double* y; /* m, equality duals */ + double* z; /* p, inequality duals */ + double* t; /* p, elastic slacks */ + double* z_t; /* p, duals of t >= 0 (INFINITY on hard rows) */ + + int exitflag; /* ELASTIQP_SOLVED, ... */ + int iter; /* inner active-set (or relax Newton) iterations */ + int outer_iter; /* proximal rounds */ + int n_active; /* rows with t = 0 and z > 0 */ + int n_saturated; /* rows with t > 0 (z at its penalty) */ + double fval; /* primal objective, penalty term included */ + double primal_res; + double dual_res; + double duality_gap; + + double setup_time; /* seconds; elastiqp_quadprog only */ + double solve_time; /* seconds */ +} ElastiQPResult; + +/* Opaque solver workspace. */ +typedef struct ElastiQPWorkspace ElastiQPWorkspace; + +ELASTIQP_C_API void elastiqp_default_settings(ElastiQPSettings* settings); + +/* Name of an exit flag ("solved", "max_iter", ..., "invalid_arg"). */ +ELASTIQP_C_API const char* elastiqp_status_name(int exitflag); + +/* One-shot solve; settings may be NULL (defaults). Returns res->exitflag. */ +ELASTIQP_C_API int elastiqp_quadprog(ElastiQPResult* res, + const ElastiQPProblem* qp, + const ElastiQPSettings* settings); + +/* + * Creates a workspace for qp (settings may be NULL). Returns 0 and sets *work, + * or a negative error and sets *work to NULL. + */ +ELASTIQP_C_API int elastiqp_setup(ElastiQPWorkspace** work, + const ElastiQPProblem* qp, + const ElastiQPSettings* settings); + +/* Solves the current problem, warm-started from the previous solve when + * settings.warm_start is set. Returns res->exitflag. */ +ELASTIQP_C_API int elastiqp_solve(ElastiQPResult* res, ElastiQPWorkspace* work); + +ELASTIQP_C_API void elastiqp_free(ElastiQPWorkspace* work); + +/* Problem dimensions of a workspace; any output may be NULL. */ +ELASTIQP_C_API void elastiqp_get_dims(const ElastiQPWorkspace* work, int* n, + int* m, int* p); + +ELASTIQP_C_API int elastiqp_get_settings(const ElastiQPWorkspace* work, + ElastiQPSettings* settings); +ELASTIQP_C_API int elastiqp_set_settings(ElastiQPWorkspace* work, + const ElastiQPSettings* settings); + +/* + * Named workspace options, independent of the settings struct layout. + * Values use double: INT requires an integer, BOOL requires 0 or 1. + * All values must be in [lower, upper]. Returns 0 or ELASTIQP_ERR_INVALID_ARG; + * rejected options remain unchanged. + */ +#define ELASTIQP_OPTION_DOUBLE 0 +#define ELASTIQP_OPTION_INT 1 +#define ELASTIQP_OPTION_BOOL 2 + +ELASTIQP_C_API int elastiqp_set_option(ElastiQPWorkspace* work, + const char* name, double value); +ELASTIQP_C_API int elastiqp_get_option(const ElastiQPWorkspace* work, + const char* name, double* value); + +/* Option count and metadata for 0 <= index < count. Outputs may be NULL. + * option_info returns 0 or ELASTIQP_ERR_INVALID_ARG; strings are static. */ +ELASTIQP_C_API int elastiqp_num_options(void); +ELASTIQP_C_API int elastiqp_option_info(int index, const char** name, + int* type, double* default_value, + double* lower, double* upper, + const char** description); + +/* + * Replace data using setup dimensions and ElastiQPProblem conventions. + * NULL Q or q means zero. Only changed A/G rows are refreshed. + * Zero-sized updates do nothing. Returns 0 or a negative error; + * rejected updates leave the workspace unchanged. + */ +ELASTIQP_C_API int elastiqp_update_Q(ElastiQPWorkspace* work, const double* Q); +ELASTIQP_C_API int elastiqp_update_q(ElastiQPWorkspace* work, const double* q); +ELASTIQP_C_API int elastiqp_update_A(ElastiQPWorkspace* work, const double* A); +ELASTIQP_C_API int elastiqp_update_b(ElastiQPWorkspace* work, const double* b); +ELASTIQP_C_API int elastiqp_update_G(ElastiQPWorkspace* work, const double* G); +ELASTIQP_C_API int elastiqp_update_h(ElastiQPWorkspace* work, const double* h); +ELASTIQP_C_API int elastiqp_update_penalty(ElastiQPWorkspace* work, + const double* penalty); + +/* + * Replace rows [first, first + count) with a count x n row-major block. + * Cost: O(count * n). + */ +ELASTIQP_C_API int elastiqp_update_A_rows(ElastiQPWorkspace* work, int first, + int count, const double* rows); +ELASTIQP_C_API int elastiqp_update_G_rows(ElastiQPWorkspace* work, int first, + int count, const double* rows); + +/* + * Seed the next solve with (x, y, z). Requires x; NULL y or z means zero. + */ +ELASTIQP_C_API int elastiqp_set_warm_start(ElastiQPWorkspace* work, + const double* x, const double* y, + const double* z); + +/* + * Compute the kappa-relaxed central point for differentiation (kappa > 0). + * Requires a successful solve on current data, else ELASTIQP_ERR_NOT_SOLVED. + * Preserves the working set and warm start; p == 0 returns the solution. + * Suggested: tol = 1e-6, max_iter = 50. Returns res->exitflag. + */ +ELASTIQP_C_API int elastiqp_relax(ElastiQPResult* res, ElastiQPWorkspace* work, + double kappa, double tol, int max_iter); + +#ifdef __cplusplus +} +#endif + +#endif /* ELASTIQP_C_H */ diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc new file mode 100644 index 0000000..8029e0e --- /dev/null +++ b/src/elastiqp_c.cc @@ -0,0 +1,599 @@ +// C interface to the dual active-set backend; see elastiqp/elastiqp_c.h. + +#include "elastiqp/elastiqp_c.h" + +#include +#include +#include +#include +#include + +#include "elastiqp/das.hpp" + +using Eigen::MatrixXd; +using Eigen::VectorXd; +using RowMajor = Eigen::Matrix; +using RowMajorMap = Eigen::Map; +using VectorMap = Eigen::Map; +namespace das = elastiqp::das; + +// Preallocated staging buffers; A/G rows are read directly from caller data. +struct ElastiQPWorkspace { + das::Solver solver; + MatrixXd Q; + VectorXd q, b, h, penalty, x, y, z; + // Successful solve on current data; required by relax(). + bool solved = false; +}; + +namespace { + +using Clock = std::chrono::steady_clock; + +double Seconds(Clock::time_point since) { + return std::chrono::duration(Clock::now() - since).count(); +} + +int ExitFlag(elastiqp::Status s) { + switch (s) { + case elastiqp::Status::kSolved: + return ELASTIQP_SOLVED; + case elastiqp::Status::kUnsolved: + return ELASTIQP_UNSOLVED; + case elastiqp::Status::kMaxIter: + return ELASTIQP_MAX_ITER; + case elastiqp::Status::kNumerics: + return ELASTIQP_NUMERICS; + case elastiqp::Status::kInfeasible: + return ELASTIQP_INFEASIBLE; + } + return ELASTIQP_ERR_INTERNAL; +} + +// Shared metadata for settings conversion, validation, and named options. +// Add new fields here and in tests/test_c_parity.cc. +struct Option { + const char* name; + int type; + double lower, upper; + const char* description; + double das::Settings::*d_cpp; + double ElastiQPSettings::*d_c; + int das::Settings::*i_cpp; + bool das::Settings::*b_cpp; + int ElastiQPSettings::*i_c; // INT and BOOL options + + double Get(const das::Settings& s) const { + if (type == ELASTIQP_OPTION_DOUBLE) return s.*d_cpp; + if (type == ELASTIQP_OPTION_INT) return s.*i_cpp; + return s.*b_cpp ? 1.0 : 0.0; + } + // value must be Valid(). + void Set(das::Settings& s, double value) const { + if (type == ELASTIQP_OPTION_DOUBLE) { + s.*d_cpp = value; + } else if (type == ELASTIQP_OPTION_INT) { + s.*i_cpp = static_cast(value); + } else { + s.*b_cpp = value != 0.0; + } + } + // C booleans are 0 / nonzero; read them as 0 / 1. + double Get(const ElastiQPSettings& c) const { + if (type == ELASTIQP_OPTION_DOUBLE) return c.*d_c; + if (type == ELASTIQP_OPTION_INT) return c.*i_c; + return c.*i_c != 0 ? 1.0 : 0.0; + } + void Set(ElastiQPSettings& c, double value) const { + if (type == ELASTIQP_OPTION_DOUBLE) { + c.*d_c = value; + } else { + c.*i_c = static_cast(value); + } + } + bool Valid(double value) const { + if (!(value >= lower && value <= upper)) return false; // NaN fails + return type == ELASTIQP_OPTION_DOUBLE || value == std::floor(value); + } +}; + +constexpr double kInf = std::numeric_limits::infinity(); +constexpr double kIntMax = std::numeric_limits::max(); +constexpr double kPositive = std::numeric_limits::min(); + +constexpr Option D(const char* name, double das::Settings::*cpp, + double ElastiQPSettings::*c, double lower, + const char* description) { + return {name, ELASTIQP_OPTION_DOUBLE, lower, kInf, description, + cpp, c, nullptr, nullptr, nullptr}; +} +constexpr Option I(const char* name, int das::Settings::*cpp, + int ElastiQPSettings::*c, const char* description) { + return {name, ELASTIQP_OPTION_INT, 0.0, kIntMax, description, + nullptr, nullptr, cpp, nullptr, c}; +} +constexpr Option B(const char* name, bool das::Settings::*cpp, + int ElastiQPSettings::*c, const char* description) { + return {name, ELASTIQP_OPTION_BOOL, 0.0, 1.0, description, + nullptr, nullptr, nullptr, cpp, c}; +} + +#define ELASTIQP_D(f, lower, desc) \ + D(#f, &das::Settings::f, &ElastiQPSettings::f, lower, desc) +#define ELASTIQP_I(f, desc) \ + I(#f, &das::Settings::f, &ElastiQPSettings::f, desc) +#define ELASTIQP_B(f, desc) \ + B(#f, &das::Settings::f, &ElastiQPSettings::f, desc) + +constexpr Option kOptions[] = { + ELASTIQP_D(eps_abs, kPositive, "Absolute tolerance."), + ELASTIQP_D(eps_rel, 0.0, "Relative tolerance."), + ELASTIQP_D(sing_tol, 0.0, + "Working-set LDL' pivot below this = dependent row."), + ELASTIQP_D(zero_tol, 0.0, + "Min/max Cholesky pivot ratio of Q_s accepted before adding " + "a prox shift."), + ELASTIQP_D(eps_prox, 0.0, + "Initial prox shift (x max|diag Q_s|) when Q_s is not PD; 0 " + "disables."), + ELASTIQP_D(eta_prox, 0.0, + "Outer-loop tolerance on eps*|x - xc|_inf; 0: eps_abs."), + ELASTIQP_D(prox_relaxation, 0.0, + "Over-relaxation of the prox center when the working set " + "stopped changing; <= 1 disables."), + ELASTIQP_I(prox_escalations, + "Max 100x prox-shift increases after an inner numerical " + "failure."), + ELASTIQP_I(max_iter, "Inner active-set iterations per outer iteration."), + ELASTIQP_I(max_outer, "Outer (proximal) iterations."), + ELASTIQP_B(warm_start, "Warm-start from the previous solve."), + ELASTIQP_B(reuse_factorization, + "Keep the Cholesky of Q_s across solves when only rows or " + "vectors changed."), + ELASTIQP_B(ruiz, "Ruiz equilibration."), + ELASTIQP_I(ruiz_max_iter, "Ruiz sweeps."), + ELASTIQP_D(ruiz_tol, 0.0, "Ruiz convergence tolerance."), + ELASTIQP_D(ruiz_refresh_ratio, 0.0, + "Re-equilibrate when scaling drift exceeds this; 0: never."), + ELASTIQP_B(check_eq_consistency, + "Return infeasible early on inconsistent equalities."), + ELASTIQP_D(progress_tol, 0.0, + "Absolute dual-objective increase counted as progress by the " + "cycle guard."), + ELASTIQP_I(cycle_tol, + "Stalled iterations after a removal before repair / " + "numerics."), + ELASTIQP_D(refactor_tol, 0.0, + "At optimality, refactor once if the smallest pivot is below " + "this."), + ELASTIQP_D(relax_reg, 0.0, "relax(): Newton regularization."), + ELASTIQP_I(relax_factor_retries, + "relax(): x100 regularization retries on a failed " + "factorization."), +}; + +#undef ELASTIQP_D +#undef ELASTIQP_I +#undef ELASTIQP_B + +constexpr int kNumOptions = sizeof(kOptions) / sizeof(kOptions[0]); + +const Option* FindOption(const char* name) { + if (name == nullptr) return nullptr; + for (const Option& o : kOptions) { + if (std::strcmp(o.name, name) == 0) return &o; + } + return nullptr; +} + +bool ValidSettings(const ElastiQPSettings& c) { + for (const Option& o : kOptions) { + if (!o.Valid(o.Get(c))) return false; + } + return true; +} + +das::Settings ToCpp(const ElastiQPSettings& c) { + das::Settings s; + for (const Option& o : kOptions) o.Set(s, o.Get(c)); + return s; +} + +ElastiQPSettings ToC(const das::Settings& s) { + ElastiQPSettings c{}; + for (const Option& o : kOptions) o.Set(c, o.Get(s)); + return c; +} + +// Positive penalties only; rejects NaN. +bool ValidPenalty(const double* w, int p) { + for (int i = 0; i < p; ++i) { + if (!(w[i] > 0.0)) return false; + } + return true; +} + +bool Valid(const ElastiQPProblem* qp) { + if (qp == nullptr || qp->n < 1 || qp->m < 0 || qp->p < 0) return false; + if (qp->m > 0 && (qp->A == nullptr || qp->b == nullptr)) return false; + if (qp->p > 0 && + (qp->G == nullptr || qp->h == nullptr || qp->penalty == nullptr || + !ValidPenalty(qp->penalty, qp->p))) { + return false; + } + return true; +} + +// Upper triangle of the row-major Q, mirrored; zero for NULL. +void StageQ(MatrixXd& dst, const double* Q, int n) { + if (Q) { + dst = RowMajorMap(Q, n, n).selfadjointView(); + } else { + dst.setZero(); + } +} + +// Penalties at or above ELASTIQP_INF become infinite (hard rows). +void StagePenalty(VectorXd& dst, const double* w) { + for (Eigen::Index i = 0; i < dst.size(); ++i) { + dst[i] = w[i] >= ELASTIQP_INF ? std::numeric_limits::infinity() + : w[i]; + } +} + +// Clear statistics on API errors. +int Fail(ElastiQPResult* res, int flag) { + res->exitflag = flag; + res->iter = res->outer_iter = res->n_active = res->n_saturated = 0; + res->fval = res->primal_res = res->dual_res = res->duality_gap = 0.0; + res->setup_time = res->solve_time = 0.0; + return flag; +} + +// Fills the caller's buffers from sol; returns the exit flag. +int Extract(ElastiQPResult* res, const elastiqp::Solution& sol) { + using Eigen::Map; + if (res->x) Map(res->x, sol.x.size()) = sol.x; + if (res->y) Map(res->y, sol.y.size()) = sol.y; + if (res->z) Map(res->z, sol.z.size()) = sol.z; + if (res->t) Map(res->t, sol.t.size()) = sol.t; + if (res->z_t) Map(res->z_t, sol.z_t.size()) = sol.z_t; + res->exitflag = ExitFlag(sol.status); + res->iter = sol.iters; + res->outer_iter = sol.outer_iters; + res->n_active = sol.n_active; + res->n_saturated = sol.n_saturated; + res->fval = sol.primal_obj; + res->primal_res = sol.primal_res; + res->dual_res = sol.dual_res; + res->duality_gap = sol.duality_gap; + return res->exitflag; +} + +// Runs f, mapping C++ exceptions to error codes. +template +int Guard(F&& f) { + try { + return f(); + } catch (const std::bad_alloc&) { + return ELASTIQP_ERR_ALLOC; + } catch (...) { + return ELASTIQP_ERR_INTERNAL; + } +} + +int SetupImpl(ElastiQPWorkspace& w, const ElastiQPProblem& qp, + const ElastiQPSettings* settings) { + const int n = qp.n, m = qp.m, p = qp.p; + w.Q.resize(n, n); + StageQ(w.Q, qp.Q, n); + if (qp.q) { + w.q = VectorMap(qp.q, n); + } else { + w.q.setZero(n); + } + MatrixXd A(m, n), G(p, n); + w.b.resize(m); + if (m > 0) { + A = RowMajorMap(qp.A, m, n); + w.b = VectorMap(qp.b, m); + } + w.h.resize(p); + w.penalty.resize(p); + if (p > 0) { + G = RowMajorMap(qp.G, p, n); + w.h = VectorMap(qp.h, p); + StagePenalty(w.penalty, qp.penalty); + } + w.x.resize(n); + w.y.resize(m); + w.z.resize(p); + if (settings) w.solver.settings = ToCpp(*settings); + w.solver.setup(w.Q, w.q, A, w.b, G, w.h, w.penalty); + return 0; +} + +// Validates a [first, first + count) row block of a matrix with `rows` rows. +bool ValidRows(int first, int count, int rows, const double* data) { + return first >= 0 && count >= 0 && first <= rows - count && + (count == 0 || data != nullptr); +} + +} // namespace + +extern "C" { + +void elastiqp_default_settings(ElastiQPSettings* settings) { + if (settings) *settings = ToC(das::Settings{}); +} + +const char* elastiqp_status_name(int exitflag) { + switch (exitflag) { + case ELASTIQP_SOLVED: + return "solved"; + case ELASTIQP_UNSOLVED: + return "unsolved"; + case ELASTIQP_MAX_ITER: + return "max_iter"; + case ELASTIQP_NUMERICS: + return "numerics"; + case ELASTIQP_INFEASIBLE: + return "infeasible"; + case ELASTIQP_ERR_INVALID_ARG: + return "invalid_arg"; + case ELASTIQP_ERR_NOT_SOLVED: + return "not_solved"; + case ELASTIQP_ERR_ALLOC: + return "alloc"; + case ELASTIQP_ERR_INTERNAL: + return "internal"; + } + return "?"; +} + +int elastiqp_quadprog(ElastiQPResult* res, const ElastiQPProblem* qp, + const ElastiQPSettings* settings) { + if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; + if (!Valid(qp) || (settings && !ValidSettings(*settings))) { + return Fail(res, ELASTIQP_ERR_INVALID_ARG); + } + double setup_time = 0.0, solve_time = 0.0; + const int flag = Guard([&] { + ElastiQPWorkspace w; + const Clock::time_point t0 = Clock::now(); + SetupImpl(w, *qp, settings); + setup_time = Seconds(t0); + const Clock::time_point t1 = Clock::now(); + const elastiqp::Solution& sol = w.solver.solve(); + solve_time = Seconds(t1); + return Extract(res, sol); + }); + if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + res->setup_time = setup_time; + res->solve_time = solve_time; + return flag; +} + +int elastiqp_setup(ElastiQPWorkspace** work, const ElastiQPProblem* qp, + const ElastiQPSettings* settings) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + *work = nullptr; + if (!Valid(qp) || (settings && !ValidSettings(*settings))) { + return ELASTIQP_ERR_INVALID_ARG; + } + ElastiQPWorkspace* w = new (std::nothrow) ElastiQPWorkspace; + if (w == nullptr) return ELASTIQP_ERR_ALLOC; + const int ret = Guard([&] { return SetupImpl(*w, *qp, settings); }); + if (ret != 0) { + delete w; + return ret; + } + *work = w; + return 0; +} + +int elastiqp_solve(ElastiQPResult* res, ElastiQPWorkspace* work) { + if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; + if (work == nullptr) return Fail(res, ELASTIQP_ERR_INVALID_ARG); + double solve_time = 0.0; + const int flag = Guard([&] { + const Clock::time_point t0 = Clock::now(); + const elastiqp::Solution& sol = work->solver.solve(); + solve_time = Seconds(t0); + return Extract(res, sol); + }); + work->solved = flag == ELASTIQP_SOLVED; + if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + res->setup_time = 0.0; + res->solve_time = solve_time; + return flag; +} + +void elastiqp_free(ElastiQPWorkspace* work) { delete work; } + +void elastiqp_get_dims(const ElastiQPWorkspace* work, int* n, int* m, int* p) { + if (n) *n = work ? work->solver.n() : 0; + if (m) *m = work ? work->solver.m() : 0; + if (p) *p = work ? work->solver.p() : 0; +} + +int elastiqp_get_settings(const ElastiQPWorkspace* work, + ElastiQPSettings* settings) { + if (work == nullptr || settings == nullptr) return ELASTIQP_ERR_INVALID_ARG; + *settings = ToC(work->solver.settings); + return 0; +} + +int elastiqp_set_settings(ElastiQPWorkspace* work, + const ElastiQPSettings* settings) { + if (work == nullptr || settings == nullptr || !ValidSettings(*settings)) { + return ELASTIQP_ERR_INVALID_ARG; + } + work->solver.settings = ToCpp(*settings); + return 0; +} + +int elastiqp_set_option(ElastiQPWorkspace* work, const char* name, + double value) { + const Option* o = FindOption(name); + if (work == nullptr || o == nullptr || !o->Valid(value)) { + return ELASTIQP_ERR_INVALID_ARG; + } + o->Set(work->solver.settings, value); + return 0; +} + +int elastiqp_get_option(const ElastiQPWorkspace* work, const char* name, + double* value) { + const Option* o = FindOption(name); + if (work == nullptr || o == nullptr || value == nullptr) { + return ELASTIQP_ERR_INVALID_ARG; + } + *value = o->Get(work->solver.settings); + return 0; +} + +int elastiqp_num_options(void) { return kNumOptions; } + +int elastiqp_option_info(int index, const char** name, int* type, + double* default_value, double* lower, double* upper, + const char** description) { + if (index < 0 || index >= kNumOptions) return ELASTIQP_ERR_INVALID_ARG; + const Option& o = kOptions[index]; + if (name) *name = o.name; + if (type) *type = o.type; + if (default_value) *default_value = o.Get(das::Settings{}); + if (lower) *lower = o.lower; + if (upper) *upper = o.upper; + if (description) *description = o.description; + return 0; +} + +int elastiqp_update_Q(ElastiQPWorkspace* work, const double* Q) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + StageQ(work->Q, Q, work->solver.n()); + work->solver.set_Q(work->Q); + work->solved = false; + return 0; +} + +int elastiqp_update_q(ElastiQPWorkspace* work, const double* q) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + if (q) { + work->q = VectorMap(q, work->solver.n()); + } else { + work->q.setZero(); + } + work->solver.set_q(work->q); + work->solved = false; + return 0; +} + +int elastiqp_update_A(ElastiQPWorkspace* work, const double* A) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + return elastiqp_update_A_rows(work, 0, work->solver.m(), A); +} + +int elastiqp_update_b(ElastiQPWorkspace* work, const double* b) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + const int m = work->solver.m(); + if (m == 0) return 0; + if (b == nullptr) return ELASTIQP_ERR_INVALID_ARG; + work->b = VectorMap(b, m); + work->solver.set_b(work->b); + work->solved = false; + return 0; +} + +int elastiqp_update_G(ElastiQPWorkspace* work, const double* G) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + return elastiqp_update_G_rows(work, 0, work->solver.p(), G); +} + +int elastiqp_update_h(ElastiQPWorkspace* work, const double* h) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + const int p = work->solver.p(); + if (p == 0) return 0; + if (h == nullptr) return ELASTIQP_ERR_INVALID_ARG; + work->h = VectorMap(h, p); + work->solver.set_h(work->h); + work->solved = false; + return 0; +} + +int elastiqp_update_penalty(ElastiQPWorkspace* work, const double* penalty) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + const int p = work->solver.p(); + if (p == 0) return 0; + if (penalty == nullptr || !ValidPenalty(penalty, p)) { + return ELASTIQP_ERR_INVALID_ARG; + } + StagePenalty(work->penalty, penalty); + work->solver.set_penalty(work->penalty); + work->solved = false; + return 0; +} + +int elastiqp_update_A_rows(ElastiQPWorkspace* work, int first, int count, + const double* rows) { + if (work == nullptr || !ValidRows(first, count, work->solver.m(), rows)) { + return ELASTIQP_ERR_INVALID_ARG; + } + if (count == 0) return 0; + work->solver.set_A_rows(first, RowMajorMap(rows, count, work->solver.n())); + work->solved = false; + return 0; +} + +int elastiqp_update_G_rows(ElastiQPWorkspace* work, int first, int count, + const double* rows) { + if (work == nullptr || !ValidRows(first, count, work->solver.p(), rows)) { + return ELASTIQP_ERR_INVALID_ARG; + } + if (count == 0) return 0; + work->solver.set_G_rows(first, RowMajorMap(rows, count, work->solver.n())); + work->solved = false; + return 0; +} + +int elastiqp_set_warm_start(ElastiQPWorkspace* work, const double* x, + const double* y, const double* z) { + if (work == nullptr || x == nullptr) return ELASTIQP_ERR_INVALID_ARG; + const int n = work->solver.n(), m = work->solver.m(), p = work->solver.p(); + work->x = VectorMap(x, n); + if (y) { + work->y = VectorMap(y, m); + } else { + work->y.setZero(); + } + if (z) { + work->z = VectorMap(z, p); + } else { + work->z.setZero(); + } + work->solver.set_warm_start(work->x, work->y, work->z); + return 0; +} + +int elastiqp_relax(ElastiQPResult* res, ElastiQPWorkspace* work, double kappa, + double tol, int max_iter) { + if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; + if (work == nullptr || !(kappa > 0.0) || !(tol > 0.0) || max_iter < 0) { + return Fail(res, ELASTIQP_ERR_INVALID_ARG); + } + if (!work->solved) return Fail(res, ELASTIQP_ERR_NOT_SOLVED); + double solve_time = 0.0; + const int flag = Guard([&] { + const Clock::time_point t0 = Clock::now(); + const elastiqp::Solution& sol = work->solver.relax(kappa, tol, max_iter); + solve_time = Seconds(t0); + return Extract(res, sol); + }); + if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + res->setup_time = 0.0; + res->solve_time = solve_time; + return flag; +} + +} // extern "C" diff --git a/tests/test_c_api.c b/tests/test_c_api.c new file mode 100644 index 0000000..97c84b3 --- /dev/null +++ b/tests/test_c_api.c @@ -0,0 +1,417 @@ +/* C API: analytic solutions, updates, warm starts, and validation. */ + +#include +#include +#include + +#include "elastiqp/elastiqp_c.h" + +static int g_all_ok = 1; + +static void check(const char* name, int ok) { + printf(" %-52s %s\n", name, ok ? "OK" : "FAIL"); + g_all_ok &= ok; +} + +static int near(const double* a, const double* b, int k, double tol) { + int i; + for (i = 0; i < k; ++i) { + if (!(fabs(a[i] - b[i]) <= tol)) return 0; + } + return 1; +} + +static ElastiQPSettings tight(void) { + ElastiQPSettings s; + elastiqp_default_settings(&s); + s.eps_abs = 1e-10; + return s; +} + +/* min 0.5|x|^2 - x1 - x2 s.t. x1 + x2 <= 1: x = (0.5, 0.5), z = 0.5. */ +static void one_shot(void) { + const double Q[] = {1, 0, 0, 1}, q[] = {-1, -1}; + const double G[] = {1, 1}, h[] = {1}, w[] = {10}; + const double x_ref[] = {0.5, 0.5}, z_ref[] = {0.5}, t_ref[] = {0}; + double x[2], z[1], t[1]; + ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; + ElastiQPResult res = {x, NULL, z, t}; + ElastiQPSettings s = tight(); + int flag = elastiqp_quadprog(&res, &qp, &s); + check("one-shot: solved", flag == ELASTIQP_SOLVED && res.exitflag == flag); + check("one-shot: x, z, t", near(x, x_ref, 2, 1e-8) && + near(z, z_ref, 1, 1e-8) && + near(t, t_ref, 1, 1e-12)); + check("one-shot: active row counted", res.n_active == 1 && + res.n_saturated == 0); + check("one-shot: objective", fabs(res.fval + 0.75) < 1e-8); + check("one-shot: default settings (NULL)", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_SOLVED && + near(x, x_ref, 2, 1e-5)); +} + +/* + * x <= -1, x >= 1, penalty = 0.5, Q = 1: x = 0, both rows saturated. + * Making x >= 1 hard gives x = 1. + */ +static void elastic_and_hard(void) { + const double Q[] = {1}, G[] = {1, -1}, h[] = {-1, -1}; + double w[] = {0.5, 0.5}; + double x[1], z[2], t[2]; + const double t_ref[] = {1, 1}, z_ref[] = {0.5, 0.5}; + const double t_hard[] = {2, 0}; + ElastiQPProblem qp = {1, 0, 2, Q, NULL, NULL, NULL, G, h, w}; + ElastiQPResult res = {x, NULL, z, t}; + ElastiQPSettings s = tight(); + int flag = elastiqp_quadprog(&res, &qp, &s); + check("conflict: solved", flag == ELASTIQP_SOLVED); + check("conflict: x = 0, t = 1, z at penalty", + fabs(x[0]) < 1e-8 && near(t, t_ref, 2, 1e-8) && + near(z, z_ref, 2, 1e-8)); + check("conflict: both rows saturated", res.n_saturated == 2); + w[1] = INFINITY; + flag = elastiqp_quadprog(&res, &qp, &s); + check("hard row: solved", flag == ELASTIQP_SOLVED); + check("hard row: x = 1, hard slack 0", + fabs(x[0] - 1) < 1e-8 && near(t, t_hard, 2, 1e-8)); +} + +/* x1 = 0.2, x2 <= h: solution x2 = min(h, 1). */ +static void workspace(void) { + const double Q[] = {1, 0, 0, 1}, q[] = {-1, -1}; + const double A[] = {1, 0}, b[] = {0.2}; + const double G[] = {0, 1}, w[] = {10}; + double h[] = {0.8}; + double x[2], y[1], z[1], t[1]; + ElastiQPProblem qp = {2, 1, 1, Q, q, A, b, G, h, w}; + ElastiQPResult res = {x, y, z, t}; + ElastiQPSettings s = tight(); + ElastiQPWorkspace* work = NULL; + int n = 0, m = 0, p = 0, flag; + const double x1[] = {0.2, 0.8}, y1[] = {0.8}, z1[] = {0.2}; + const double x2[] = {0.2, 1.0}, z2[] = {0.0}; + const double x3[] = {-0.5, 1.0}; + + check("setup: returns 0", elastiqp_setup(&work, &qp, &s) == 0 && work); + elastiqp_get_dims(work, &n, &m, &p); + check("setup: dims", n == 2 && m == 1 && p == 1); + flag = elastiqp_solve(&res, work); + check("solve: solved", flag == ELASTIQP_SOLVED); + check("solve: x, y, z", near(x, x1, 2, 1e-8) && near(y, y1, 1, 1e-8) && + near(z, z1, 1, 1e-8)); + + h[0] = 2.0; /* workspace retains its copy */ + elastiqp_solve(&res, work); + check("data copied at setup", near(x, x1, 2, 1e-8)); + check("update_h", elastiqp_update_h(work, h) == 0); + elastiqp_solve(&res, work); + check("update_h: constraint released", + near(x, x2, 2, 1e-8) && near(z, z2, 1, 1e-8)); + + { + const double b2[] = {-0.5}; + elastiqp_update_b(work, b2); + elastiqp_solve(&res, work); + check("update_b", res.exitflag == ELASTIQP_SOLVED && near(x, x3, 2, 1e-8)); + } + { + const double q2[] = {-1, -3}, x4[] = {-0.5, 2.0}; + elastiqp_update_q(work, q2); + elastiqp_solve(&res, work); + check("update_q (row active at h = 2)", near(x, x4, 2, 1e-8)); + } + { + /* Penalty 0.5 < marginal cost 1 at x2 = 2: the row gives way. */ + const double w2[] = {0.5}, x5[] = {-0.5, 2.5}, t5[] = {0.5}; + elastiqp_update_penalty(work, w2); + elastiqp_solve(&res, work); + check("update_penalty: row saturates", + near(x, x5, 2, 1e-8) && near(t, t5, 1, 1e-8) && + res.n_saturated == 1); + } + { + const double Q2[] = {2, 0, 0, 2}, A2[] = {0, 1}, G2[] = {1, 0}; + const double x6[] = {0.5, -0.5}; + /* min |x|^2 - x1 - 3 x2 s.t. x2 = -0.5, x1 <= 2 (inactive). */ + elastiqp_update_Q(work, Q2); + elastiqp_update_A(work, A2); + elastiqp_update_G(work, G2); + elastiqp_solve(&res, work); + check("update_Q/A/G", near(x, x6, 2, 1e-8)); + } + { + ElastiQPSettings got; + const double xw[] = {0.5, -0.5}; + elastiqp_get_settings(work, &got); + check("get_settings", got.eps_abs == 1e-10 && got.warm_start); + got.max_iter = 0; + elastiqp_set_settings(work, &got); + elastiqp_set_warm_start(work, xw, NULL, NULL); + elastiqp_solve(&res, work); + check("set_settings: max_iter = 0 fails", + res.exitflag == ELASTIQP_MAX_ITER && res.exitflag < 0); + elastiqp_set_settings(work, &s); + elastiqp_solve(&res, work); + check("set_settings: restored", res.exitflag == ELASTIQP_SOLVED); + } + elastiqp_free(work); +} + +/* + * Row-major layout: G is 2 x 3. min 0.5|x|^2 + q'x with + * G = [[1, 2, 3], [0, 0, 1]], h = [0, 10] (second row inactive). The + * minimizer of 0.5|x - c|^2 over g'x <= 0 is c - (g'c / g'g) g. + */ +static void row_major(void) { + const double Q[] = {1, 0, 0, 0, 1, 0, 0, 0, 1}, q[] = {-1, -1, -1}; + const double G[] = {1, 2, 3, 0, 0, 1}, h[] = {0, 10}, w[] = {100, 100}; + /* c = (1,1,1), g'c = 6, g'g = 14 */ + const double x_ref[] = {1 - 6.0 / 14, 1 - 12.0 / 14, 1 - 18.0 / 14}; + double x[3]; + ElastiQPProblem qp = {3, 0, 2, Q, q, NULL, NULL, G, h, w}; + ElastiQPResult res = {x, NULL, NULL, NULL}; + ElastiQPSettings s = tight(); + elastiqp_quadprog(&res, &qp, &s); + check("row-major G", res.exitflag == ELASTIQP_SOLVED && + near(x, x_ref, 3, 1e-8)); +} + +/* Only the upper triangle of Q is read: full and upper-only Q agree. */ +static void q_upper_triangle(void) { + const double Qfull[] = {2, 1, 1, 2}, Qupper[] = {2, 1, -99, 2}; + const double q[] = {-1, 0}; + double x1[2], x2[2]; + ElastiQPProblem qp = {2, 0, 0, Qfull, q, NULL, NULL, NULL, NULL, NULL}; + ElastiQPResult r1 = {x1}, r2 = {x2}; + /* 2 x1 + x2 = 1, x1 + 2 x2 = 0 */ + const double x_ref[] = {2.0 / 3, -1.0 / 3}; + elastiqp_quadprog(&r1, &qp, NULL); + qp.Q = Qupper; + elastiqp_quadprog(&r2, &qp, NULL); + check("Q: full matrix", near(x1, x_ref, 2, 1e-10)); + check("Q: upper triangle only, lower ignored", near(x2, x_ref, 2, 1e-10)); +} + +/* ELASTIQP_INF enforces a hard row. */ +static void inf_convention(void) { + const double Q[] = {1}, G[] = {1, -1}, h[] = {-1, -1}; + const double w[] = {0.5, ELASTIQP_INF}; + double x[1], t[2], z_t[2]; + ElastiQPProblem qp = {1, 0, 2, Q, NULL, NULL, NULL, G, h, w}; + ElastiQPResult res = {x, NULL, NULL, t, z_t}; + elastiqp_quadprog(&res, &qp, NULL); + check("ELASTIQP_INF penalty is a hard row", + res.exitflag == ELASTIQP_SOLVED && fabs(x[0] - 1) < 1e-8 && + fabs(t[1]) < 1e-12 && isinf(z_t[1])); +} + +/* Rows of a 3 x 2 G replaced in blocks, against fresh one-shot solves. */ +static void row_updates(void) { + const double Q[] = {1, 0, 0, 1}, q[] = {-2, -2}, h[] = {1, 1, 1}; + const double w[] = {10, 10, 10}; + double G[] = {1, 0, 0, 1, 1, 1}; + const double rows[] = {0.5, 1, 2, 0.5}; + double x[2], x_ref[2]; + ElastiQPProblem qp = {2, 0, 3, Q, q, NULL, NULL, G, h, w}; + ElastiQPResult res = {x}, ref = {x_ref}; + ElastiQPWorkspace* work = NULL; + elastiqp_setup(&work, &qp, NULL); + elastiqp_solve(&res, work); + check("update_G_rows: rows 1..2", + elastiqp_update_G_rows(work, 1, 2, rows) == 0); + G[2] = 0.5, G[3] = 1, G[4] = 2, G[5] = 0.5; + elastiqp_solve(&res, work); + elastiqp_quadprog(&ref, &qp, NULL); + check("update_G_rows: matches a fresh solve", + res.exitflag == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-8)); + check("update_G_rows: out of range rejected", + elastiqp_update_G_rows(work, 2, 2, rows) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_update_G_rows(work, -1, 1, rows) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_update_A_rows(work, 0, 1, rows) == + ELASTIQP_ERR_INVALID_ARG); + check("update_G_rows: count 0 is a no-op", + elastiqp_update_G_rows(work, 3, 0, NULL) == 0); + elastiqp_free(work); +} + +static void relax(void) { + const double Q[] = {1, 0, 0, 1}, q[] = {-1, -1}, q2[] = {-1, -0.5}; + const double G[] = {1, 1}, h[] = {1}, w[] = {10}; + double x[2], z[1], t[1], z_t[1]; + ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; + ElastiQPResult res = {x, NULL, z, t, z_t}; + ElastiQPWorkspace* work = NULL; + elastiqp_setup(&work, &qp, NULL); + check("relax before solve: not_solved", + elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_ERR_NOT_SOLVED); + elastiqp_solve(&res, work); + check("solve: z_t = penalty - z", fabs(z_t[0] - (10 - z[0])) < 1e-12); + check("relax: kappa <= 0 rejected", + elastiqp_relax(&res, work, 0.0, 1e-8, 50) == + ELASTIQP_ERR_INVALID_ARG && + res.iter == 0); + elastiqp_relax(&res, work, 1e-3, 1e-8, 50); + /* At the relaxed point the inequality slack s = h - Gx + t is O(kappa). */ + check("relax: solved, interior point", + res.exitflag == ELASTIQP_SOLVED && 1 - x[0] - x[1] + t[0] > 0 && + z[0] > 0 && z[0] < 10 && z_t[0] > 0); + elastiqp_update_q(work, q2); + check("relax after an update: not_solved", + elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_ERR_NOT_SOLVED); + elastiqp_solve(&res, work); + check("relax after re-solve", + elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_SOLVED); + elastiqp_free(work); +} + +/* Options by name: listing, defaults, validation, effect on a solve. */ +static void options(void) { + const double Q[] = {1, 0, 0, 1}, q[] = {-1, -1}; + const double G[] = {1, 1}, h[] = {1}, w[] = {10}; + double x[2], v = 0, def = 0, lo = 0, hi = 0; + ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; + ElastiQPResult res = {x}; + ElastiQPSettings s; + ElastiQPWorkspace* work = NULL; + const char *name = NULL, *name2 = NULL, *desc = NULL; + int i, j, type = -1, ok = 1, unique = 1, n_opt = elastiqp_num_options(); + + elastiqp_setup(&work, &qp, NULL); + check("options: 22 listed", n_opt == 22); + for (i = 0; i < n_opt; ++i) { + ok &= elastiqp_option_info(i, &name, &type, &def, &lo, &hi, &desc) == 0; + ok &= name != NULL && desc != NULL && desc[0] != '\0'; + ok &= type >= ELASTIQP_OPTION_DOUBLE && type <= ELASTIQP_OPTION_BOOL; + ok &= lo <= def && def <= hi; + ok &= elastiqp_get_option(work, name, &v) == 0 && v == def; + for (j = 0; j < i; ++j) { + elastiqp_option_info(j, &name2, NULL, NULL, NULL, NULL, NULL); + unique &= strcmp(name, name2) != 0; + } + } + check("options: info consistent with workspace defaults", ok); + check("options: names unique", unique); + check("options: index out of range rejected", + elastiqp_option_info(n_opt, &name, NULL, NULL, NULL, NULL, NULL) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_option_info(-1, NULL, NULL, NULL, NULL, NULL, NULL) == + ELASTIQP_ERR_INVALID_ARG); + + check("set_option: unknown name rejected", + elastiqp_set_option(work, "eps_absolute", 1e-8) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_option(work, NULL, 1e-8) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_get_option(work, "nope", &v) == + ELASTIQP_ERR_INVALID_ARG); + check("set_option: non-integral int rejected", + elastiqp_set_option(work, "max_iter", 2.5) == + ELASTIQP_ERR_INVALID_ARG); + check("set_option: bool must be 0 or 1", + elastiqp_set_option(work, "ruiz", 2) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_option(work, "ruiz", 0) == 0 && + elastiqp_get_option(work, "ruiz", &v) == 0 && v == 0); + check("set_option: out of range / NaN rejected", + elastiqp_set_option(work, "eps_abs", 0) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_option(work, "eps_rel", -1) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_option(work, "max_iter", -1) == + ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_option(work, "eps_abs", NAN) == + ELASTIQP_ERR_INVALID_ARG); + check("set_option: rejected value leaves the option", + elastiqp_get_option(work, "eps_abs", &v) == 0 && v == 1e-6); + + elastiqp_set_option(work, "max_iter", 0); + check("set_option: max_iter = 0 takes effect", + elastiqp_solve(&res, work) == ELASTIQP_MAX_ITER); + elastiqp_set_option(work, "max_iter", 10000); + elastiqp_set_option(work, "eps_abs", 1e-10); + elastiqp_get_settings(work, &s); + check("set_option: visible through get_settings", + s.eps_abs == 1e-10 && s.ruiz == 0 && s.max_iter == 10000); + check("set_option: solves", elastiqp_solve(&res, work) == ELASTIQP_SOLVED && + fabs(x[0] - 0.5) < 1e-8); + + s.eps_abs = -1; + check("set_settings: out-of-range struct rejected", + elastiqp_set_settings(work, &s) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_get_option(work, "eps_abs", &v) == 0 && v == 1e-10); + elastiqp_free(work); + check("setup: out-of-range settings rejected", + elastiqp_setup(&work, &qp, &s) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_quadprog(&res, &qp, &s) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_default_settings(&s); + s.warm_start = 5; /* struct booleans are 0 / nonzero */ + elastiqp_setup(&work, &qp, &s); + check("struct bool nonzero accepted, reads as 1", + work != NULL && elastiqp_get_option(work, "warm_start", &v) == 0 && + v == 1); + elastiqp_free(work); +} + +static void invalid_args(void) { + const double Q[] = {1}, G[] = {1}, h[] = {0}; + double x[1]; + double bad[] = {0}; + ElastiQPProblem qp = {1, 0, 1, Q, NULL, NULL, NULL, G, h, NULL}; + ElastiQPResult res = {x}; + ElastiQPWorkspace* work = (ElastiQPWorkspace*)&qp; + res.iter = 7; + check("missing penalty rejected, stats cleared", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG && + res.exitflag == ELASTIQP_ERR_INVALID_ARG && res.iter == 0); + check("setup: rejects and clears *work", + elastiqp_setup(&work, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG && + work == NULL); + qp.penalty = bad; + check("penalty 0 rejected", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + bad[0] = -1; + check("negative penalty rejected", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + bad[0] = NAN; + check("NaN penalty rejected", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + { + const double w[] = {1}, x_ref[] = {0}; + qp.penalty = w; + elastiqp_setup(&work, &qp, NULL); + check("update_penalty: bad penalty rejected", + elastiqp_update_penalty(work, bad) == ELASTIQP_ERR_INVALID_ARG); + check("update_penalty: rejected update changed nothing", + elastiqp_solve(&res, work) == ELASTIQP_SOLVED && + near(x, x_ref, 1, 1e-10)); + elastiqp_free(work); + } + qp.n = 0; + check("n = 0 rejected", + elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + check("NULL workspace rejected", + elastiqp_solve(&res, NULL) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_update_q(NULL, x) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_free(NULL); + check("status names", + strcmp(elastiqp_status_name(ELASTIQP_SOLVED), "solved") == 0 && + strcmp(elastiqp_status_name(ELASTIQP_INFEASIBLE), + "infeasible") == 0 && + strcmp(elastiqp_status_name(ELASTIQP_ERR_NOT_SOLVED), + "not_solved") == 0); +} + +int main(void) { + printf("C interface\n"); + one_shot(); + elastic_and_hard(); + workspace(); + row_major(); + q_upper_triangle(); + inf_convention(); + row_updates(); + relax(); + options(); + invalid_args(); + printf("%s\n", g_all_ok ? "ALL OK" : "FAILURES"); + return g_all_ok ? 0 : 1; +} diff --git a/tests/test_c_parity.cc b/tests/test_c_parity.cc new file mode 100644 index 0000000..2965f78 --- /dev/null +++ b/tests/test_c_parity.cc @@ -0,0 +1,245 @@ +// C/C++ parity for solves, updates, warm starts, and relaxation. + +#include +#include +#include +#include + +#include "elastiqp/elastiqp.hpp" +#include "elastiqp/elastiqp_c.h" +#include "problem_gen.hpp" +#include "test_util.hpp" + +using Eigen::MatrixXd; +using Eigen::VectorXd; +using problem_gen::QPData; +using test_util::Check; +namespace das = elastiqp::das; + +namespace { + +using RowMajor = + Eigen::Matrix; + +// Row-major copies of the problem data, kept alive for the C calls. +struct CData { + RowMajor Q, A, G; + VectorXd q, b, h, w; + ElastiQPProblem Problem() const { + return {static_cast(q.size()), + static_cast(b.size()), + static_cast(h.size()), + Q.data(), + q.data(), + A.data(), + b.data(), + G.data(), + h.data(), + w.data()}; + } +}; + +// Ensure both APIs read a bitwise-identical symmetric Q. +void Symmetrize(MatrixXd& Q) { Q = (0.5 * (Q + Q.transpose())).eval(); } + +CData ToC(const QPData& qp, const VectorXd& w) { + return {qp.Q, qp.A, qp.G, qp.q, qp.b, qp.h, w}; +} + +struct CSolution { + VectorXd x, y, z, t, z_t; + ElastiQPResult res{}; + CSolution(int n, int m, int p) : x(n), y(m), z(p), t(p), z_t(p) { + res.x = x.data(); + res.y = y.data(); + res.z = z.data(); + res.t = t.data(); + res.z_t = z_t.data(); + } +}; + +// Largest difference between the C and C++ results (0 when identical). +double Diff(const CSolution& c, const elastiqp::Solution& s) { + double d = 0.0; + d = std::max(d, test_util::InfNorm(c.x - s.x)); + d = std::max(d, test_util::InfNorm(c.y - s.y)); + d = std::max(d, test_util::InfNorm(c.z - s.z)); + d = std::max(d, test_util::InfNorm(c.t - s.t)); + for (Eigen::Index i = 0; i < s.z_t.size(); ++i) { // inf on hard rows + if (c.z_t[i] != s.z_t[i]) d = std::max(d, std::abs(c.z_t[i] - s.z_t[i])); + } + const bool same_meta = std::strcmp(elastiqp_status_name(c.res.exitflag), + elastiqp::status_name(s.status)) == 0 && + c.res.iter == s.iters && + c.res.outer_iter == s.outer_iters && + c.res.n_active == s.n_active && + c.res.n_saturated == s.n_saturated && + c.res.fval == s.primal_obj; + return same_meta ? d : std::max(d, 1.0); +} + +void OneShot(std::mt19937& rng) { + double worst = 0.0; + for (int trial = 0; trial < 20; ++trial) { + const int n = 5 + trial % 7, m = trial % 3, p = 3 * n; + QPData qp = problem_gen::InfeasibleEq(rng, n, m, p, trial % 4); + Symmetrize(qp.Q); + VectorXd w = VectorXd::Constant(p, 5.0); + if (trial % 5 == 0) w[0] = std::numeric_limits::infinity(); + das::Settings st; + st.eps_abs = 1e-9; + const elastiqp::Solution ref = + das::Solve(qp.Q, qp.q, qp.A, qp.b, qp.G, qp.h, w, st); + + ElastiQPSettings cs; + elastiqp_default_settings(&cs); + cs.eps_abs = 1e-9; + const CData cd = ToC(qp, w); + const ElastiQPProblem cp = cd.Problem(); + CSolution sol(n, m, p); + elastiqp_quadprog(&sol.res, &cp, &cs); + worst = std::max(worst, Diff(sol, ref)); + } + Check("quadprog matches das::Solve (20 problems)", worst == 0.0, worst, + "diff"); +} + +void Chain(std::mt19937& rng) { + const int n = 12, m = 3, p = 40; + QPData qp = problem_gen::InfeasibleEq(rng, n, m, p, 4); + Symmetrize(qp.Q); + VectorXd w = VectorXd::Constant(p, 10.0); + das::Solver ref; + ref.setup(qp.Q, qp.q, qp.A, qp.b, qp.G, qp.h, w); + + CData cd = ToC(qp, w); + const ElastiQPProblem cp = cd.Problem(); + ElastiQPWorkspace* work = nullptr; + const int ret = elastiqp_setup(&work, &cp, nullptr); + CSolution sol(n, m, p); + + double worst = 0.0; + int solved = 0; + for (int tick = 0; tick < 60 && ret == 0; ++tick) { + qp.q += 0.05 * problem_gen::Randn(rng, n, 1); + qp.h += 0.05 * problem_gen::Randn(rng, p, 1); + qp.b += 0.01 * problem_gen::Randn(rng, m, 1); + cd.q = qp.q; + cd.h = qp.h; + cd.b = qp.b; + ref.set_q(qp.q); + ref.set_h(qp.h); + ref.set_b(qp.b); + elastiqp_update_q(work, cd.q.data()); + elastiqp_update_h(work, cd.h.data()); + elastiqp_update_b(work, cd.b.data()); + if (tick % 4 == 1) { // a few rows of G and A move, in full + qp.G.row(tick % p) += 0.1 * problem_gen::Randn(rng, 1, n); + qp.A.row(tick % m) += 0.01 * problem_gen::Randn(rng, 1, n); + cd.G = qp.G; + cd.A = qp.A; + ref.set_G(qp.G); + ref.set_A(qp.A); + elastiqp_update_G(work, cd.G.data()); + elastiqp_update_A(work, cd.A.data()); + } + if (tick % 4 == 3) { // a block of G rows, through update_G_rows + const int first = tick % (p - 3); + qp.G.middleRows(first, 3) += 0.1 * problem_gen::Randn(rng, 3, n); + cd.G = qp.G; + ref.set_G(qp.G); + elastiqp_update_G_rows(work, first, 3, cd.G.row(first).data()); + } + if (tick % 10 == 3) { + const MatrixXd dQ = 0.05 * problem_gen::Randn(rng, n, n); + qp.Q += dQ * dQ.transpose(); + Symmetrize(qp.Q); + cd.Q = qp.Q; + ref.set_Q(qp.Q); + elastiqp_update_Q(work, cd.Q.data()); + } + if (tick % 10 == 7) { + w = VectorXd::Constant(p, 2.0 + tick % 3); + cd.w = w; + ref.set_penalty(w); + elastiqp_update_penalty(work, cd.w.data()); + } + if (tick % 15 == 9) { // explicit warm start from the previous solution + const VectorXd x = sol.x, y = sol.y, z = sol.z; + ref.set_warm_start(x, y, z); + elastiqp_set_warm_start(work, x.data(), y.data(), z.data()); + } + const elastiqp::Solution& r = ref.solve(); + elastiqp_solve(&sol.res, work); + worst = std::max(worst, Diff(sol, r)); + solved += sol.res.exitflag == ELASTIQP_SOLVED; + } + Check("setup", ret == 0, ret, "ret"); + Check("update chain matches das::Solver (60 ticks)", worst == 0.0, worst, + "diff"); + Check("update chain solved", solved == 60, solved, "solved"); + + const elastiqp::Solution& r = ref.relax(1e-4); + elastiqp_relax(&sol.res, work, 1e-4, 1e-6, 50); + Check("relax matches", Diff(sol, r) == 0.0, Diff(sol, r), "diff"); + elastiqp_free(work); +} + +// Check defaults and named access for every settings field. +void SettingsTable() { + const das::Settings d; + ElastiQPSettings c; + elastiqp_default_settings(&c); + bool ok = true; + int fields = 0; + ElastiQPWorkspace* work = nullptr; + const double Q[] = {1}, q[] = {0}; + const ElastiQPProblem qp = {1, 0, 0, Q, q, + nullptr, nullptr, nullptr, nullptr, nullptr}; + elastiqp_setup(&work, &qp, nullptr); + const auto field = [&](const char* name, double c_value, double cpp_value) { + double v = -1; + ok &= c_value == cpp_value; + ok &= elastiqp_get_option(work, name, &v) == 0 && v == cpp_value; + ++fields; + }; +#define FIELD(f) field(#f, c.f, d.f) + FIELD(eps_abs); + FIELD(eps_rel); + FIELD(sing_tol); + FIELD(zero_tol); + FIELD(eps_prox); + FIELD(eta_prox); + FIELD(prox_relaxation); + FIELD(prox_escalations); + FIELD(max_iter); + FIELD(max_outer); + FIELD(warm_start); + FIELD(reuse_factorization); + FIELD(ruiz); + FIELD(ruiz_max_iter); + FIELD(ruiz_tol); + FIELD(ruiz_refresh_ratio); + FIELD(check_eq_consistency); + FIELD(progress_tol); + FIELD(cycle_tol); + FIELD(refactor_tol); + FIELD(relax_reg); + FIELD(relax_factor_retries); +#undef FIELD + elastiqp_free(work); + Check("option table covers every settings field", + ok && fields == elastiqp_num_options(), fields, "fields"); +} + +} // namespace + +int main() { + std::printf("C interface vs C++ API\n"); + std::mt19937 rng(7); + SettingsTable(); + OneShot(rng); + Chain(rng); + std::printf("%s\n", test_util::g_all_ok ? "ALL OK" : "FAILURES"); + return test_util::g_all_ok ? 0 : 1; +} From 9df93f6ae073ac0a193c771410e606345d4fe5da Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Fri, 9 Oct 2026 08:28:59 -0700 Subject: [PATCH 2/8] Update naming conventions to better align with the C++ Co-authored-by: Joris Gillis --- README.md | 10 +-- include/elastiqp/elastiqp_c.h | 70 +++++++++--------- src/elastiqp_c.cc | 62 ++++++++-------- tests/test_c_api.c | 134 +++++++++++++++++----------------- tests/test_c_parity.cc | 39 +++++----- 5 files changed, 156 insertions(+), 159 deletions(-) diff --git a/README.md b/README.md index 5527130..7a7b7a4 100644 --- a/README.md +++ b/README.md @@ -106,7 +106,7 @@ while (running) { ### C -The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Exit flags: `> 0` success, `0` unsolved, `< 0` failure. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. +The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Status codes: `> 0` success, `0` unsolved, `< 0` failure. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. ```c #include "elastiqp/elastiqp_c.h" @@ -114,15 +114,15 @@ The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP // 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}; -ElastiQPResult res = {.x = x, .z = z}; // caller-owned outputs; NULL ones are skipped -elastiqp_quadprog(&res, &qp, NULL); // NULL: default settings +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_update_q(work, q_k); elastiqp_update_h(work, h_k); elastiqp_update_b(work, b_k); - elastiqp_update_G_rows(work, first, count, rows_k); // changed rows + elastiqp_set_q(work, q_k); elastiqp_set_h(work, h_k); elastiqp_set_b(work, b_k); + elastiqp_set_G_rows(work, first, count, rows_k); // changed rows if (elastiqp_solve(&res, work) < 0) { /* ... */ } } elastiqp_free(work); diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h index c16b9d7..3f57e1a 100644 --- a/include/elastiqp/elastiqp_c.h +++ b/include/elastiqp/elastiqp_c.h @@ -11,11 +11,11 @@ * Q is symmetric; only its upper triangle is read. * Two-sided constraints and variable bounds use two rows of G (g and -g). * Penalties must be > 0; values >= ELASTIQP_INF make the row hard (t = 0). - * Exit flags: > 0 success, 0 unsolved, < 0 failure. + * Status codes: > 0 success, 0 unsolved, < 0 failure. * - * Use elastiqp_quadprog() for one-shot solves, or setup/update/solve/free + * Use elastiqp_solve_once() for one-shot solves, or setup/set_X/solve/free * for repeated solves. Input data is copied; caller buffers can be reused. - * Allocation-free update/solve coverage: tests/test_das_alloc.cc. + * Allocation-free set_X/solve coverage: tests/test_das_alloc.cc. * * Windows DLL users outside CMake: define ELASTIQP_C_SHARED before including. */ @@ -42,7 +42,7 @@ extern "C" { /* Hard-row penalty threshold. */ #define ELASTIQP_INF 1e30 -/* Exit flags. */ +/* Status codes. */ #define ELASTIQP_SOLVED 1 #define ELASTIQP_UNSOLVED 0 /* no solve yet */ #define ELASTIQP_MAX_ITER (-1) /* iteration limit reached */ @@ -109,7 +109,7 @@ typedef struct { * Caller-owned output arrays; NULL arrays are skipped. * KKT: Q x + q + A'y + G'z = 0, 0 <= z <= penalty, t = max(G x - h, 0). * z_t = penalty - z after solve; relax() computes it independently. - * On API errors, only exitflag is valid; statistics are zeroed. + * On API errors, only status is valid; statistics are zeroed. */ typedef struct { double* x; /* n */ @@ -118,32 +118,32 @@ typedef struct { double* t; /* p, elastic slacks */ double* z_t; /* p, duals of t >= 0 (INFINITY on hard rows) */ - int exitflag; /* ELASTIQP_SOLVED, ... */ - int iter; /* inner active-set (or relax Newton) iterations */ - int outer_iter; /* proximal rounds */ + int status; /* ELASTIQP_SOLVED, ... */ + int iters; /* inner active-set (or relax Newton) iterations */ + int outer_iters; /* proximal rounds */ int n_active; /* rows with t = 0 and z > 0 */ int n_saturated; /* rows with t > 0 (z at its penalty) */ - double fval; /* primal objective, penalty term included */ + double primal_obj; /* penalty term included */ double primal_res; double dual_res; double duality_gap; - double setup_time; /* seconds; elastiqp_quadprog only */ + double setup_time; /* seconds; elastiqp_solve_once only */ double solve_time; /* seconds */ -} ElastiQPResult; +} ElastiQPSolution; /* Opaque solver workspace. */ typedef struct ElastiQPWorkspace ElastiQPWorkspace; ELASTIQP_C_API void elastiqp_default_settings(ElastiQPSettings* settings); -/* Name of an exit flag ("solved", "max_iter", ..., "invalid_arg"). */ -ELASTIQP_C_API const char* elastiqp_status_name(int exitflag); +/* Name of a status code ("solved", "max_iter", ..., "invalid_arg"). */ +ELASTIQP_C_API const char* elastiqp_status_name(int status); -/* One-shot solve; settings may be NULL (defaults). Returns res->exitflag. */ -ELASTIQP_C_API int elastiqp_quadprog(ElastiQPResult* res, - const ElastiQPProblem* qp, - const ElastiQPSettings* settings); +/* One-shot solve; settings may be NULL (defaults). Returns res->status. */ +ELASTIQP_C_API int elastiqp_solve_once(ElastiQPSolution* res, + const ElastiQPProblem* qp, + const ElastiQPSettings* settings); /* * Creates a workspace for qp (settings may be NULL). Returns 0 and sets *work, @@ -154,8 +154,9 @@ ELASTIQP_C_API int elastiqp_setup(ElastiQPWorkspace** work, const ElastiQPSettings* settings); /* Solves the current problem, warm-started from the previous solve when - * settings.warm_start is set. Returns res->exitflag. */ -ELASTIQP_C_API int elastiqp_solve(ElastiQPResult* res, ElastiQPWorkspace* work); + * settings.warm_start is set. Returns res->status. */ +ELASTIQP_C_API int elastiqp_solve(ElastiQPSolution* res, + ElastiQPWorkspace* work); ELASTIQP_C_API void elastiqp_free(ElastiQPWorkspace* work); @@ -197,23 +198,23 @@ ELASTIQP_C_API int elastiqp_option_info(int index, const char** name, * Zero-sized updates do nothing. Returns 0 or a negative error; * rejected updates leave the workspace unchanged. */ -ELASTIQP_C_API int elastiqp_update_Q(ElastiQPWorkspace* work, const double* Q); -ELASTIQP_C_API int elastiqp_update_q(ElastiQPWorkspace* work, const double* q); -ELASTIQP_C_API int elastiqp_update_A(ElastiQPWorkspace* work, const double* A); -ELASTIQP_C_API int elastiqp_update_b(ElastiQPWorkspace* work, const double* b); -ELASTIQP_C_API int elastiqp_update_G(ElastiQPWorkspace* work, const double* G); -ELASTIQP_C_API int elastiqp_update_h(ElastiQPWorkspace* work, const double* h); -ELASTIQP_C_API int elastiqp_update_penalty(ElastiQPWorkspace* work, - const double* penalty); +ELASTIQP_C_API int elastiqp_set_Q(ElastiQPWorkspace* work, const double* Q); +ELASTIQP_C_API int elastiqp_set_q(ElastiQPWorkspace* work, const double* q); +ELASTIQP_C_API int elastiqp_set_A(ElastiQPWorkspace* work, const double* A); +ELASTIQP_C_API int elastiqp_set_b(ElastiQPWorkspace* work, const double* b); +ELASTIQP_C_API int elastiqp_set_G(ElastiQPWorkspace* work, const double* G); +ELASTIQP_C_API int elastiqp_set_h(ElastiQPWorkspace* work, const double* h); +ELASTIQP_C_API int elastiqp_set_penalty(ElastiQPWorkspace* work, + const double* penalty); /* * Replace rows [first, first + count) with a count x n row-major block. * Cost: O(count * n). */ -ELASTIQP_C_API int elastiqp_update_A_rows(ElastiQPWorkspace* work, int first, - int count, const double* rows); -ELASTIQP_C_API int elastiqp_update_G_rows(ElastiQPWorkspace* work, int first, - int count, const double* rows); +ELASTIQP_C_API int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, + int count, const double* rows); +ELASTIQP_C_API int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, + int count, const double* rows); /* * Seed the next solve with (x, y, z). Requires x; NULL y or z means zero. @@ -226,10 +227,11 @@ ELASTIQP_C_API int elastiqp_set_warm_start(ElastiQPWorkspace* work, * Compute the kappa-relaxed central point for differentiation (kappa > 0). * Requires a successful solve on current data, else ELASTIQP_ERR_NOT_SOLVED. * Preserves the working set and warm start; p == 0 returns the solution. - * Suggested: tol = 1e-6, max_iter = 50. Returns res->exitflag. + * Suggested: tol = 1e-6, max_iter = 50. Returns res->status. */ -ELASTIQP_C_API int elastiqp_relax(ElastiQPResult* res, ElastiQPWorkspace* work, - double kappa, double tol, int max_iter); +ELASTIQP_C_API int elastiqp_relax(ElastiQPSolution* res, + ElastiQPWorkspace* work, double kappa, + double tol, int max_iter); #ifdef __cplusplus } diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index 8029e0e..4644a0f 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -35,7 +35,7 @@ double Seconds(Clock::time_point since) { return std::chrono::duration(Clock::now() - since).count(); } -int ExitFlag(elastiqp::Status s) { +int StatusCode(elastiqp::Status s) { switch (s) { case elastiqp::Status::kSolved: return ELASTIQP_SOLVED; @@ -243,32 +243,32 @@ void StagePenalty(VectorXd& dst, const double* w) { } // Clear statistics on API errors. -int Fail(ElastiQPResult* res, int flag) { - res->exitflag = flag; - res->iter = res->outer_iter = res->n_active = res->n_saturated = 0; - res->fval = res->primal_res = res->dual_res = res->duality_gap = 0.0; +int Fail(ElastiQPSolution* res, int flag) { + res->status = flag; + res->iters = res->outer_iters = res->n_active = res->n_saturated = 0; + res->primal_obj = res->primal_res = res->dual_res = res->duality_gap = 0.0; res->setup_time = res->solve_time = 0.0; return flag; } -// Fills the caller's buffers from sol; returns the exit flag. -int Extract(ElastiQPResult* res, const elastiqp::Solution& sol) { +// Fills the caller's buffers from sol; returns the status code. +int Extract(ElastiQPSolution* res, const elastiqp::Solution& sol) { using Eigen::Map; if (res->x) Map(res->x, sol.x.size()) = sol.x; if (res->y) Map(res->y, sol.y.size()) = sol.y; if (res->z) Map(res->z, sol.z.size()) = sol.z; if (res->t) Map(res->t, sol.t.size()) = sol.t; if (res->z_t) Map(res->z_t, sol.z_t.size()) = sol.z_t; - res->exitflag = ExitFlag(sol.status); - res->iter = sol.iters; - res->outer_iter = sol.outer_iters; + res->status = StatusCode(sol.status); + res->iters = sol.iters; + res->outer_iters = sol.outer_iters; res->n_active = sol.n_active; res->n_saturated = sol.n_saturated; - res->fval = sol.primal_obj; + res->primal_obj = sol.primal_obj; res->primal_res = sol.primal_res; res->dual_res = sol.dual_res; res->duality_gap = sol.duality_gap; - return res->exitflag; + return res->status; } // Runs f, mapping C++ exceptions to error codes. @@ -328,8 +328,8 @@ void elastiqp_default_settings(ElastiQPSettings* settings) { if (settings) *settings = ToC(das::Settings{}); } -const char* elastiqp_status_name(int exitflag) { - switch (exitflag) { +const char* elastiqp_status_name(int status) { + switch (status) { case ELASTIQP_SOLVED: return "solved"; case ELASTIQP_UNSOLVED: @@ -352,8 +352,8 @@ const char* elastiqp_status_name(int exitflag) { return "?"; } -int elastiqp_quadprog(ElastiQPResult* res, const ElastiQPProblem* qp, - const ElastiQPSettings* settings) { +int elastiqp_solve_once(ElastiQPSolution* res, const ElastiQPProblem* qp, + const ElastiQPSettings* settings) { if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; if (!Valid(qp) || (settings && !ValidSettings(*settings))) { return Fail(res, ELASTIQP_ERR_INVALID_ARG); @@ -393,7 +393,7 @@ int elastiqp_setup(ElastiQPWorkspace** work, const ElastiQPProblem* qp, return 0; } -int elastiqp_solve(ElastiQPResult* res, ElastiQPWorkspace* work) { +int elastiqp_solve(ElastiQPSolution* res, ElastiQPWorkspace* work) { if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; if (work == nullptr) return Fail(res, ELASTIQP_ERR_INVALID_ARG); double solve_time = 0.0; @@ -470,7 +470,7 @@ int elastiqp_option_info(int index, const char** name, int* type, return 0; } -int elastiqp_update_Q(ElastiQPWorkspace* work, const double* Q) { +int elastiqp_set_Q(ElastiQPWorkspace* work, const double* Q) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; StageQ(work->Q, Q, work->solver.n()); work->solver.set_Q(work->Q); @@ -478,7 +478,7 @@ int elastiqp_update_Q(ElastiQPWorkspace* work, const double* Q) { return 0; } -int elastiqp_update_q(ElastiQPWorkspace* work, const double* q) { +int elastiqp_set_q(ElastiQPWorkspace* work, const double* q) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; if (q) { work->q = VectorMap(q, work->solver.n()); @@ -490,12 +490,12 @@ int elastiqp_update_q(ElastiQPWorkspace* work, const double* q) { return 0; } -int elastiqp_update_A(ElastiQPWorkspace* work, const double* A) { +int elastiqp_set_A(ElastiQPWorkspace* work, const double* A) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; - return elastiqp_update_A_rows(work, 0, work->solver.m(), A); + return elastiqp_set_A_rows(work, 0, work->solver.m(), A); } -int elastiqp_update_b(ElastiQPWorkspace* work, const double* b) { +int elastiqp_set_b(ElastiQPWorkspace* work, const double* b) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; const int m = work->solver.m(); if (m == 0) return 0; @@ -506,12 +506,12 @@ int elastiqp_update_b(ElastiQPWorkspace* work, const double* b) { return 0; } -int elastiqp_update_G(ElastiQPWorkspace* work, const double* G) { +int elastiqp_set_G(ElastiQPWorkspace* work, const double* G) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; - return elastiqp_update_G_rows(work, 0, work->solver.p(), G); + return elastiqp_set_G_rows(work, 0, work->solver.p(), G); } -int elastiqp_update_h(ElastiQPWorkspace* work, const double* h) { +int elastiqp_set_h(ElastiQPWorkspace* work, const double* h) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; const int p = work->solver.p(); if (p == 0) return 0; @@ -522,7 +522,7 @@ int elastiqp_update_h(ElastiQPWorkspace* work, const double* h) { return 0; } -int elastiqp_update_penalty(ElastiQPWorkspace* work, const double* penalty) { +int elastiqp_set_penalty(ElastiQPWorkspace* work, const double* penalty) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; const int p = work->solver.p(); if (p == 0) return 0; @@ -535,8 +535,8 @@ int elastiqp_update_penalty(ElastiQPWorkspace* work, const double* penalty) { return 0; } -int elastiqp_update_A_rows(ElastiQPWorkspace* work, int first, int count, - const double* rows) { +int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, int count, + const double* rows) { if (work == nullptr || !ValidRows(first, count, work->solver.m(), rows)) { return ELASTIQP_ERR_INVALID_ARG; } @@ -546,8 +546,8 @@ int elastiqp_update_A_rows(ElastiQPWorkspace* work, int first, int count, return 0; } -int elastiqp_update_G_rows(ElastiQPWorkspace* work, int first, int count, - const double* rows) { +int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, int count, + const double* rows) { if (work == nullptr || !ValidRows(first, count, work->solver.p(), rows)) { return ELASTIQP_ERR_INVALID_ARG; } @@ -576,7 +576,7 @@ int elastiqp_set_warm_start(ElastiQPWorkspace* work, const double* x, return 0; } -int elastiqp_relax(ElastiQPResult* res, ElastiQPWorkspace* work, double kappa, +int elastiqp_relax(ElastiQPSolution* res, ElastiQPWorkspace* work, double kappa, double tol, int max_iter) { if (res == nullptr) return ELASTIQP_ERR_INVALID_ARG; if (work == nullptr || !(kappa > 0.0) || !(tol > 0.0) || max_iter < 0) { diff --git a/tests/test_c_api.c b/tests/test_c_api.c index 97c84b3..a644a67 100644 --- a/tests/test_c_api.c +++ b/tests/test_c_api.c @@ -35,18 +35,18 @@ static void one_shot(void) { const double x_ref[] = {0.5, 0.5}, z_ref[] = {0.5}, t_ref[] = {0}; double x[2], z[1], t[1]; ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; - ElastiQPResult res = {x, NULL, z, t}; + ElastiQPSolution res = {x, NULL, z, t}; ElastiQPSettings s = tight(); - int flag = elastiqp_quadprog(&res, &qp, &s); - check("one-shot: solved", flag == ELASTIQP_SOLVED && res.exitflag == flag); + int flag = elastiqp_solve_once(&res, &qp, &s); + check("one-shot: solved", flag == ELASTIQP_SOLVED && res.status == flag); check("one-shot: x, z, t", near(x, x_ref, 2, 1e-8) && near(z, z_ref, 1, 1e-8) && near(t, t_ref, 1, 1e-12)); check("one-shot: active row counted", res.n_active == 1 && res.n_saturated == 0); - check("one-shot: objective", fabs(res.fval + 0.75) < 1e-8); + check("one-shot: objective", fabs(res.primal_obj + 0.75) < 1e-8); check("one-shot: default settings (NULL)", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_SOLVED && + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-5)); } @@ -61,16 +61,16 @@ static void elastic_and_hard(void) { const double t_ref[] = {1, 1}, z_ref[] = {0.5, 0.5}; const double t_hard[] = {2, 0}; ElastiQPProblem qp = {1, 0, 2, Q, NULL, NULL, NULL, G, h, w}; - ElastiQPResult res = {x, NULL, z, t}; + ElastiQPSolution res = {x, NULL, z, t}; ElastiQPSettings s = tight(); - int flag = elastiqp_quadprog(&res, &qp, &s); + int flag = elastiqp_solve_once(&res, &qp, &s); check("conflict: solved", flag == ELASTIQP_SOLVED); check("conflict: x = 0, t = 1, z at penalty", fabs(x[0]) < 1e-8 && near(t, t_ref, 2, 1e-8) && near(z, z_ref, 2, 1e-8)); check("conflict: both rows saturated", res.n_saturated == 2); w[1] = INFINITY; - flag = elastiqp_quadprog(&res, &qp, &s); + flag = elastiqp_solve_once(&res, &qp, &s); check("hard row: solved", flag == ELASTIQP_SOLVED); check("hard row: x = 1, hard slack 0", fabs(x[0] - 1) < 1e-8 && near(t, t_hard, 2, 1e-8)); @@ -84,7 +84,7 @@ static void workspace(void) { double h[] = {0.8}; double x[2], y[1], z[1], t[1]; ElastiQPProblem qp = {2, 1, 1, Q, q, A, b, G, h, w}; - ElastiQPResult res = {x, y, z, t}; + ElastiQPSolution res = {x, y, z, t}; ElastiQPSettings s = tight(); ElastiQPWorkspace* work = NULL; int n = 0, m = 0, p = 0, flag; @@ -103,41 +103,40 @@ static void workspace(void) { h[0] = 2.0; /* workspace retains its copy */ elastiqp_solve(&res, work); check("data copied at setup", near(x, x1, 2, 1e-8)); - check("update_h", elastiqp_update_h(work, h) == 0); + check("set_h", elastiqp_set_h(work, h) == 0); elastiqp_solve(&res, work); - check("update_h: constraint released", + check("set_h: constraint released", near(x, x2, 2, 1e-8) && near(z, z2, 1, 1e-8)); { const double b2[] = {-0.5}; - elastiqp_update_b(work, b2); + elastiqp_set_b(work, b2); elastiqp_solve(&res, work); - check("update_b", res.exitflag == ELASTIQP_SOLVED && near(x, x3, 2, 1e-8)); + check("set_b", res.status == ELASTIQP_SOLVED && near(x, x3, 2, 1e-8)); } { const double q2[] = {-1, -3}, x4[] = {-0.5, 2.0}; - elastiqp_update_q(work, q2); + elastiqp_set_q(work, q2); elastiqp_solve(&res, work); - check("update_q (row active at h = 2)", near(x, x4, 2, 1e-8)); + check("set_q (row active at h = 2)", near(x, x4, 2, 1e-8)); } { /* Penalty 0.5 < marginal cost 1 at x2 = 2: the row gives way. */ const double w2[] = {0.5}, x5[] = {-0.5, 2.5}, t5[] = {0.5}; - elastiqp_update_penalty(work, w2); + elastiqp_set_penalty(work, w2); elastiqp_solve(&res, work); - check("update_penalty: row saturates", - near(x, x5, 2, 1e-8) && near(t, t5, 1, 1e-8) && - res.n_saturated == 1); + check("set_penalty: row saturates", + near(x, x5, 2, 1e-8) && near(t, t5, 1, 1e-8) && res.n_saturated == 1); } { const double Q2[] = {2, 0, 0, 2}, A2[] = {0, 1}, G2[] = {1, 0}; const double x6[] = {0.5, -0.5}; /* min |x|^2 - x1 - 3 x2 s.t. x2 = -0.5, x1 <= 2 (inactive). */ - elastiqp_update_Q(work, Q2); - elastiqp_update_A(work, A2); - elastiqp_update_G(work, G2); + elastiqp_set_Q(work, Q2); + elastiqp_set_A(work, A2); + elastiqp_set_G(work, G2); elastiqp_solve(&res, work); - check("update_Q/A/G", near(x, x6, 2, 1e-8)); + check("set_Q/A/G", near(x, x6, 2, 1e-8)); } { ElastiQPSettings got; @@ -149,10 +148,10 @@ static void workspace(void) { elastiqp_set_warm_start(work, xw, NULL, NULL); elastiqp_solve(&res, work); check("set_settings: max_iter = 0 fails", - res.exitflag == ELASTIQP_MAX_ITER && res.exitflag < 0); + res.status == ELASTIQP_MAX_ITER && res.status < 0); elastiqp_set_settings(work, &s); elastiqp_solve(&res, work); - check("set_settings: restored", res.exitflag == ELASTIQP_SOLVED); + check("set_settings: restored", res.status == ELASTIQP_SOLVED); } elastiqp_free(work); } @@ -169,11 +168,11 @@ static void row_major(void) { const double x_ref[] = {1 - 6.0 / 14, 1 - 12.0 / 14, 1 - 18.0 / 14}; double x[3]; ElastiQPProblem qp = {3, 0, 2, Q, q, NULL, NULL, G, h, w}; - ElastiQPResult res = {x, NULL, NULL, NULL}; + ElastiQPSolution res = {x, NULL, NULL, NULL}; ElastiQPSettings s = tight(); - elastiqp_quadprog(&res, &qp, &s); - check("row-major G", res.exitflag == ELASTIQP_SOLVED && - near(x, x_ref, 3, 1e-8)); + elastiqp_solve_once(&res, &qp, &s); + check("row-major G", + res.status == ELASTIQP_SOLVED && near(x, x_ref, 3, 1e-8)); } /* Only the upper triangle of Q is read: full and upper-only Q agree. */ @@ -182,12 +181,12 @@ static void q_upper_triangle(void) { const double q[] = {-1, 0}; double x1[2], x2[2]; ElastiQPProblem qp = {2, 0, 0, Qfull, q, NULL, NULL, NULL, NULL, NULL}; - ElastiQPResult r1 = {x1}, r2 = {x2}; + ElastiQPSolution r1 = {x1}, r2 = {x2}; /* 2 x1 + x2 = 1, x1 + 2 x2 = 0 */ const double x_ref[] = {2.0 / 3, -1.0 / 3}; - elastiqp_quadprog(&r1, &qp, NULL); + elastiqp_solve_once(&r1, &qp, NULL); qp.Q = Qupper; - elastiqp_quadprog(&r2, &qp, NULL); + elastiqp_solve_once(&r2, &qp, NULL); check("Q: full matrix", near(x1, x_ref, 2, 1e-10)); check("Q: upper triangle only, lower ignored", near(x2, x_ref, 2, 1e-10)); } @@ -198,10 +197,10 @@ static void inf_convention(void) { const double w[] = {0.5, ELASTIQP_INF}; double x[1], t[2], z_t[2]; ElastiQPProblem qp = {1, 0, 2, Q, NULL, NULL, NULL, G, h, w}; - ElastiQPResult res = {x, NULL, NULL, t, z_t}; - elastiqp_quadprog(&res, &qp, NULL); + ElastiQPSolution res = {x, NULL, NULL, t, z_t}; + elastiqp_solve_once(&res, &qp, NULL); check("ELASTIQP_INF penalty is a hard row", - res.exitflag == ELASTIQP_SOLVED && fabs(x[0] - 1) < 1e-8 && + res.status == ELASTIQP_SOLVED && fabs(x[0] - 1) < 1e-8 && fabs(t[1]) < 1e-12 && isinf(z_t[1])); } @@ -213,25 +212,23 @@ static void row_updates(void) { const double rows[] = {0.5, 1, 2, 0.5}; double x[2], x_ref[2]; ElastiQPProblem qp = {2, 0, 3, Q, q, NULL, NULL, G, h, w}; - ElastiQPResult res = {x}, ref = {x_ref}; + ElastiQPSolution res = {x}, ref = {x_ref}; ElastiQPWorkspace* work = NULL; elastiqp_setup(&work, &qp, NULL); elastiqp_solve(&res, work); - check("update_G_rows: rows 1..2", - elastiqp_update_G_rows(work, 1, 2, rows) == 0); + check("set_G_rows: rows 1..2", elastiqp_set_G_rows(work, 1, 2, rows) == 0); G[2] = 0.5, G[3] = 1, G[4] = 2, G[5] = 0.5; elastiqp_solve(&res, work); - elastiqp_quadprog(&ref, &qp, NULL); - check("update_G_rows: matches a fresh solve", - res.exitflag == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-8)); - check("update_G_rows: out of range rejected", - elastiqp_update_G_rows(work, 2, 2, rows) == ELASTIQP_ERR_INVALID_ARG && - elastiqp_update_G_rows(work, -1, 1, rows) == - ELASTIQP_ERR_INVALID_ARG && - elastiqp_update_A_rows(work, 0, 1, rows) == - ELASTIQP_ERR_INVALID_ARG); - check("update_G_rows: count 0 is a no-op", - elastiqp_update_G_rows(work, 3, 0, NULL) == 0); + elastiqp_solve_once(&ref, &qp, NULL); + check("set_G_rows: matches a fresh solve", + res.status == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-8)); + check( + "set_G_rows: out of range rejected", + elastiqp_set_G_rows(work, 2, 2, rows) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_G_rows(work, -1, 1, rows) == ELASTIQP_ERR_INVALID_ARG && + elastiqp_set_A_rows(work, 0, 1, rows) == ELASTIQP_ERR_INVALID_ARG); + check("set_G_rows: count 0 is a no-op", + elastiqp_set_G_rows(work, 3, 0, NULL) == 0); elastiqp_free(work); } @@ -240,7 +237,7 @@ static void relax(void) { const double G[] = {1, 1}, h[] = {1}, w[] = {10}; double x[2], z[1], t[1], z_t[1]; ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; - ElastiQPResult res = {x, NULL, z, t, z_t}; + ElastiQPSolution res = {x, NULL, z, t, z_t}; ElastiQPWorkspace* work = NULL; elastiqp_setup(&work, &qp, NULL); check("relax before solve: not_solved", @@ -248,15 +245,14 @@ static void relax(void) { elastiqp_solve(&res, work); check("solve: z_t = penalty - z", fabs(z_t[0] - (10 - z[0])) < 1e-12); check("relax: kappa <= 0 rejected", - elastiqp_relax(&res, work, 0.0, 1e-8, 50) == - ELASTIQP_ERR_INVALID_ARG && - res.iter == 0); + elastiqp_relax(&res, work, 0.0, 1e-8, 50) == ELASTIQP_ERR_INVALID_ARG && + res.iters == 0); elastiqp_relax(&res, work, 1e-3, 1e-8, 50); /* At the relaxed point the inequality slack s = h - Gx + t is O(kappa). */ check("relax: solved, interior point", - res.exitflag == ELASTIQP_SOLVED && 1 - x[0] - x[1] + t[0] > 0 && + res.status == ELASTIQP_SOLVED && 1 - x[0] - x[1] + t[0] > 0 && z[0] > 0 && z[0] < 10 && z_t[0] > 0); - elastiqp_update_q(work, q2); + elastiqp_set_q(work, q2); check("relax after an update: not_solved", elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_ERR_NOT_SOLVED); elastiqp_solve(&res, work); @@ -271,7 +267,7 @@ static void options(void) { const double G[] = {1, 1}, h[] = {1}, w[] = {10}; double x[2], v = 0, def = 0, lo = 0, hi = 0; ElastiQPProblem qp = {2, 0, 1, Q, q, NULL, NULL, G, h, w}; - ElastiQPResult res = {x}; + ElastiQPSolution res = {x}; ElastiQPSettings s; ElastiQPWorkspace* work = NULL; const char *name = NULL, *name2 = NULL, *desc = NULL; @@ -341,7 +337,7 @@ static void options(void) { elastiqp_free(work); check("setup: out-of-range settings rejected", elastiqp_setup(&work, &qp, &s) == ELASTIQP_ERR_INVALID_ARG && - elastiqp_quadprog(&res, &qp, &s) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_solve_once(&res, &qp, &s) == ELASTIQP_ERR_INVALID_ARG); elastiqp_default_settings(&s); s.warm_start = 5; /* struct booleans are 0 / nonzero */ elastiqp_setup(&work, &qp, &s); @@ -356,41 +352,41 @@ static void invalid_args(void) { double x[1]; double bad[] = {0}; ElastiQPProblem qp = {1, 0, 1, Q, NULL, NULL, NULL, G, h, NULL}; - ElastiQPResult res = {x}; + ElastiQPSolution res = {x}; ElastiQPWorkspace* work = (ElastiQPWorkspace*)&qp; - res.iter = 7; + res.iters = 7; check("missing penalty rejected, stats cleared", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG && - res.exitflag == ELASTIQP_ERR_INVALID_ARG && res.iter == 0); + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG && + res.status == ELASTIQP_ERR_INVALID_ARG && res.iters == 0); check("setup: rejects and clears *work", elastiqp_setup(&work, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG && work == NULL); qp.penalty = bad; check("penalty 0 rejected", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); bad[0] = -1; check("negative penalty rejected", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); bad[0] = NAN; check("NaN penalty rejected", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); { const double w[] = {1}, x_ref[] = {0}; qp.penalty = w; elastiqp_setup(&work, &qp, NULL); - check("update_penalty: bad penalty rejected", - elastiqp_update_penalty(work, bad) == ELASTIQP_ERR_INVALID_ARG); - check("update_penalty: rejected update changed nothing", + check("set_penalty: bad penalty rejected", + elastiqp_set_penalty(work, bad) == ELASTIQP_ERR_INVALID_ARG); + check("set_penalty: rejected update changed nothing", elastiqp_solve(&res, work) == ELASTIQP_SOLVED && near(x, x_ref, 1, 1e-10)); elastiqp_free(work); } qp.n = 0; check("n = 0 rejected", - elastiqp_quadprog(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); check("NULL workspace rejected", elastiqp_solve(&res, NULL) == ELASTIQP_ERR_INVALID_ARG && - elastiqp_update_q(NULL, x) == ELASTIQP_ERR_INVALID_ARG); + elastiqp_set_q(NULL, x) == ELASTIQP_ERR_INVALID_ARG); elastiqp_free(NULL); check("status names", strcmp(elastiqp_status_name(ELASTIQP_SOLVED), "solved") == 0 && diff --git a/tests/test_c_parity.cc b/tests/test_c_parity.cc index 2965f78..0c60505 100644 --- a/tests/test_c_parity.cc +++ b/tests/test_c_parity.cc @@ -48,7 +48,7 @@ CData ToC(const QPData& qp, const VectorXd& w) { struct CSolution { VectorXd x, y, z, t, z_t; - ElastiQPResult res{}; + ElastiQPSolution res{}; CSolution(int n, int m, int p) : x(n), y(m), z(p), t(p), z_t(p) { res.x = x.data(); res.y = y.data(); @@ -68,13 +68,12 @@ double Diff(const CSolution& c, const elastiqp::Solution& s) { for (Eigen::Index i = 0; i < s.z_t.size(); ++i) { // inf on hard rows if (c.z_t[i] != s.z_t[i]) d = std::max(d, std::abs(c.z_t[i] - s.z_t[i])); } - const bool same_meta = std::strcmp(elastiqp_status_name(c.res.exitflag), - elastiqp::status_name(s.status)) == 0 && - c.res.iter == s.iters && - c.res.outer_iter == s.outer_iters && - c.res.n_active == s.n_active && - c.res.n_saturated == s.n_saturated && - c.res.fval == s.primal_obj; + const bool same_meta = + std::strcmp(elastiqp_status_name(c.res.status), + elastiqp::status_name(s.status)) == 0 && + c.res.iters == s.iters && c.res.outer_iters == s.outer_iters && + c.res.n_active == s.n_active && c.res.n_saturated == s.n_saturated && + c.res.primal_obj == s.primal_obj; return same_meta ? d : std::max(d, 1.0); } @@ -97,10 +96,10 @@ void OneShot(std::mt19937& rng) { const CData cd = ToC(qp, w); const ElastiQPProblem cp = cd.Problem(); CSolution sol(n, m, p); - elastiqp_quadprog(&sol.res, &cp, &cs); + elastiqp_solve_once(&sol.res, &cp, &cs); worst = std::max(worst, Diff(sol, ref)); } - Check("quadprog matches das::Solve (20 problems)", worst == 0.0, worst, + Check("solve_once matches das::Solve (20 problems)", worst == 0.0, worst, "diff"); } @@ -130,9 +129,9 @@ void Chain(std::mt19937& rng) { ref.set_q(qp.q); ref.set_h(qp.h); ref.set_b(qp.b); - elastiqp_update_q(work, cd.q.data()); - elastiqp_update_h(work, cd.h.data()); - elastiqp_update_b(work, cd.b.data()); + elastiqp_set_q(work, cd.q.data()); + elastiqp_set_h(work, cd.h.data()); + elastiqp_set_b(work, cd.b.data()); if (tick % 4 == 1) { // a few rows of G and A move, in full qp.G.row(tick % p) += 0.1 * problem_gen::Randn(rng, 1, n); qp.A.row(tick % m) += 0.01 * problem_gen::Randn(rng, 1, n); @@ -140,15 +139,15 @@ void Chain(std::mt19937& rng) { cd.A = qp.A; ref.set_G(qp.G); ref.set_A(qp.A); - elastiqp_update_G(work, cd.G.data()); - elastiqp_update_A(work, cd.A.data()); + elastiqp_set_G(work, cd.G.data()); + elastiqp_set_A(work, cd.A.data()); } - if (tick % 4 == 3) { // a block of G rows, through update_G_rows + if (tick % 4 == 3) { // a block of G rows, through set_G_rows const int first = tick % (p - 3); qp.G.middleRows(first, 3) += 0.1 * problem_gen::Randn(rng, 3, n); cd.G = qp.G; ref.set_G(qp.G); - elastiqp_update_G_rows(work, first, 3, cd.G.row(first).data()); + elastiqp_set_G_rows(work, first, 3, cd.G.row(first).data()); } if (tick % 10 == 3) { const MatrixXd dQ = 0.05 * problem_gen::Randn(rng, n, n); @@ -156,13 +155,13 @@ void Chain(std::mt19937& rng) { Symmetrize(qp.Q); cd.Q = qp.Q; ref.set_Q(qp.Q); - elastiqp_update_Q(work, cd.Q.data()); + elastiqp_set_Q(work, cd.Q.data()); } if (tick % 10 == 7) { w = VectorXd::Constant(p, 2.0 + tick % 3); cd.w = w; ref.set_penalty(w); - elastiqp_update_penalty(work, cd.w.data()); + elastiqp_set_penalty(work, cd.w.data()); } if (tick % 15 == 9) { // explicit warm start from the previous solution const VectorXd x = sol.x, y = sol.y, z = sol.z; @@ -172,7 +171,7 @@ void Chain(std::mt19937& rng) { const elastiqp::Solution& r = ref.solve(); elastiqp_solve(&sol.res, work); worst = std::max(worst, Diff(sol, r)); - solved += sol.res.exitflag == ELASTIQP_SOLVED; + solved += sol.res.status == ELASTIQP_SOLVED; } Check("setup", ret == 0, ret, "ret"); Check("update chain matches das::Solver (60 ticks)", worst == 0.0, worst, From bcccfd0d3448d8f7b564b450e65c174a2edd08cb Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Fri, 9 Oct 2026 08:44:06 -0700 Subject: [PATCH 3/8] Align status code values between c/c++ Co-authored-by: Joris Gillis --- README.md | 4 ++-- include/elastiqp/elastiqp_c.h | 22 ++++++++++++---------- src/elastiqp_c.cc | 30 +++++++++++------------------- tests/test_c_api.c | 2 +- 4 files changed, 26 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index 7a7b7a4..1729411 100644 --- a/README.md +++ b/README.md @@ -106,7 +106,7 @@ while (running) { ### C -The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Status codes: `> 0` success, `0` unsolved, `< 0` failure. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. +The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Status codes match `elastiqp::Status` (`ELASTIQP_SOLVED`, `ELASTIQP_MAX_ITER`, ...); API errors are negative. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. ```c #include "elastiqp/elastiqp_c.h" @@ -123,7 +123,7 @@ 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_rows(work, first, count, rows_k); // changed rows - if (elastiqp_solve(&res, work) < 0) { /* ... */ } + if (elastiqp_solve(&res, work) != ELASTIQP_SOLVED) { /* ... */ } } elastiqp_free(work); ``` diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h index 3f57e1a..af4d420 100644 --- a/include/elastiqp/elastiqp_c.h +++ b/include/elastiqp/elastiqp_c.h @@ -11,7 +11,7 @@ * Q is symmetric; only its upper triangle is read. * Two-sided constraints and variable bounds use two rows of G (g and -g). * Penalties must be > 0; values >= ELASTIQP_INF make the row hard (t = 0). - * Status codes: > 0 success, 0 unsolved, < 0 failure. + * Solver status codes match elastiqp::Status (0-4); API errors are < 0. * * Use elastiqp_solve_once() for one-shot solves, or setup/set_X/solve/free * for repeated solves. Input data is copied; caller buffers can be reused. @@ -42,16 +42,18 @@ extern "C" { /* Hard-row penalty threshold. */ #define ELASTIQP_INF 1e30 -/* Status codes. */ +/* Solver status codes, equal to elastiqp::Status. */ +#define ELASTIQP_UNSOLVED 0 /* no solve yet */ #define ELASTIQP_SOLVED 1 -#define ELASTIQP_UNSOLVED 0 /* no solve yet */ -#define ELASTIQP_MAX_ITER (-1) /* iteration limit reached */ -#define ELASTIQP_NUMERICS (-2) /* numerical failure */ -#define ELASTIQP_INFEASIBLE (-3) /* inconsistent equalities or hard rows */ -#define ELASTIQP_ERR_INVALID_ARG (-10) /* NULL, bad dimension, value, name */ -#define ELASTIQP_ERR_NOT_SOLVED (-11) /* relax() without a current solve */ -#define ELASTIQP_ERR_ALLOC (-12) /* out of memory */ -#define ELASTIQP_ERR_INTERNAL (-13) /* unexpected C++ exception */ +#define ELASTIQP_MAX_ITER 2 /* iteration limit reached */ +#define ELASTIQP_NUMERICS 3 /* numerical failure */ +#define ELASTIQP_INFEASIBLE 4 /* inconsistent equalities or hard rows */ + +/* API errors. */ +#define ELASTIQP_ERR_INVALID_ARG (-1) /* NULL, bad dimension, value, name */ +#define ELASTIQP_ERR_NOT_SOLVED (-2) /* relax() without a current solve */ +#define ELASTIQP_ERR_ALLOC (-3) /* out of memory */ +#define ELASTIQP_ERR_INTERNAL (-4) /* unexpected C++ exception */ typedef struct { int n; /* number of variables, >= 1 */ diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index 4644a0f..9687c95 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -35,21 +35,13 @@ double Seconds(Clock::time_point since) { return std::chrono::duration(Clock::now() - since).count(); } -int StatusCode(elastiqp::Status s) { - switch (s) { - case elastiqp::Status::kSolved: - return ELASTIQP_SOLVED; - case elastiqp::Status::kUnsolved: - return ELASTIQP_UNSOLVED; - case elastiqp::Status::kMaxIter: - return ELASTIQP_MAX_ITER; - case elastiqp::Status::kNumerics: - return ELASTIQP_NUMERICS; - case elastiqp::Status::kInfeasible: - return ELASTIQP_INFEASIBLE; - } - return ELASTIQP_ERR_INTERNAL; -} +static_assert( + ELASTIQP_UNSOLVED == static_cast(elastiqp::Status::kUnsolved) && + ELASTIQP_SOLVED == static_cast(elastiqp::Status::kSolved) && + ELASTIQP_MAX_ITER == static_cast(elastiqp::Status::kMaxIter) && + ELASTIQP_NUMERICS == static_cast(elastiqp::Status::kNumerics) && + ELASTIQP_INFEASIBLE == static_cast(elastiqp::Status::kInfeasible), + "C status codes must equal elastiqp::Status"); // Shared metadata for settings conversion, validation, and named options. // Add new fields here and in tests/test_c_parity.cc. @@ -259,7 +251,7 @@ int Extract(ElastiQPSolution* res, const elastiqp::Solution& sol) { if (res->z) Map(res->z, sol.z.size()) = sol.z; if (res->t) Map(res->t, sol.t.size()) = sol.t; if (res->z_t) Map(res->z_t, sol.z_t.size()) = sol.z_t; - res->status = StatusCode(sol.status); + res->status = static_cast(sol.status); res->iters = sol.iters; res->outer_iters = sol.outer_iters; res->n_active = sol.n_active; @@ -369,7 +361,7 @@ int elastiqp_solve_once(ElastiQPSolution* res, const ElastiQPProblem* qp, solve_time = Seconds(t1); return Extract(res, sol); }); - if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + if (flag < 0) return Fail(res, flag); res->setup_time = setup_time; res->solve_time = solve_time; return flag; @@ -404,7 +396,7 @@ int elastiqp_solve(ElastiQPSolution* res, ElastiQPWorkspace* work) { return Extract(res, sol); }); work->solved = flag == ELASTIQP_SOLVED; - if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + if (flag < 0) return Fail(res, flag); res->setup_time = 0.0; res->solve_time = solve_time; return flag; @@ -590,7 +582,7 @@ int elastiqp_relax(ElastiQPSolution* res, ElastiQPWorkspace* work, double kappa, solve_time = Seconds(t0); return Extract(res, sol); }); - if (flag <= ELASTIQP_ERR_INVALID_ARG) return Fail(res, flag); + if (flag < 0) return Fail(res, flag); res->setup_time = 0.0; res->solve_time = solve_time; return flag; diff --git a/tests/test_c_api.c b/tests/test_c_api.c index a644a67..8a57ee9 100644 --- a/tests/test_c_api.c +++ b/tests/test_c_api.c @@ -148,7 +148,7 @@ static void workspace(void) { elastiqp_set_warm_start(work, xw, NULL, NULL); elastiqp_solve(&res, work); check("set_settings: max_iter = 0 fails", - res.status == ELASTIQP_MAX_ITER && res.status < 0); + res.status == ELASTIQP_MAX_ITER); elastiqp_set_settings(work, &s); elastiqp_solve(&res, work); check("set_settings: restored", res.status == ELASTIQP_SOLVED); From 4f901d2bab6c9bc9c959ef659184c52cb6eb660f Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Fri, 9 Oct 2026 14:31:03 -0700 Subject: [PATCH 4/8] Add can_relax helper in das cpp + use in C interface --- include/elastiqp/das.hpp | 23 +++++++++++++++++------ src/elastiqp_c.cc | 12 +----------- tests/test_c_api.c | 11 +++++++++++ 3 files changed, 29 insertions(+), 17 deletions(-) diff --git a/include/elastiqp/das.hpp b/include/elastiqp/das.hpp index 81159ae..406a629 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -219,6 +219,17 @@ class Solver { 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(m_ + i)]; } @@ -355,13 +366,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; diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index 9687c95..e83b3a0 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -23,8 +23,6 @@ struct ElastiQPWorkspace { das::Solver solver; MatrixXd Q; VectorXd q, b, h, penalty, x, y, z; - // Successful solve on current data; required by relax(). - bool solved = false; }; namespace { @@ -395,7 +393,6 @@ int elastiqp_solve(ElastiQPSolution* res, ElastiQPWorkspace* work) { solve_time = Seconds(t0); return Extract(res, sol); }); - work->solved = flag == ELASTIQP_SOLVED; if (flag < 0) return Fail(res, flag); res->setup_time = 0.0; res->solve_time = solve_time; @@ -466,7 +463,6 @@ int elastiqp_set_Q(ElastiQPWorkspace* work, const double* Q) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; StageQ(work->Q, Q, work->solver.n()); work->solver.set_Q(work->Q); - work->solved = false; return 0; } @@ -478,7 +474,6 @@ int elastiqp_set_q(ElastiQPWorkspace* work, const double* q) { work->q.setZero(); } work->solver.set_q(work->q); - work->solved = false; return 0; } @@ -494,7 +489,6 @@ int elastiqp_set_b(ElastiQPWorkspace* work, const double* b) { if (b == nullptr) return ELASTIQP_ERR_INVALID_ARG; work->b = VectorMap(b, m); work->solver.set_b(work->b); - work->solved = false; return 0; } @@ -510,7 +504,6 @@ int elastiqp_set_h(ElastiQPWorkspace* work, const double* h) { if (h == nullptr) return ELASTIQP_ERR_INVALID_ARG; work->h = VectorMap(h, p); work->solver.set_h(work->h); - work->solved = false; return 0; } @@ -523,7 +516,6 @@ int elastiqp_set_penalty(ElastiQPWorkspace* work, const double* penalty) { } StagePenalty(work->penalty, penalty); work->solver.set_penalty(work->penalty); - work->solved = false; return 0; } @@ -534,7 +526,6 @@ int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, int count, } if (count == 0) return 0; work->solver.set_A_rows(first, RowMajorMap(rows, count, work->solver.n())); - work->solved = false; return 0; } @@ -545,7 +536,6 @@ int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, int count, } if (count == 0) return 0; work->solver.set_G_rows(first, RowMajorMap(rows, count, work->solver.n())); - work->solved = false; return 0; } @@ -574,7 +564,7 @@ int elastiqp_relax(ElastiQPSolution* res, ElastiQPWorkspace* work, double kappa, if (work == nullptr || !(kappa > 0.0) || !(tol > 0.0) || max_iter < 0) { return Fail(res, ELASTIQP_ERR_INVALID_ARG); } - if (!work->solved) return Fail(res, ELASTIQP_ERR_NOT_SOLVED); + if (!work->solver.can_relax()) return Fail(res, ELASTIQP_ERR_NOT_SOLVED); double solve_time = 0.0; const int flag = Guard([&] { const Clock::time_point t0 = Clock::now(); diff --git a/tests/test_c_api.c b/tests/test_c_api.c index 8a57ee9..3c6f329 100644 --- a/tests/test_c_api.c +++ b/tests/test_c_api.c @@ -258,6 +258,17 @@ static void relax(void) { elastiqp_solve(&res, work); check("relax after re-solve", elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_SOLVED); + elastiqp_set_G(work, G); + check("relax after re-sending identical rows", + elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_SOLVED); + ElastiQPSettings s; + elastiqp_get_settings(work, &s); + s.max_iter = 0; + elastiqp_set_settings(work, &s); + elastiqp_set_q(work, q); + check("solve hits max_iter", elastiqp_solve(&res, work) == ELASTIQP_MAX_ITER); + check("relax after a failed solve: not_solved", + elastiqp_relax(&res, work, 1e-3, 1e-8, 50) == ELASTIQP_ERR_NOT_SOLVED); elastiqp_free(work); } From 2e7c87389a23bc394ce900ee082aa72ccebceb1a Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Mon, 12 Oct 2026 21:48:37 -0700 Subject: [PATCH 5/8] Column-major C interface --- README.md | 2 -- include/elastiqp/elastiqp_c.h | 10 +++++----- src/elastiqp_c.cc | 16 +++++++--------- tests/test_c_api.c | 28 +++++++++++++++++----------- tests/test_c_parity.cc | 10 ++++------ 5 files changed, 33 insertions(+), 33 deletions(-) diff --git a/README.md b/README.md index 1729411..1dcb14b 100644 --- a/README.md +++ b/README.md @@ -106,8 +106,6 @@ while (running) { ### C -The C interface uses the active-set backend. Build `elastiqp_c` with `-DELASTIQP_BUILD_C=ON` (default for top-level builds); add `-DBUILD_SHARED_LIBS=ON` for a shared library. Matrices are dense row-major; only the upper triangle of `Q` is read. Status codes match `elastiqp::Status` (`ELASTIQP_SOLVED`, `ELASTIQP_MAX_ITER`, ...); API errors are negative. See [`elastiqp_c.h`](include/elastiqp/elastiqp_c.h) for the API. - ```c #include "elastiqp/elastiqp_c.h" diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h index af4d420..75fa81c 100644 --- a/include/elastiqp/elastiqp_c.h +++ b/include/elastiqp/elastiqp_c.h @@ -7,8 +7,8 @@ * subject to A x = b * G x - t <= h, t >= 0 * - * Dense row-major matrices: M(i, j) = M[i * cols + j]. - * Q is symmetric; only its upper triangle is read. + * Dense column-major matrices: M(i, j) = M[i + j * rows]. + * Q is symmetric; only its upper triangle (i <= j) is read. * Two-sided constraints and variable bounds use two rows of G (g and -g). * Penalties must be > 0; values >= ELASTIQP_INF make the row hard (t = 0). * Solver status codes match elastiqp::Status (0-4); API errors are < 0. @@ -62,9 +62,9 @@ typedef struct { const double* Q; /* n x n, upper triangle read; NULL: zero (LP) */ const double* q; /* n; NULL: zero */ - const double* A; /* m x n, row-major */ + const double* A; /* m x n, column-major */ const double* b; /* m */ - const double* G; /* p x n, row-major */ + const double* G; /* p x n, column-major */ const double* h; /* p */ const double* penalty; /* p, > 0; >= ELASTIQP_INF: hard row */ } ElastiQPProblem; @@ -210,7 +210,7 @@ ELASTIQP_C_API int elastiqp_set_penalty(ElastiQPWorkspace* work, const double* penalty); /* - * Replace rows [first, first + count) with a count x n row-major block. + * Replace rows [first, first + count) with a count x n column-major block. * Cost: O(count * n). */ ELASTIQP_C_API int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index e83b3a0..bf343f4 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -12,9 +12,7 @@ using Eigen::MatrixXd; using Eigen::VectorXd; -using RowMajor = Eigen::Matrix; -using RowMajorMap = Eigen::Map; +using MatrixMap = Eigen::Map; using VectorMap = Eigen::Map; namespace das = elastiqp::das; @@ -215,10 +213,10 @@ bool Valid(const ElastiQPProblem* qp) { return true; } -// Upper triangle of the row-major Q, mirrored; zero for NULL. +// Upper triangle of the column-major Q, mirrored; zero for NULL. void StageQ(MatrixXd& dst, const double* Q, int n) { if (Q) { - dst = RowMajorMap(Q, n, n).selfadjointView(); + dst = MatrixMap(Q, n, n).selfadjointView(); } else { dst.setZero(); } @@ -286,13 +284,13 @@ int SetupImpl(ElastiQPWorkspace& w, const ElastiQPProblem& qp, MatrixXd A(m, n), G(p, n); w.b.resize(m); if (m > 0) { - A = RowMajorMap(qp.A, m, n); + A = MatrixMap(qp.A, m, n); w.b = VectorMap(qp.b, m); } w.h.resize(p); w.penalty.resize(p); if (p > 0) { - G = RowMajorMap(qp.G, p, n); + G = MatrixMap(qp.G, p, n); w.h = VectorMap(qp.h, p); StagePenalty(w.penalty, qp.penalty); } @@ -525,7 +523,7 @@ int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, int count, return ELASTIQP_ERR_INVALID_ARG; } if (count == 0) return 0; - work->solver.set_A_rows(first, RowMajorMap(rows, count, work->solver.n())); + work->solver.set_A_rows(first, MatrixMap(rows, count, work->solver.n())); return 0; } @@ -535,7 +533,7 @@ int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, int count, return ELASTIQP_ERR_INVALID_ARG; } if (count == 0) return 0; - work->solver.set_G_rows(first, RowMajorMap(rows, count, work->solver.n())); + work->solver.set_G_rows(first, MatrixMap(rows, count, work->solver.n())); return 0; } diff --git a/tests/test_c_api.c b/tests/test_c_api.c index 3c6f329..b73886c 100644 --- a/tests/test_c_api.c +++ b/tests/test_c_api.c @@ -157,13 +157,13 @@ static void workspace(void) { } /* - * Row-major layout: G is 2 x 3. min 0.5|x|^2 + q'x with + * Column-major layout: G is 2 x 3. min 0.5|x|^2 + q'x with * G = [[1, 2, 3], [0, 0, 1]], h = [0, 10] (second row inactive). The * minimizer of 0.5|x - c|^2 over g'x <= 0 is c - (g'c / g'g) g. */ -static void row_major(void) { +static void column_major(void) { const double Q[] = {1, 0, 0, 0, 1, 0, 0, 0, 1}, q[] = {-1, -1, -1}; - const double G[] = {1, 2, 3, 0, 0, 1}, h[] = {0, 10}, w[] = {100, 100}; + const double G[] = {1, 0, 2, 0, 3, 1}, h[] = {0, 10}, w[] = {100, 100}; /* c = (1,1,1), g'c = 6, g'g = 14 */ const double x_ref[] = {1 - 6.0 / 14, 1 - 12.0 / 14, 1 - 18.0 / 14}; double x[3]; @@ -171,13 +171,16 @@ static void row_major(void) { ElastiQPSolution res = {x, NULL, NULL, NULL}; ElastiQPSettings s = tight(); elastiqp_solve_once(&res, &qp, &s); - check("row-major G", + check("column-major G", res.status == ELASTIQP_SOLVED && near(x, x_ref, 3, 1e-8)); } -/* Only the upper triangle of Q is read: full and upper-only Q agree. */ +/* + * Only the upper triangle of Q is read: full and upper-only Q agree. Column- + * major, Q(0, 1) is Q[2]; reading the lower triangle would pick up the -99. + */ static void q_upper_triangle(void) { - const double Qfull[] = {2, 1, 1, 2}, Qupper[] = {2, 1, -99, 2}; + const double Qfull[] = {2, 1, 1, 2}, Qupper[] = {2, -99, 1, 2}; const double q[] = {-1, 0}; double x1[2], x2[2]; ElastiQPProblem qp = {2, 0, 0, Qfull, q, NULL, NULL, NULL, NULL, NULL}; @@ -204,12 +207,15 @@ static void inf_convention(void) { fabs(t[1]) < 1e-12 && isinf(z_t[1])); } -/* Rows of a 3 x 2 G replaced in blocks, against fresh one-shot solves. */ +/* + * Rows of a 3 x 2 G replaced in blocks, against fresh one-shot solves. + * G = [[1, 0], [0, 1], [1, 1]]; rows 1..2 become [[0.5, 1], [2, 0.5]]. + */ static void row_updates(void) { const double Q[] = {1, 0, 0, 1}, q[] = {-2, -2}, h[] = {1, 1, 1}; const double w[] = {10, 10, 10}; - double G[] = {1, 0, 0, 1, 1, 1}; - const double rows[] = {0.5, 1, 2, 0.5}; + double G[] = {1, 0, 1, 0, 1, 1}; + const double rows[] = {0.5, 2, 1, 0.5}; double x[2], x_ref[2]; ElastiQPProblem qp = {2, 0, 3, Q, q, NULL, NULL, G, h, w}; ElastiQPSolution res = {x}, ref = {x_ref}; @@ -217,7 +223,7 @@ static void row_updates(void) { elastiqp_setup(&work, &qp, NULL); elastiqp_solve(&res, work); check("set_G_rows: rows 1..2", elastiqp_set_G_rows(work, 1, 2, rows) == 0); - G[2] = 0.5, G[3] = 1, G[4] = 2, G[5] = 0.5; + G[1] = 0.5, G[2] = 2, G[4] = 1, G[5] = 0.5; elastiqp_solve(&res, work); elastiqp_solve_once(&ref, &qp, NULL); check("set_G_rows: matches a fresh solve", @@ -412,7 +418,7 @@ int main(void) { one_shot(); elastic_and_hard(); workspace(); - row_major(); + column_major(); q_upper_triangle(); inf_convention(); row_updates(); diff --git a/tests/test_c_parity.cc b/tests/test_c_parity.cc index 0c60505..f19ca09 100644 --- a/tests/test_c_parity.cc +++ b/tests/test_c_parity.cc @@ -18,12 +18,9 @@ namespace das = elastiqp::das; namespace { -using RowMajor = - Eigen::Matrix; - -// Row-major copies of the problem data, kept alive for the C calls. +// Copies of the problem data, kept alive for the C calls. struct CData { - RowMajor Q, A, G; + MatrixXd Q, A, G; VectorXd q, b, h, w; ElastiQPProblem Problem() const { return {static_cast(q.size()), @@ -145,9 +142,10 @@ void Chain(std::mt19937& rng) { if (tick % 4 == 3) { // a block of G rows, through set_G_rows const int first = tick % (p - 3); qp.G.middleRows(first, 3) += 0.1 * problem_gen::Randn(rng, 3, n); + const MatrixXd block = qp.G.middleRows(first, 3); cd.G = qp.G; ref.set_G(qp.G); - elastiqp_set_G_rows(work, first, 3, cd.G.row(first).data()); + elastiqp_set_G_rows(work, first, 3, block.data()); } if (tick % 10 == 3) { const MatrixXd dQ = 0.05 * problem_gen::Randn(rng, n, n); From 41b9b0ec9551e2258f8dc6bcdc03462fb9d1fec2 Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Mon, 12 Oct 2026 22:38:40 -0700 Subject: [PATCH 6/8] Use Refs to prevent extra copies of solver data (helps simplify the C interface) --- include/elastiqp/das.hpp | 37 ++++++++++++----------- src/elastiqp_c.cc | 65 ++++++++++++---------------------------- tests/test_das_alloc.cc | 37 +++++++++++++++++++++++ 3 files changed, 76 insertions(+), 63 deletions(-) diff --git a/include/elastiqp/das.hpp b/include/elastiqp/das.hpp index 406a629..e1ef7c4 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -16,6 +16,9 @@ namespace elastiqp::das { +using MatRef = Eigen::Ref; +using VecRef = Eigen::Ref; + struct Settings { double eps_abs = 1e-6; double eps_rel = 0.0; @@ -76,9 +79,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(q.size()); m_ = static_cast(b.size()); p_ = static_cast(h.size()); @@ -158,17 +161,17 @@ 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)); } @@ -176,20 +179,20 @@ class Solver { 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_G(const MatRef& G) { set_rows(G, m_); } // Replace rows [first, first + R.rows()); read R in place. template void set_A_rows(int first, const Eigen::MatrixBase& R) { @@ -199,11 +202,11 @@ class Solver { void set_G_rows(int first, const Eigen::MatrixBase& R) { set_rows(R, m_ + first); } - void set_h(const VectorXd& h) { + 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; } @@ -212,7 +215,7 @@ 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; diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index bf343f4..2513223 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -2,6 +2,7 @@ #include "elastiqp/elastiqp_c.h" +#include #include #include #include @@ -16,11 +17,10 @@ using MatrixMap = Eigen::Map; using VectorMap = Eigen::Map; namespace das = elastiqp::das; -// Preallocated staging buffers; A/G rows are read directly from caller data. struct ElastiQPWorkspace { das::Solver solver; MatrixXd Q; - VectorXd q, b, h, penalty, x, y, z; + VectorXd penalty, zero; }; namespace { @@ -271,34 +271,24 @@ int Guard(F&& f) { } } +// v, or zeros when v is NULL. +VectorMap OrZero(const ElastiQPWorkspace& w, const double* v, int size) { + return VectorMap(v ? v : w.zero.data(), size); +} + int SetupImpl(ElastiQPWorkspace& w, const ElastiQPProblem& qp, const ElastiQPSettings* settings) { const int n = qp.n, m = qp.m, p = qp.p; w.Q.resize(n, n); StageQ(w.Q, qp.Q, n); - if (qp.q) { - w.q = VectorMap(qp.q, n); - } else { - w.q.setZero(n); - } - MatrixXd A(m, n), G(p, n); - w.b.resize(m); - if (m > 0) { - A = MatrixMap(qp.A, m, n); - w.b = VectorMap(qp.b, m); - } - w.h.resize(p); w.penalty.resize(p); - if (p > 0) { - G = MatrixMap(qp.G, p, n); - w.h = VectorMap(qp.h, p); - StagePenalty(w.penalty, qp.penalty); - } - w.x.resize(n); - w.y.resize(m); - w.z.resize(p); + if (p > 0) StagePenalty(w.penalty, qp.penalty); + w.zero.setZero(std::max({n, m, p})); if (settings) w.solver.settings = ToCpp(*settings); - w.solver.setup(w.Q, w.q, A, w.b, G, w.h, w.penalty); + // NULL A/G/b/h only occur with m or p zero; Maps of size 0 never read them. + w.solver.setup(w.Q, OrZero(w, qp.q, n), MatrixMap(qp.A, m, n), + VectorMap(qp.b, m), MatrixMap(qp.G, p, n), VectorMap(qp.h, p), + w.penalty); return 0; } @@ -466,12 +456,7 @@ int elastiqp_set_Q(ElastiQPWorkspace* work, const double* Q) { int elastiqp_set_q(ElastiQPWorkspace* work, const double* q) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; - if (q) { - work->q = VectorMap(q, work->solver.n()); - } else { - work->q.setZero(); - } - work->solver.set_q(work->q); + work->solver.set_q(OrZero(*work, q, work->solver.n())); return 0; } @@ -485,8 +470,7 @@ int elastiqp_set_b(ElastiQPWorkspace* work, const double* b) { const int m = work->solver.m(); if (m == 0) return 0; if (b == nullptr) return ELASTIQP_ERR_INVALID_ARG; - work->b = VectorMap(b, m); - work->solver.set_b(work->b); + work->solver.set_b(VectorMap(b, m)); return 0; } @@ -500,8 +484,7 @@ int elastiqp_set_h(ElastiQPWorkspace* work, const double* h) { const int p = work->solver.p(); if (p == 0) return 0; if (h == nullptr) return ELASTIQP_ERR_INVALID_ARG; - work->h = VectorMap(h, p); - work->solver.set_h(work->h); + work->solver.set_h(VectorMap(h, p)); return 0; } @@ -540,19 +523,9 @@ int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, int count, int elastiqp_set_warm_start(ElastiQPWorkspace* work, const double* x, const double* y, const double* z) { if (work == nullptr || x == nullptr) return ELASTIQP_ERR_INVALID_ARG; - const int n = work->solver.n(), m = work->solver.m(), p = work->solver.p(); - work->x = VectorMap(x, n); - if (y) { - work->y = VectorMap(y, m); - } else { - work->y.setZero(); - } - if (z) { - work->z = VectorMap(z, p); - } else { - work->z.setZero(); - } - work->solver.set_warm_start(work->x, work->y, work->z); + work->solver.set_warm_start(VectorMap(x, work->solver.n()), + OrZero(*work, y, work->solver.m()), + OrZero(*work, z, work->solver.p())); return 0; } diff --git a/tests/test_das_alloc.cc b/tests/test_das_alloc.cc index 96ce686..62d738d 100644 --- a/tests/test_das_alloc.cc +++ b/tests/test_das_alloc.cc @@ -305,6 +305,42 @@ void ExplicitWarm(std::mt19937& rng) { r.Report("explicit warm start"); } +// Updates through Eigen::Map over caller buffers (as in the C interface) bind +// to the setters' Ref arguments without a copy. +void MappedUpdates(std::mt19937& rng) { + QPData qp = problem_gen::InfeasibleEq(rng, 20, 4, 40, 4); + VectorXd w = VectorXd::Constant(40, 10.0); + using MatMap = Eigen::Map; + using VecMap = Eigen::Map; + const int n = 20, m = 4, p = 40; + Run r; + r.s.setup(qp.Q, qp.q, qp.A, qp.b, qp.G, qp.h, w); + r.Solve(); + for (int tick = 0; tick < 6; ++tick) { + Perturb(rng, qp.q, 0.05); + Perturb(rng, qp.h, 0.05); + Perturb(rng, qp.b, 0.01); + Perturb(rng, qp.A, 0.01); + Perturb(rng, qp.G, 0.05); + w *= 1.1; + const elastiqp::Solution& prev = r.s.solution(); + r.Update([&] { + r.s.set_q(VecMap(qp.q.data(), n)); + r.s.set_b(VecMap(qp.b.data(), m)); + r.s.set_h(VecMap(qp.h.data(), p)); + r.s.set_A(MatMap(qp.A.data(), m, n)); + r.s.set_G(MatMap(qp.G.data(), p, n)); + r.s.set_penalty(VecMap(w.data(), p)); + if (tick % 2 == 1) { + r.s.set_warm_start(VecMap(prev.x.data(), n), VecMap(prev.y.data(), m), + VecMap(prev.z.data(), p)); + } + }); + r.Solve(); + } + r.Report("updates through Eigen::Map"); +} + // relax() sizes its workspace on the first call after setup(); later calls // (and the solves in between) are allocation-free. void Relax(std::mt19937& rng) { @@ -351,6 +387,7 @@ int main() { Proximal(rng); RuizRefresh(rng); ExplicitWarm(rng); + MappedUpdates(rng); Relax(rng); Unconstrained(rng); for (const auto& [site, k] : alloc_guard::Sites()) From 9181f9631043e8789c1ff6a120067efacfe4eb09 Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Mon, 12 Oct 2026 22:59:21 -0700 Subject: [PATCH 7/8] Remove A/G row setters, stick to just A/G setters felt like this was unnecessary complexity now that we have the column-major structure in the C api and the Ref inputs --- README.md | 2 +- include/elastiqp/das.hpp | 12 +---------- include/elastiqp/elastiqp_c.h | 9 --------- src/elastiqp_c.cc | 38 +++++++++-------------------------- tests/test_c_api.c | 25 +++++++++-------------- tests/test_c_parity.cc | 8 -------- 6 files changed, 22 insertions(+), 72 deletions(-) diff --git a/README.md b/README.md index 1dcb14b..33cc42b 100644 --- a/README.md +++ b/README.md @@ -120,7 +120,7 @@ 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_rows(work, first, count, rows_k); // changed rows + elastiqp_set_G(work, G_k); // only rows that differ are refreshed if (elastiqp_solve(&res, work) != ELASTIQP_SOLVED) { /* ... */ } } elastiqp_free(work); diff --git a/include/elastiqp/das.hpp b/include/elastiqp/das.hpp index e1ef7c4..93f678e 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -193,15 +193,6 @@ class Solver { rhs_dirty_ = true; } void set_G(const MatRef& G) { set_rows(G, m_); } - // Replace rows [first, first + R.rows()); read R in place. - template - void set_A_rows(int first, const Eigen::MatrixBase& R) { - set_rows(R, first); - } - template - void set_G_rows(int first, const Eigen::MatrixBase& R) { - set_rows(R, m_ + first); - } void set_h(const VecRef& h) { rhs_.tail(p_) = h; rhs_dirty_ = true; @@ -693,8 +684,7 @@ class Solver { } // Marks only rows that actually changed. - template - void set_rows(const Eigen::MatrixBase& 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); diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h index 75fa81c..a741346 100644 --- a/include/elastiqp/elastiqp_c.h +++ b/include/elastiqp/elastiqp_c.h @@ -209,15 +209,6 @@ ELASTIQP_C_API int elastiqp_set_h(ElastiQPWorkspace* work, const double* h); ELASTIQP_C_API int elastiqp_set_penalty(ElastiQPWorkspace* work, const double* penalty); -/* - * Replace rows [first, first + count) with a count x n column-major block. - * Cost: O(count * n). - */ -ELASTIQP_C_API int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, - int count, const double* rows); -ELASTIQP_C_API int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, - int count, const double* rows); - /* * Seed the next solve with (x, y, z). Requires x; NULL y or z means zero. */ diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc index 2513223..e356ea4 100644 --- a/src/elastiqp_c.cc +++ b/src/elastiqp_c.cc @@ -292,12 +292,6 @@ int SetupImpl(ElastiQPWorkspace& w, const ElastiQPProblem& qp, return 0; } -// Validates a [first, first + count) row block of a matrix with `rows` rows. -bool ValidRows(int first, int count, int rows, const double* data) { - return first >= 0 && count >= 0 && first <= rows - count && - (count == 0 || data != nullptr); -} - } // namespace extern "C" { @@ -462,7 +456,11 @@ int elastiqp_set_q(ElastiQPWorkspace* work, const double* q) { int elastiqp_set_A(ElastiQPWorkspace* work, const double* A) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; - return elastiqp_set_A_rows(work, 0, work->solver.m(), A); + const int m = work->solver.m(); + if (m == 0) return 0; + if (A == nullptr) return ELASTIQP_ERR_INVALID_ARG; + work->solver.set_A(MatrixMap(A, m, work->solver.n())); + return 0; } int elastiqp_set_b(ElastiQPWorkspace* work, const double* b) { @@ -476,7 +474,11 @@ int elastiqp_set_b(ElastiQPWorkspace* work, const double* b) { int elastiqp_set_G(ElastiQPWorkspace* work, const double* G) { if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; - return elastiqp_set_G_rows(work, 0, work->solver.p(), G); + const int p = work->solver.p(); + if (p == 0) return 0; + if (G == nullptr) return ELASTIQP_ERR_INVALID_ARG; + work->solver.set_G(MatrixMap(G, p, work->solver.n())); + return 0; } int elastiqp_set_h(ElastiQPWorkspace* work, const double* h) { @@ -500,26 +502,6 @@ int elastiqp_set_penalty(ElastiQPWorkspace* work, const double* penalty) { return 0; } -int elastiqp_set_A_rows(ElastiQPWorkspace* work, int first, int count, - const double* rows) { - if (work == nullptr || !ValidRows(first, count, work->solver.m(), rows)) { - return ELASTIQP_ERR_INVALID_ARG; - } - if (count == 0) return 0; - work->solver.set_A_rows(first, MatrixMap(rows, count, work->solver.n())); - return 0; -} - -int elastiqp_set_G_rows(ElastiQPWorkspace* work, int first, int count, - const double* rows) { - if (work == nullptr || !ValidRows(first, count, work->solver.p(), rows)) { - return ELASTIQP_ERR_INVALID_ARG; - } - if (count == 0) return 0; - work->solver.set_G_rows(first, MatrixMap(rows, count, work->solver.n())); - return 0; -} - int elastiqp_set_warm_start(ElastiQPWorkspace* work, const double* x, const double* y, const double* z) { if (work == nullptr || x == nullptr) return ELASTIQP_ERR_INVALID_ARG; diff --git a/tests/test_c_api.c b/tests/test_c_api.c index b73886c..cb9255a 100644 --- a/tests/test_c_api.c +++ b/tests/test_c_api.c @@ -208,33 +208,28 @@ static void inf_convention(void) { } /* - * Rows of a 3 x 2 G replaced in blocks, against fresh one-shot solves. - * G = [[1, 0], [0, 1], [1, 1]]; rows 1..2 become [[0.5, 1], [2, 0.5]]. + * set_G with one changed row, against a fresh one-shot solve. + * G = [[1, 0], [0, 1], [1, 1]]; row 1 becomes [2, 0.5]. */ -static void row_updates(void) { +static void set_G_partial(void) { const double Q[] = {1, 0, 0, 1}, q[] = {-2, -2}, h[] = {1, 1, 1}; const double w[] = {10, 10, 10}; double G[] = {1, 0, 1, 0, 1, 1}; - const double rows[] = {0.5, 2, 1, 0.5}; double x[2], x_ref[2]; ElastiQPProblem qp = {2, 0, 3, Q, q, NULL, NULL, G, h, w}; ElastiQPSolution res = {x}, ref = {x_ref}; ElastiQPWorkspace* work = NULL; elastiqp_setup(&work, &qp, NULL); elastiqp_solve(&res, work); - check("set_G_rows: rows 1..2", elastiqp_set_G_rows(work, 1, 2, rows) == 0); - G[1] = 0.5, G[2] = 2, G[4] = 1, G[5] = 0.5; + G[1] = 2, G[4] = 0.5; + check("set_G: one row changed", elastiqp_set_G(work, G) == 0); elastiqp_solve(&res, work); elastiqp_solve_once(&ref, &qp, NULL); - check("set_G_rows: matches a fresh solve", + check("set_G: matches a fresh solve", res.status == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-8)); - check( - "set_G_rows: out of range rejected", - elastiqp_set_G_rows(work, 2, 2, rows) == ELASTIQP_ERR_INVALID_ARG && - elastiqp_set_G_rows(work, -1, 1, rows) == ELASTIQP_ERR_INVALID_ARG && - elastiqp_set_A_rows(work, 0, 1, rows) == ELASTIQP_ERR_INVALID_ARG); - check("set_G_rows: count 0 is a no-op", - elastiqp_set_G_rows(work, 3, 0, NULL) == 0); + check("set_G: NULL rejected", + elastiqp_set_G(work, NULL) == ELASTIQP_ERR_INVALID_ARG); + check("set_A: m = 0 is a no-op", elastiqp_set_A(work, NULL) == 0); elastiqp_free(work); } @@ -421,7 +416,7 @@ int main(void) { column_major(); q_upper_triangle(); inf_convention(); - row_updates(); + set_G_partial(); relax(); options(); invalid_args(); diff --git a/tests/test_c_parity.cc b/tests/test_c_parity.cc index f19ca09..1136588 100644 --- a/tests/test_c_parity.cc +++ b/tests/test_c_parity.cc @@ -139,14 +139,6 @@ void Chain(std::mt19937& rng) { elastiqp_set_G(work, cd.G.data()); elastiqp_set_A(work, cd.A.data()); } - if (tick % 4 == 3) { // a block of G rows, through set_G_rows - const int first = tick % (p - 3); - qp.G.middleRows(first, 3) += 0.1 * problem_gen::Randn(rng, 3, n); - const MatrixXd block = qp.G.middleRows(first, 3); - cd.G = qp.G; - ref.set_G(qp.G); - elastiqp_set_G_rows(work, first, 3, block.data()); - } if (tick % 10 == 3) { const MatrixXd dQ = 0.05 * problem_gen::Randn(rng, n, n); qp.Q += dQ * dQ.transpose(); From d06226c56c90c110202b86fbb0f1186dfd035a29 Mon Sep 17 00:00:00 2001 From: Daniel Morton Date: Mon, 12 Oct 2026 23:10:01 -0700 Subject: [PATCH 8/8] Consistent use of MatRef/VecRef across backends --- include/elastiqp/common.hpp | 7 ++-- include/elastiqp/das.hpp | 3 -- include/elastiqp/ipm.hpp | 54 +++++++++++++++++------------- include/elastiqp/pdal.hpp | 67 +++++++++++++++++++++---------------- 4 files changed, 74 insertions(+), 57 deletions(-) diff --git a/include/elastiqp/common.hpp b/include/elastiqp/common.hpp index c9df705..799bf87 100644 --- a/include/elastiqp/common.hpp +++ b/include/elastiqp/common.hpp @@ -12,6 +12,9 @@ namespace elastiqp { using Eigen::MatrixXd; using Eigen::VectorXd; +using MatRef = Eigen::Ref; +using VecRef = Eigen::Ref; + 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(m_)); diff --git a/include/elastiqp/das.hpp b/include/elastiqp/das.hpp index 93f678e..eeed7f0 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -16,9 +16,6 @@ namespace elastiqp::das { -using MatRef = Eigen::Ref; -using VecRef = Eigen::Ref; - struct Settings { double eps_abs = 1e-6; double eps_rel = 0.0; diff --git a/include/elastiqp/ipm.hpp b/include/elastiqp/ipm.hpp index 689d023..6c667ff 100644 --- a/include/elastiqp/ipm.hpp +++ b/include/elastiqp/ipm.hpp @@ -54,9 +54,9 @@ class Solver { public: Settings settings; - 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_ = q.size(); m_ = b.size(); p_ = h.size(); @@ -138,19 +138,19 @@ class Solver { matrix_dirty_ = 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)); } @@ -161,34 +161,40 @@ class Solver { // Setters store data in the current Ruiz frame; solve() re-equilibrates // after a matrix change (every solve is cold, so there is no iterate to // carry across the new scaling). - void set_Q(const MatrixXd& Q) { + void set_Q(const MatRef& Q) { Q_ = 0.5 * (Q + Q.transpose()); if (ruiz_) Q_ = c_s_ * dx_s_.asDiagonal() * Q_ * dx_s_.asDiagonal(); matrix_dirty_ = true; } - void set_q(const VectorXd& q) { - q_ = ruiz_ ? VectorXd(c_s_ * q.cwiseProduct(dx_s_)) : q; + void set_q(const VecRef& q) { + q_ = q; + if (ruiz_) q_ = c_s_ * q_.cwiseProduct(dx_s_); } - void set_A(const MatrixXd& A) { + void set_A(const MatRef& A) { check_eq_A(A); check_eq_b(ruiz_ ? VectorXd(b_.cwiseProduct(inv_de_)) : b_); - A_ = ruiz_ ? MatrixXd(de_s_.asDiagonal() * A * dx_s_.asDiagonal()) : A; + A_ = A; + if (ruiz_) A_ = de_s_.asDiagonal() * A_ * dx_s_.asDiagonal(); compute_AtA(); matrix_dirty_ = true; } - void set_b(const VectorXd& b) { + void set_b(const VecRef& b) { check_eq_b(b); - b_ = ruiz_ ? VectorXd(b.cwiseProduct(de_s_)) : b; + b_ = b; + if (ruiz_) b_ = b_.cwiseProduct(de_s_); } - void set_G(const MatrixXd& G) { - G_ = ruiz_ ? MatrixXd(di_s_.asDiagonal() * G * dx_s_.asDiagonal()) : G; + void set_G(const MatRef& G) { + G_ = G; + if (ruiz_) G_ = di_s_.asDiagonal() * G_ * dx_s_.asDiagonal(); matrix_dirty_ = true; } - void set_h(const VectorXd& h) { - h_ = ruiz_ ? VectorXd(h.cwiseProduct(di_s_)) : h; + void set_h(const VecRef& h) { + h_ = h; + if (ruiz_) h_ = h_.cwiseProduct(di_s_); } - void set_penalty(const VectorXd& penalty) { - penalty_ = ruiz_ ? VectorXd(c_s_ * penalty.cwiseQuotient(di_s_)) : penalty; + void set_penalty(const VecRef& penalty) { + penalty_ = penalty; + if (ruiz_) penalty_ = c_s_ * penalty_.cwiseQuotient(di_s_); } const Solution& solution() const { return sol_; } @@ -495,12 +501,12 @@ class Solver { } // Consistency certificate for Ax = b, refreshed whenever A or b changes. - void check_eq_A(const MatrixXd& A) { + void check_eq_A(const MatRef& A) { eq_infeas_lb_ = 0.0; if (m_ == 0 || !settings.check_eq_consistency) return; eq_cert_.set_A(A); } - void check_eq_b(const VectorXd& b) { + void check_eq_b(const VecRef& b) { eq_infeas_lb_ = 0.0; if (m_ == 0 || !settings.check_eq_consistency) return; eq_infeas_lb_ = eq_cert_.bound(b); diff --git a/include/elastiqp/pdal.hpp b/include/elastiqp/pdal.hpp index 5ad1195..450761b 100644 --- a/include/elastiqp/pdal.hpp +++ b/include/elastiqp/pdal.hpp @@ -79,9 +79,9 @@ class Solver { // Inner-loop row classification by the unclamped dual estimate ztilde. enum class RowState : unsigned char { kInactive, kActive, kSaturated }; - 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_ = q.size(); m_ = b.size(); p_ = h.size(); @@ -161,19 +161,19 @@ class Solver { compute_AtA(); } - 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)); } @@ -183,34 +183,40 @@ class Solver { // Setters store data in the current Ruiz frame; matrix changes invalidate // the factorization. - void set_Q(const MatrixXd& Q) { + void set_Q(const MatRef& Q) { Q_ = 0.5 * (Q + Q.transpose()); if (ruiz_) Q_ = c_s_ * dx_s_.asDiagonal() * Q_ * dx_s_.asDiagonal(); matrix_dirty_ = true; } - void set_q(const VectorXd& q) { - q_ = ruiz_ ? VectorXd(c_s_ * q.cwiseProduct(dx_s_)) : q; + void set_q(const VecRef& q) { + q_ = q; + if (ruiz_) q_ = c_s_ * q_.cwiseProduct(dx_s_); } - void set_A(const MatrixXd& A) { + void set_A(const MatRef& A) { check_eq_A(A); check_eq_b(ruiz_ ? VectorXd(b_.cwiseProduct(inv_de_)) : b_); - A_ = ruiz_ ? MatrixXd(de_s_.asDiagonal() * A * dx_s_.asDiagonal()) : A; + A_ = A; + if (ruiz_) A_ = de_s_.asDiagonal() * A_ * dx_s_.asDiagonal(); compute_AtA(); matrix_dirty_ = true; } - void set_b(const VectorXd& b) { + void set_b(const VecRef& b) { check_eq_b(b); - b_ = ruiz_ ? VectorXd(b.cwiseProduct(de_s_)) : b; + b_ = b; + if (ruiz_) b_ = b_.cwiseProduct(de_s_); } - void set_G(const MatrixXd& G) { - G_ = ruiz_ ? MatrixXd(di_s_.asDiagonal() * G * dx_s_.asDiagonal()) : G; + void set_G(const MatRef& G) { + G_ = G; + if (ruiz_) G_ = di_s_.asDiagonal() * G_ * dx_s_.asDiagonal(); matrix_dirty_ = true; } - void set_h(const VectorXd& h) { - h_ = ruiz_ ? VectorXd(h.cwiseProduct(di_s_)) : h; + void set_h(const VecRef& h) { + h_ = h; + if (ruiz_) h_ = h_.cwiseProduct(di_s_); } - void set_penalty(const VectorXd& penalty) { - penalty_ = ruiz_ ? VectorXd(c_s_ * penalty.cwiseQuotient(di_s_)) : penalty; + void set_penalty(const VecRef& penalty) { + penalty_ = penalty; + if (ruiz_) penalty_ = c_s_ * penalty_.cwiseQuotient(di_s_); } // Recompute the Ruiz scaling from the user frame and carry the warm-start @@ -266,12 +272,17 @@ class Solver { } // Seed the next solve from user-frame (x, y, z); 0 keeps the default rho/mu. - 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, double rho = 0.0, double mu_eq = 0.0, double mu_in = 0.0) { - x_ = ruiz_ ? VectorXd(x.cwiseQuotient(dx_s_)) : x; - if (m_ > 0) y_ = ruiz_ ? VectorXd(c_s_ * y.cwiseQuotient(de_s_)) : y; - z_ = ruiz_ ? VectorXd(c_s_ * z.cwiseQuotient(di_s_)) : z; + x_ = x; + if (ruiz_) x_ = x_.cwiseQuotient(dx_s_); + if (m_ > 0) { + y_ = y; + if (ruiz_) y_ = c_s_ * y_.cwiseQuotient(de_s_); + } + z_ = z; + if (ruiz_) z_ = c_s_ * z_.cwiseQuotient(di_s_); rho_ = rho > 0 ? rho : settings.rho; mu_eq_ = mu_eq > 0 ? mu_eq : settings.mu_eq_init; mu_in_ = mu_in > 0 ? mu_in : settings.mu_in_init; @@ -456,12 +467,12 @@ class Solver { } // Consistency certificate for Ax = b, refreshed whenever A or b changes. - void check_eq_A(const MatrixXd& A) { + void check_eq_A(const MatRef& A) { eq_infeas_lb_ = 0.0; if (m_ == 0 || !settings.check_eq_consistency) return; eq_cert_.set_A(A); } - void check_eq_b(const VectorXd& b) { + void check_eq_b(const VecRef& b) { eq_infeas_lb_ = 0.0; if (m_ == 0 || !settings.check_eq_consistency) return; eq_infeas_lb_ = eq_cert_.bound(b);