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..33cc42b 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,32 @@ while (running) { } ``` +### C + +```c +#include "elastiqp/elastiqp_c.h" + +// One-shot solve; penalty >= ELASTIQP_INF makes a row hard. +ElastiQPProblem qp = {.n = n, .m = m, .p = p, .Q = Q, .q = q, .A = A, .b = b, + .G = G, .h = h, .penalty = penalty}; +ElastiQPSolution res = {.x = x, .z = z}; // caller-owned outputs; NULL ones are skipped +elastiqp_solve_once(&res, &qp, NULL); // NULL: default settings + +// Repeated solves +ElastiQPWorkspace* work; +elastiqp_setup(&work, &qp, NULL); +while (running) { + elastiqp_set_q(work, q_k); elastiqp_set_h(work, h_k); elastiqp_set_b(work, b_k); + elastiqp_set_G(work, G_k); // only rows that differ are refreshed + if (elastiqp_solve(&res, work) != ELASTIQP_SOLVED) { /* ... */ } +} +elastiqp_free(work); +``` + +Pass settings as an `ElastiQPSettings` struct or set them by name: `elastiqp_set_option(work, "eps_abs", 1e-8)`. Use `elastiqp_num_options` and `elastiqp_option_info` to list names, types, defaults, ranges, and descriptions. + +With CMake: `find_package(elastiqp)` and link `elastiqp::elastiqp_c`. + ### Python ```python 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/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 3dba2f9..eeed7f0 100644 --- a/include/elastiqp/das.hpp +++ b/include/elastiqp/das.hpp @@ -76,9 +76,9 @@ class Solver { kDropped }; - void setup(const MatrixXd& Q, const VectorXd& q, const MatrixXd& A, - const VectorXd& b, const MatrixXd& G, const VectorXd& h, - const VectorXd& penalty) { + void setup(const MatRef& Q, const VecRef& q, const MatRef& A, + const VecRef& b, const MatRef& G, const VecRef& h, + const VecRef& penalty) { n_ = static_cast(q.size()); m_ = static_cast(b.size()); p_ = static_cast(h.size()); @@ -158,17 +158,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,25 +176,25 @@ 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_h(const VectorXd& h) { + void set_G(const MatRef& G) { set_rows(G, m_); } + void set_h(const VecRef& h) { rhs_.tail(p_) = h; rhs_dirty_ = true; } - void set_penalty(const VectorXd& penalty) { + void set_penalty(const VecRef& penalty) { penalty_ = penalty; penalty_dirty_ = true; } @@ -203,13 +203,24 @@ class Solver { // Seeds states and multipliers from (x, y, z) on the next solve(), replacing // the internal warm start. - void set_warm_start(const VectorXd& x, const VectorXd& y, const VectorXd& z) { + void set_warm_start(const VecRef& x, const VecRef& y, const VecRef& z) { warm_x_ = x; warm_y_ = y; warm_z_ = z; explicit_warm_ = true; } + // Can relax/differentiate if the last solve succeeded and no data changed + bool can_relax() const { + if (!have_solution_ || Q_dirty_ || penalty_dirty_ || rhs_dirty_) { + return false; + } + for (char c : col_dirty_) { + if (c != 0) return false; + } + return true; + } + RowState row_state(int i) const { return state_[static_cast(m_ + i)]; } @@ -346,13 +357,13 @@ class Solver { // start untouched. The first call after setup() allocates its workspace. const Solution& relax(double kappa, double tol = 1e-6, int max_iter = 50) { if (p_ == 0 || kappa <= 0.0) return sol_; - bool stale = Q_dirty_ || penalty_dirty_ || rhs_dirty_; - for (char c : col_dirty_) stale |= c != 0; - if (stale && have_solution_) { - sol_.status = Status::kUnsolved; - sol_.converged = 0; + if (!can_relax()) { + if (have_solution_) { // data changed since the solve + sol_.status = Status::kUnsolved; + sol_.converged = 0; + } + return sol_; } - if (stale || !have_solution_) return sol_; if (!rw_.ready) relax_alloc(); RelaxWork& w = rw_; const double kappa_s = c_ * kappa; @@ -670,7 +681,7 @@ class Solver { } // Marks only rows that actually changed. - void set_rows(const MatrixXd& R, int offset) { + void set_rows(const MatRef& R, int offset) { for (int i = 0; i < R.rows(); ++i) { bool same = true; for (int k = 0; k < n_ && same; ++k) same = Ct_(k, offset + i) == R(i, k); diff --git a/include/elastiqp/elastiqp_c.h b/include/elastiqp/elastiqp_c.h new file mode 100644 index 0000000..a741346 --- /dev/null +++ b/include/elastiqp/elastiqp_c.h @@ -0,0 +1,233 @@ +/* + * 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 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. + * + * 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 set_X/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 + +/* Solver status codes, equal to elastiqp::Status. */ +#define ELASTIQP_UNSOLVED 0 /* no solve yet */ +#define ELASTIQP_SOLVED 1 +#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 */ + 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, column-major */ + const double* b; /* m */ + const double* G; /* p x n, column-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 status 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 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 primal_obj; /* penalty term included */ + double primal_res; + double dual_res; + double duality_gap; + + double setup_time; /* seconds; elastiqp_solve_once only */ + double solve_time; /* seconds */ +} ElastiQPSolution; + +/* Opaque solver workspace. */ +typedef struct ElastiQPWorkspace ElastiQPWorkspace; + +ELASTIQP_C_API void elastiqp_default_settings(ElastiQPSettings* settings); + +/* 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->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, + * 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->status. */ +ELASTIQP_C_API int elastiqp_solve(ElastiQPSolution* 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_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); + +/* + * 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->status. + */ +ELASTIQP_C_API int elastiqp_relax(ElastiQPSolution* res, + ElastiQPWorkspace* work, double kappa, + double tol, int max_iter); + +#ifdef __cplusplus +} +#endif + +#endif /* ELASTIQP_C_H */ 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); diff --git a/src/elastiqp_c.cc b/src/elastiqp_c.cc new file mode 100644 index 0000000..e356ea4 --- /dev/null +++ b/src/elastiqp_c.cc @@ -0,0 +1,534 @@ +// C interface to the dual active-set backend; see elastiqp/elastiqp_c.h. + +#include "elastiqp/elastiqp_c.h" + +#include +#include +#include +#include +#include +#include + +#include "elastiqp/das.hpp" + +using Eigen::MatrixXd; +using Eigen::VectorXd; +using MatrixMap = Eigen::Map; +using VectorMap = Eigen::Map; +namespace das = elastiqp::das; + +struct ElastiQPWorkspace { + das::Solver solver; + MatrixXd Q; + VectorXd penalty, zero; +}; + +namespace { + +using Clock = std::chrono::steady_clock; + +double Seconds(Clock::time_point since) { + return std::chrono::duration(Clock::now() - since).count(); +} + +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. +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 column-major Q, mirrored; zero for NULL. +void StageQ(MatrixXd& dst, const double* Q, int n) { + if (Q) { + dst = MatrixMap(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(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 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->status = static_cast(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->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->status; +} + +// 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; + } +} + +// 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); + w.penalty.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); + // 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; +} + +} // namespace + +extern "C" { + +void elastiqp_default_settings(ElastiQPSettings* settings) { + if (settings) *settings = ToC(das::Settings{}); +} + +const char* elastiqp_status_name(int status) { + switch (status) { + 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_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); + } + 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 < 0) 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(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; + 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); + }); + if (flag < 0) 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_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); + return 0; +} + +int elastiqp_set_q(ElastiQPWorkspace* work, const double* q) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + work->solver.set_q(OrZero(*work, q, work->solver.n())); + return 0; +} + +int elastiqp_set_A(ElastiQPWorkspace* work, const double* A) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + 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) { + 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->solver.set_b(VectorMap(b, m)); + return 0; +} + +int elastiqp_set_G(ElastiQPWorkspace* work, const double* G) { + if (work == nullptr) return ELASTIQP_ERR_INVALID_ARG; + 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) { + 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->solver.set_h(VectorMap(h, p)); + return 0; +} + +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; + if (penalty == nullptr || !ValidPenalty(penalty, p)) { + return ELASTIQP_ERR_INVALID_ARG; + } + StagePenalty(work->penalty, penalty); + work->solver.set_penalty(work->penalty); + 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; + work->solver.set_warm_start(VectorMap(x, work->solver.n()), + OrZero(*work, y, work->solver.m()), + OrZero(*work, z, work->solver.p())); + return 0; +} + +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) { + return Fail(res, ELASTIQP_ERR_INVALID_ARG); + } + 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(); + const elastiqp::Solution& sol = work->solver.relax(kappa, tol, max_iter); + solve_time = Seconds(t0); + return Extract(res, sol); + }); + if (flag < 0) 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..cb9255a --- /dev/null +++ b/tests/test_c_api.c @@ -0,0 +1,425 @@ +/* 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}; + ElastiQPSolution res = {x, NULL, z, t}; + ElastiQPSettings s = tight(); + 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.primal_obj + 0.75) < 1e-8); + check("one-shot: default settings (NULL)", + elastiqp_solve_once(&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}; + ElastiQPSolution res = {x, NULL, z, t}; + ElastiQPSettings s = tight(); + 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_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)); +} + +/* 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}; + ElastiQPSolution 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("set_h", elastiqp_set_h(work, h) == 0); + elastiqp_solve(&res, work); + check("set_h: constraint released", + near(x, x2, 2, 1e-8) && near(z, z2, 1, 1e-8)); + + { + const double b2[] = {-0.5}; + elastiqp_set_b(work, b2); + elastiqp_solve(&res, work); + check("set_b", res.status == ELASTIQP_SOLVED && near(x, x3, 2, 1e-8)); + } + { + const double q2[] = {-1, -3}, x4[] = {-0.5, 2.0}; + elastiqp_set_q(work, q2); + elastiqp_solve(&res, work); + 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_set_penalty(work, w2); + elastiqp_solve(&res, work); + 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_set_Q(work, Q2); + elastiqp_set_A(work, A2); + elastiqp_set_G(work, G2); + elastiqp_solve(&res, work); + check("set_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.status == ELASTIQP_MAX_ITER); + elastiqp_set_settings(work, &s); + elastiqp_solve(&res, work); + check("set_settings: restored", res.status == ELASTIQP_SOLVED); + } + elastiqp_free(work); +} + +/* + * 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 column_major(void) { + const double Q[] = {1, 0, 0, 0, 1, 0, 0, 0, 1}, q[] = {-1, -1, -1}; + 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]; + ElastiQPProblem qp = {3, 0, 2, Q, q, NULL, NULL, G, h, w}; + ElastiQPSolution res = {x, NULL, NULL, NULL}; + ElastiQPSettings s = tight(); + elastiqp_solve_once(&res, &qp, &s); + 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. 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, -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}; + ElastiQPSolution r1 = {x1}, r2 = {x2}; + /* 2 x1 + x2 = 1, x1 + 2 x2 = 0 */ + const double x_ref[] = {2.0 / 3, -1.0 / 3}; + elastiqp_solve_once(&r1, &qp, NULL); + qp.Q = Qupper; + 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)); +} + +/* 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}; + ElastiQPSolution res = {x, NULL, NULL, t, z_t}; + elastiqp_solve_once(&res, &qp, NULL); + check("ELASTIQP_INF penalty is a hard row", + res.status == ELASTIQP_SOLVED && fabs(x[0] - 1) < 1e-8 && + fabs(t[1]) < 1e-12 && isinf(z_t[1])); +} + +/* + * 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 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}; + 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); + 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: matches a fresh solve", + res.status == ELASTIQP_SOLVED && near(x, x_ref, 2, 1e-8)); + 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); +} + +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}; + ElastiQPSolution 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.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.status == ELASTIQP_SOLVED && 1 - x[0] - x[1] + t[0] > 0 && + z[0] > 0 && z[0] < 10 && z_t[0] > 0); + 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); + 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); +} + +/* 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}; + ElastiQPSolution 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_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); + 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}; + ElastiQPSolution res = {x}; + ElastiQPWorkspace* work = (ElastiQPWorkspace*)&qp; + res.iters = 7; + check("missing penalty rejected, stats cleared", + 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_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + bad[0] = -1; + check("negative penalty rejected", + elastiqp_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + bad[0] = NAN; + check("NaN penalty rejected", + 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("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_solve_once(&res, &qp, NULL) == ELASTIQP_ERR_INVALID_ARG); + check("NULL workspace rejected", + elastiqp_solve(&res, NULL) == 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 && + 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(); + column_major(); + q_upper_triangle(); + inf_convention(); + set_G_partial(); + 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..1136588 --- /dev/null +++ b/tests/test_c_parity.cc @@ -0,0 +1,234 @@ +// 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 { + +// Copies of the problem data, kept alive for the C calls. +struct CData { + MatrixXd 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; + 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(); + 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.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); +} + +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_solve_once(&sol.res, &cp, &cs); + worst = std::max(worst, Diff(sol, ref)); + } + Check("solve_once 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_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); + cd.G = qp.G; + cd.A = qp.A; + ref.set_G(qp.G); + ref.set_A(qp.A); + elastiqp_set_G(work, cd.G.data()); + elastiqp_set_A(work, cd.A.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_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_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; + 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.status == 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; +} 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())