A Dependency-Free Mathematical Optimization Engine
PrimalSolver is a dependency-free mathematical optimization solver written from scratch in C99, supporting LP, MILP, QP, QCQP, SOCP, SDP and mixed-integer conic optimization.
- 48 C ports of MOSEK examples (47 distinct) reproduced with matching optimal values
- 4242 reliability checks
- Full primal/dual solution certificates
- Hybrid engine: presolve + sparse normal-equations IPM for large LPs (90000-var transport in ~0.6 s) and large sparse QPs (n=20000 in ~0.06 s)
- Zero external dependencies (libc + libm only)
PrimalSolver is an independent mathematical optimization engine, implemented from scratch in C99, with extensive compatibility validation against the public MOSEK 11.0 examples. It is not a wrapper, and it is not a teaching implementation of the simplex method: it is a general-purpose solver built around a single internal conic representation, with optimality certificates exposed through a public API.
Most solvers return a solution and ask you to trust it. PrimalSolver returns a certificate:
solver
│
▼
primal solution + dual solution (y, slc/suc, slx/sux)
│
├── stationarity c + Qx + A'y + z = 0
├── primal feasibility
├── dual feasibility
└── complementarity
│
▼
independently verifiable certificate (via public getters)
Every solution comes with duals that satisfy stationarity, primal/dual feasibility and complementarity, and strong duality is reported as |pobj − dobj|. All of it is retrievable through public getters, so the caller can verify the result deterministically instead of trusting the solver.
Infeasibility and unboundedness are certificates as well. PRIMAL_getdualray
returns y (numcon entries) with Σ_i y_i·b_i^act > 0 and A'y ≤ 0 on every
nonnegative variable — a linear combination of the constraints that no
nonnegative point can satisfy, so the model has no feasible solution — and
PRIMAL_getprimalray returns ρ (numvar entries), a recession direction of the
feasible set along which the objective improves without bound (c'ρ < 0 for a
minimization, > 0 for a maximization). Both are scaled to max |·| = 1. A
vector is handed out only after it has been measured in the standard form the
solver actually solved, and the PRIM_INFEAS_CER / DUAL_INFEAS_CER statuses are
reported only together with the matching ray: a status names a certificate, so it
does not outlive one. Where there is no vector to give, the getters answer
PRIMAL_RES_ERR_ARG — the same code a solution getter uses for a problem that has
no solution.
That rule now holds on every route, and it did not before T98. Three places
declared a certificate they never produced — the tangent-cut route in both
directions, and a mixed-integer search whose relaxation is unbounded — and a
fourth declared a primal-infeasibility certificate from a basis the caller had
handed in, which is a statement about the basis. None of them can carry a ray:
the cut route solves in cut space, whose rows are tangents and not the
model's constraints, so a witness of it is not one of it (lifting it is not
sound), and branch-and-bound never lifts the relaxation's direction into the
integer model. All four now publish UNKNOWN and carry the model verdict in
prosta, which answers the question about the model. The basis case was also a
crash: it raised has_sol without allocating a solution — opt_prepare is the
only place t->x is allocated — so the next PRIMAL_getxx copied from a null
pointer. T98 measures both halves, including a positive control: a genuinely
infeasible LP still answers PRIM_INFEAS_CER with a ray that measures.
That fix exposed a second conflation, closed by T99: has_sol was answering
two questions at once — is there a point to hand out, and was a verdict
reached. Because prosta and solsta were read through it, every route that
concluded without a point had to raise it, and the solution getters then
answered OK on the all-zero buffer opt_prepare left behind. On an infeasible
LP that is self-contradictory output: PRIMAL_getxx returns x = 0 while
PRIMAL_getprimalinfeas, reading that same buffer, returns by how much x = 0
violates the model. has_sol now means only the first question, and a verdict
travels in rc plus prosta (solsta too where a certificate has its ray
beside it), so the four routes that had to raise has_sol to be readable — the
LP/QP non-optimal branch, the cut route in both directions, and a MIP stopped by
the node cap with no incumbent — no longer do, and every solution getter
refuses. Certificates are untouched: the ray getters gate on their own
has_dray/has_pray, so an infeasible LP still reports PRIM_INFEAS_CER and
still hands out the ray that measures, with PRIMAL_getxx refusing beside it.
T99 asserts exactly that pairing on eight cases, two of them positive
controls: an optimal LP and a node-capped MIP that does have an incumbent.
A third conflation was between a model and its own representation. The
tangent-cut route caps every bar entry at ±SDP_BIGM (1e6) so that its very
first LP — which has no cuts yet — has a finite answer. A model unbounded
through a bar stopped on that cap and the answer was published as a statement
about the model: rc = OK, solsta = OPTIMAL, pobj = dobj = -2e6, with a dual
bar Z = C − Σ y_i A^i outside the dual cone because the cap's multipliers never
enter that sum — the hole T89 had warned about. Both feasibility getters were
already right (getprimalinfeas = 0: the rows are satisfied;
getdualinfeas = 1: λmin(Z) = −1); the verdict was the part that did not
listen. What decides now is a recession direction measured on the model as
written — its rows, its bounds, and the cone relaxed to a nonnegative
diagonal, the widest relaxation an LP can write for a PSD matrix. The LP only
supplies the candidate, so widening that relaxation can lose a ray and cannot
invent one. A direction that measures means the model has no finite value:
rc = ERR_UNBOUNDED, prosta = DUAL_INFEAS, and no point published — not even
getbarxj/getbarsj — because the direction lives partly inside a bar while
PRIMAL_getprimalray has numvar scalars as its object, the deviation T85
declared for the conic routes. A direction that does not measure means no claim
at all: TRM_MAX_ITER, prosta = UNKNOWN. That verdict lives in one
function, bar_cap_verdict, which both capping routes call after their own loop:
the tangent-cut route (± entry ≤ SDP_BIGM as R_+ rows in its LP) and the
conic build (the same rows, same constant, in optimize_conic). A bar-plus-cone
task never reaches the cut route — numcones > 0 skips it — so on that shape the
conic build was the only publisher, and before this it answered -2e6 for one bar
and -4e6 for two, with getdualinfeas already saying 1. The ray test is now
cone-aware: a candidate direction must lie in every cone of the task, so a
model that escapes through a cone as well as through a bar is still recognized,
and the old numcones == 0 refusal to conclude on that path is gone. T100
asserts the verdict on eleven models: five run on both routes (A..E, because
what is claimed is the verdict and not the path), a scalar unboundedness that
still hands out its ray beside the certificate (F), and five on the conic build —
an unbounded bar, which also proves from the public log line that it was the conic
build asking about the cap (G), the same model starved to one iteration (L), a
bounded bar at 2 (H), a bar at r = 1000 answering 2000 (I), and two
bars with distinct off-diagonals whose published blocks come out [[2,2],[2,2]]
and [[3,3],[3,3]] and swap when the two rows are swapped — the placement control
for the compressed-triangle column blocks, where a shared base would make the
model infeasible rather than merely inaccurate. Two consequences are declared
rather than hidden: a model whose infimum is merely not attained inside a bar
(max −x0 under [[x0,1],[1,x1]] ⪰ 0, which needs x1 → ∞) parks on the cap
with no recession direction and is answered "no verdict" instead of with a value;
and the placement control above has no broken-build behind it — re-introducing the
aliased offsets was refused by this session's permission level, so that one
assertion is a measured invariant, not a measured negative control. On the same
path the PSD tangent rows are now a row kind of their own (CR_PSDCUT): the dual
assembly used to walk them as nonlinear cuts and index cutcol slots that were
never written for them, reading uninitialised memory into the reduced-cost vector.
No published figure moved — the whole 68-sample corpus, every [route]/[cones]
trace line included, is byte-identical before and after.
A cone adds a side the triple cannot see, and a stalled conic run is not an
answer either way. A cone has no rows to leave a residual on: the vertex of a QUAD
block satisfies every equation the model has, so rel_pri, rel_dual and
rel_gap — the three figures every other guard of this family reads — can all
close there on a point that is nowhere near the optimum. What says so is the
reduced cost of the block, which has to lie in the dual cone K*, and it is
judged in PRIMAL_optimize rather than in a route because that is the only place
the model the user wrote is answered (a shadow task calls a route directly, and
its cones belong to the encoder). Measured before the guard: min −x0 on one free
3-D QUAD block with no rows at all returned rc = OK, solsta = OPTIMAL,
pobj ≈ 5e-9 at x = 0, with a cone-dual violation of exactly 1. A conic run that
exhausts its iterations has likewise said something about its own trajectory and
nothing about the model, and TRM_MAX_ITER is the answer that leaves the user with;
two questions are settled there too, each on the relaxation that makes its own
answer sound — infeasibility on a set containing the cones (only linear
consequences of membership: t ≥ 0, t ≥ ±u_i for a QUAD, t ≥ 0, u ≥ 0 for the
rest), because an LP infeasible on a subset of the model proves nothing;
unboundedness on a set inside them (a polyhedral inner approximation), because a
direction of a wider set is not a direction of the model. The two relaxations are
opposite and neither substitutes for the other. Neither verdict carries a vector:
the conic routes publish no Farkas ray (T85), so prosta says PRIM_INFEAS/
DUAL_INFEAS, solsta says UNKNOWN, and nothing about the point is published —
the family T98 established. Reading the block back also needed the sense algebra
written down: the route's min form has cost s·c and slack s·c − A'π, and y is
published as −s·π, so the user-space dual is d = s·(c + A'y). The reading
c + s·A'y that was there first is identical on a MINIMIZE model and differs by
exactly 2c on a MAXIMIZE one — a sign flip of the cone, not a rounding — and a
MAX model whose cone members carry no cost (where this was first pinned) cannot
discriminate them. T101 G is the model that can: max −(x0+x1+x2) with
x0 = 5, x1 = 4 on a QUAD block, whose correct dual is the boundary point
(5/3, −4/3, 1) and whose wrong reading is (−1/3, −10/3, −1). Reverting just
that factor does not move a digit there — it flips the whole verdict to
rc = 1007 with every getter refusing, and T101 is green on dobj = pobj = −6.
Fixing the ray test on that shape exposed a third reading error, of the NaN > viol family of T90: the u = 0 face of a PEXP was taken as "no violation",
where the cone in fact closes onto {u = 0, v = 0, t ≥ 0} — so a candidate
direction escaping along v passed the membership test and a bounded model was
answered unbounded. T101 asserts nine models (A..I) with the dual block, its
boundary identity t = |(u,v)| and its complementarity d·x = 0 recomputed in the
test from the published y, and keeps two declared deviations visible: a cone
member with a finite bound leaves its block unmeasured rather than measured
wrong (H), and the bar route publishes slx/sux = 0 on the cone's own variables
where the true block dual is (5/3, −4/3, 1) — getdualinfeas therefore reads that
dual from c and y directly, not from the published slacks, so the zeros do not
become a violation (I).
Two limits are declared rather than hidden. A dual ray whose support needs the
cap row of a ranged variable (x_j ≤ ux_j), or both sides of a ranged
constraint, has no one-entry-per-constraint image in the user's model, so it is
not published even though the reduced-space argument behind it is sound. And
explicit witnesses are written only by the tableau simplex — an interior-point
route has none, so its last iterate is offered to the same measurement and the
measurement decides. The conic routes now publish both rays of a bar-free
model. The primal ray (a recession direction) comes from
model_lp_witness's inner relaxation and is checked by sdp_ray_measures —
ρ ∈ K for every cone, Aρ = 0 on every row, the bound signs, c'ρ < 0,
normalized to max|·| = 1 — with solsta = DUAL_INFEAS_CER and getprimalray
answering (T101 A/B). The dual ray of an infeasibility comes from the
outer relaxation's Farkas witness and is checked against the model's own
cones — d = A'y ∈ K* on every block, the row multipliers signed, the RHS
combination negative — with solsta = PRIM_INFEAS_CER and getdualray
answering (T101 C/J, the §3.7 shape). A candidate that does not measure
proves nothing and the verdict stays in prosta. With bars the compressed block
still has no image as numvar scalars, so no ray is published there. The same
pass also refuses a conic point that lies outside the cones: the route's
three figures do not see membership, and min x0 with x0 = -1 over QUAD
answered OPTIMAL at x = (-1,0,0), outside by 0.5, before conic_primal_cone_worst
put the fourth figure into the verdict. MIP rays remain open, declared.
Feasibility of the cones is a fourth, separate number. rel_pri counts the
equality rows only, so none of the three published figures says by how much the
point is inside K; membership used to be guaranteed by the backtrack test and
never measured. GMB_DBG=1 now prints one line per solve that published into a
conic, SDP or putqconk model:
[cones] task rel_slack=-2.88e-09 pri_slack=-5.77e-09 where=cone 0 PPOW a=0.666666666667{...}
rel_slack is the worst signed slack over bounds, cones, quadratic rows and
PSD blocks, divided by 1 + (size of that object) — max |·| of the members for
a cone, the spectral norm for a bar, the bound itself for a bound — because the
tolerance the task declares is a relative one. Reading a block in its own
homogeneous units matters: the t-space form of PPOW(a) shrinks a geometric
violation by 1/(a·t^(a-1)), which is 6.4 on regression_regularized and 3333 at
a=1/3, t=1e-6, and outside {t >= 0, u >= 0} pow() returns NaN, which
compares false against any violation and so reads as "no violation". Measured on
the -O2 build: 100 published solutions in the suite (94 on HEAD) and 31
across the samples, worst 4.5e-09. Five of the hundred sit outside the 1e-8
the task declares — {4.8e-05, 2.7e-04, 1.7e-03, 2.6e-02, 0.67}, a quadratic row,
a 2×2 bar, two bounds and one QUAD cone — and all five print identically on
HEAD, so the six publications the new guard added are all inside. The
tangent-cut route now accepts on
that same quantity (T90): it used to stop as soon as the t-space residual
fitted tol·(1+|x0|), which published t = -4.7e-11 — a point where the cone has
no value at all — with PRIMAL_RES_OK.
Rather than shipping independent solvers per problem family, PrimalSolver maps every supported class onto one shared internal representation:
LP → QP → SOCP → exponential/power cones → SDP → mixed-integer conic
QCQP constraints are encoded exactly through an eigendecomposition into rotated quadratic cones; semidefinite programs are solved by a primal-dual interior point on the PSD blocks (with a λ_min tangent-cut fallback); SOC and SDP are handled in the same conic loop; mixed-integer conic problems are solved by branch & bound over conic relaxations. This keeps the solver surface small and the behaviour consistent across families.
- Optimality certificates: primal + duals satisfying stationarity, feasibility and complementarity, fully checkable through public getters
- Infeasibility certificates: Farkas rays as vectors (getdualray/getprimalray), normalized tomax |·| = 1and published only when they measure in the form that was solved
- Basis solve: exact evaluation of a given basis (solvebasis), MPS-BAS read/write
- Sensitivity analysis: cost ranges for which the current solution stays optimal (costsensitivity/rhssensitivity)
- The model reads back as it was written: PRIMAL_getarow/PRIMAL_getacoland their*sliceforms hand back the entries of a row or column as subscripts and values, answer how many were written, and refuse without touching the caller's buffers when the space does not fit;getnumanzs/getmaxnumanzs/getnumqobjnzcount what is stored. The scalar read of one entry answers something else — the operator, i.e. the sum of the entries stored at that(i,j), which is the coefficient every route solved — and that split is the subject of its own deviation note below (T107). Variables and constraints carry names —putvarname/putconname,getvarnameidx/getconnameidx, and a lookup by name (getvarname/getidxvar,getconname/getidxcon) — held in two independent tables, so a constraint and a variable may share a name the way an MPS file does (T102). Each table also reads back whole —getallvarname/getallconnamewrite exactlynumvar(numcon) borrowed strings, and an unnamed slot reads""there and fromgetvarnameidxalike, because two readings of one table that disagree are two rules (T103). The two quadratic forms read back whole the same way:PRIMAL_getqobj/PRIMAL_getqconkenumerate the table their own counter counts and therefore answer in different shapes, because the two stores are different — the objective is the user's triplet list in write order, stored zeros included (a written0.0is a write, not a hole) with a cross term stated once, while a constraint'sQis the upper triangle of the dense store, ascending, so a row whose entries cancel reads empty rather than refused (T104). Bar variables are the third name table —putbarname/getbarnameidx/getbarname/getidxbarvar/getallbarname— and they add no rule, only a table: the five rules live in the same two functions and the capacity is the samebarcapas the blocks themselves. The one thing the two scalar tables could not decide is decided here: a bar's name lives in a namespace of its own, because a bar is not one of thenumvarscalars — it has its own index space, its own dimension and its own cone block. So one string may name a variable, a constraint and a bar at the same time, while two bars may not share it (T105). Cone blocks are the fourth table —putconename/getconenameidx/getconename/getidxcone/getallconename— sharing the one capacity ofcone_type/cone_nmem/cone_mem/cone_paramrather than carrying a fifth to keep in step, and its own namespace too. Becausegetconname(a constraint) andgetconename(a cone) differ by one letter, the independence is not asserted in prose but measured: one string naming a variable, a constraint, a bar and a cone at once answers four lookups with four different indices, and a duplicate cone name is refused with the constraint's name still alive (T106C). The objective is the same mechanism over a table of one —putobjname/getobjname— which is why renaming it replaces instead of refusing: withn = 1the uniqueness loop has no second slot to collide with, and that difference from a table is asserted rather than inferred (T106G). The last two getters that still broke that contract were the bar-term list readers:PRIMAL_getbaraidxij/PRIMAL_getbarcidxjstopped at the caller'smaxnumand answeredPRIMAL_RES_OKwith a truncated list, which is a reading that lies — an⟨A,X⟩built from a prefix is a different model. They now count first and refuse without writing, withNULL/NULLas the door that asks for the number alone (T108)
- Solution quality: max primal/dual violation of the current solution, plus
the interior-point relative criteria in MOSEK's form — pfeas/(1+‖b‖∞),dfeas/(1+‖c‖∞),|pobj−dobj|/(1+|pobj|+|dobj|)— declared throughPRIMAL_DPAR_INTPNT_TOL_{PFEAS,DFEAS,REL_GAP}(1e-8 each by default), which is what the conic route is accepted or rejected on, andPRIMAL_DPAR_INTPNT_TOL_NEAR_REL(1000), the reference's near-optimal acceptance factor, which scales all three before the verdict where this IPM is the only conic route
- A "not solved" answer publishes no solution: solstastaysUNKNOWNand every solution getter answersPRIMAL_RES_ERR_ARG, rather than the zeroed vectors the pre-solve left behind. That covers the answers which do reach a verdict and still have no point: an infeasible model, an unbounded one, a node-capped integer search that never found an incumbent (T99), and an answer that stopped on the solver's own cap on the bar entries (T100)
- Feasibility repair: elastic LP — minimum-total-violation point for infeasible problems
- Progress callback: per-iteration/per-round notifications observable externally
- Warm start: user-provided starting point for the IPM
- Reported problem and solution status: PRIMAL_getprostaandPRIMAL_getsolstapublish the reference's two numbers, with the pairing its tables fix, soprostais derived fromsolstarather than stored beside it
- Declared parameters: the 27 PRIMAL_IPAR_*/PRIMAL_DPAR_*ids live in one table (PRIMAL_PARAMSinprimal.c) holding kind, target field, default and inclusive range. Defaults, the range checks of both setters, the getters and the copy into shadow tasks all read that table — there is no second list of parameters to keep in sync — andPRIMAL_getparaminforeports a row's default and bounds. Integers and doubles are separate id namespaces sharing numbers, soPRIMAL_getparaminfotakes the kind as an input. Wired:IPAR_SCALING,IPAR_MIP_MAX_NODES,DPAR_MIP_TOL_{ABS_GAP,REL_GAP,INTHER,FEAS}(the last four used to be compile-time constants;FEASdid not exist at all — it is the tolerance an incumbent is re-verified against, see below)
- LP presolve: empty/singleton-row and empty-column reductions with exact primal+dual postsolve (default on, PRIMAL_IPAR_PRESOLVE)
- I/O: free-format MPS (RHS/RANGES/BOUNDS/MARKER/QSECTION), CPLEX LP, CBF v4 (cones + SDP + CHANGE section).
The MPS and LP readers populate the name tables they read through: what a file
labels is what the task calls it — putvarname/putconname/putobjnameat the index the label names (T112). A file that gives one name to two rows is refusedPRIMAL_RES_ERR_FILEbefore anything reaches the task, because every later section looks rows up by name and the second row would receive no entries at all. CBF keeps its group labels out of those tables: a group namesszvariables at once, and a name table admits one index per name.
- Bound slices and solution slices (T114,T115):PRIMAL_getvarboundslice/getconboundsliceand theput*counterparts,[first, last), reference-shaped; a refused read writes nothing and a refused write validates the whole slice first.PRIMAL_getxxslice/getyslice/getslcslice/getsucslice/getslxslice/getsuxslice/getskxslice/getskcslice,getxc/getxcslice(the constraint activity, read by the samerow_activityasgetpviolcon) andgetreducedcosts((s_l^x)_j − (s_u^x)_j = −(slx+sux)) are a second reading of the solution vectors, asserted against the whole vector entry by entry. The counters are named as the reference names them (getnumanz/getmaxnumanz).T117adds the block-triplet form of A-bar and C-bar (getbarablocktriplet/getbarcblocktripletand their counts): one row per stored lower-triangle entry of every block.T118adds the objective vector (getc/getcslice,putclist/putcslice) and the scalar quadratic entryputqobjij(q_ij = q_ji, lower triangle, replaces the pair and removes it at 0).T119adds the A store in triplet form (getatrip), the row/column nonzero counters and their slices, andputaij(the scalar entry that replaces and collapses the pair).T120adds the list forms (putvarboundlist/putconboundlist,putvartypelist/getvartypelist) and the constant-bound slices (putvarboundsliceconst/putconboundsliceconst), all of which validate the whole input before applying.T121adds the task name (puttaskname/gettaskname/gettasknamelen) and the name lengths (getvarnamelen/getconnamelen/getobjnamelen/getmaxnamelen; a length never counts the terminator, andgettasknamerefuses when the buffer cannot hold it).T122addsputqcon(replace all quadratic terms of all constraints from one triplet list, lower triangle, duplicates accumulate, empty list clears) and the utilitiesgetversion/isinfinity/getresponseclass.T123fixes a branch-and-bound soundness bug: a semi-continuous variable's root box is now relaxed to[0,u], so the LP relaxation keeps the deactivationx=0and its optimum is again a valid lower bound (min xover{0}∪[2,5]is 0, not 2). The MPS writer's round-trip was fixed too: for a file-read model its row 0 is theNobjective row, so the writer emitsNfor it and rows1..numcon-1as constraints, preservingnumcon(T112G). A hand-built model whose row 0 is a real constraint keeps the old path (one extra row on re-read), and names stay positional -- both declared.T124adds the general callback (putcallbackfunc/getcallbackfunc: BEGIN/END of OPTIMIZER/READ/WRITE with theMSKcallbackcodeenumbers) and the response callback (putresponsefunc, invoked when a solve ends with a non-OK code); the detail vectors are NULL, declared.T255completes the phase callbacks: the route that answers also emits BEGIN/END of SIMPLEX or INTPNT (LP/QP), MIO (branch & bound) and CONIC (conic/SDP), so the callback reports which algorithm ran.T256adds the OPF format (reference 16.3) -- reader and writer for the linear/quadratic/conic/ integer core, dispatched by the.opfextension, with the data callbacksREAD_OPF/READ_OPF_SECTION/WRITE_OPFemitted.T257adds the "middle" solution-update codes (CONIC,PRIMAL_SIMPLEX,INTPNT) at each cut round / basis / relaxation.T258closes the per-iteration gap: the engine loops (simplex.c93/36,ipm.c90,socp.c/sdp.c34) emit their code at every internal iteration through a hook (primal_cb_iter, guarded byprimal_cb_iter_onso the common case with no callback is a plain load+branch).T260adds the sub-step codes (IM_LUat every LU factorisation,IM_ORDERat the AMD ordering). Declared deviation: withIPAR_NUM_THREADS > 1the concurrent LP path interleaves the events.T259reads and writes OPF quadratic constraints (the writer emits the Q part, so they round-trip): a convexq(x) <= buses the UP convention andq(x) >= bis rewritten as-q(x) <= -b; an equality or ranged quadratic is refused (declared). OPF declared deviations: thezerocone,[solutions]and[vendor]are not read; the OPFdpowmaps to our RPOW.T125addsgetsolution(the reference's one-call reader: sta/solsta, basis keys,xc,xx,y,slc/suc/slx/sux,snx) plusgetskn/getsnx;sknis SK_UNDEF for cones andsnxis 0, both declared.T126addsputaijlist,chgvarbound/chgconbound(change one side and recompute the key), the matrix store in bulk (getsparsesymmat/appendsparsesymmatlist) andgetprobtype(LO/QO/QCQO/CONIC/MIXED).T127adds the parameter resets (resetintparam/resetdouparam/resetparameters, from the declarative table's defaults) and the 64-bit counter variants (getnumanz64/getmaxnumanz64/getnumqobjnz64/getnumqconknz64).T128adds the contiguous-member cone appends (appendconeseq/appendconesseq) and the bulk row/column setters (putarowlist/putacollist/putarowslice/putacolslice, CSR data validated before writing).T129adds the enum symbolic names (prostatostr/solstatostr/bktostr/conetypetostr/probtypetostr/sktostr/rescodetostr); the text is this solver's own, a declared deviation (the reference's exact strings were not read).T130addscheckversion/getbuildinfo/getcodedesc/getlasterror(the last response code and message follow the last solve) andget/putatruncatetol(stored, not applied: this solver does not truncate A).T131adds the removalsremovecones(compacts the cone list) andremovebarvars(drops bar variables and remaps the A-bar/C-bar terms).T132adds file streams (linkfiletotaskstream/linkfiletoenvstream), function streams on the env (linkfunctoenvstream, inherited by tasks created afterwards), the unlinks and the echo helpers (echotask/echoenv/echointro).T133adds the bar sparsity and per-block index info (getbarasparsity/getbaraidxinfo/getbaraidxand the C variants;idxis this solver's vectorisationi*numbarvar+jfor A,jfor C, a declared deviation).T134addsinputdata/inputdata64(load the linear part in one call; the task must be empty) andgetmemusagetask(a byte estimate).removecons(T134) removes constraints and compacts A/qcon/barA, remapping the remaining rows.removevars(T135) removes variables and reshapes the dense qcon blocks at the new stride, remapping qobj and the cone member lists.T136addsstrtoconetype/strtosk(the inverse of the symbolic names) andgetnumparam(parameters per type in the declarative table).T137starts the AFE layer:appendafes/getnumafeand the sparse-per-row storage ofFandg(putafefentry/putafefrow/putafegand the getters).T138adds the conic domains (append*domainandgetnumdomain/getdomaintype/getdomainn, with MSKdomaintypee's values) that an affine conic constraint is a member of.T139closes the layer with the reference-styleappendacc(domidx, numafeidx, afeidxlist, b)(AFE + domain, adapter over the internal encoder), its gettersgetnumacc/getaccn/getaccdomain/getaccafeidxlist/getaccb, and the linear domainsR/RZERO/RPLUS/RMINUSencoded as rows. The reference convention thatbis subtracted (F x + g - b, as documented for bothappendaccandputdjc) is applied here; the previous+breading was a sign error, andT139passesb = +2for the same(x0+x1-2) ∈ R+model.T140adds the reference-style DJC layer:appenddjcs(empty slots) +putdjc(djcidx, numdomidx, domidxlist, numafeidx, afeidxlist, b, numterms, termsizelist)(an OR of clauses, each a conjunction of domains over affine expressions),putdjcslice, the gettersgetnumdjc/getdjcnumdomain/getdjcnumafe/getdjcnumterm/getdjcdomainidxlist/getdjcafeidxlist/getdjcb/getdjctermsizelistand the totals, plusputdjcname/getdjcname/getdjcnamelen. The big-M encoder is now the internaldjc_encode, and only linear domains are representable in a DJC (a conic domain isERR_ARG, a declared deviation).T141ports the transaction-cost LP of the StackOverflow document (free variables, absolute-value epigraph written as two rows, a cost-free degenerate column the presolve removes) and asserts the face optimum — value0.8plusx_i >= x0_iande'x = 0, not the posted vector.T142covers the two parity hooks of arXiv:2508.03704 / Seven Sins:||w||_1 <= 2on sign-freew, and the non-integer powers|x|^{3/2}/|x|^{4/3}asPPOW(2/3)/PPOW(3/4)(§3.6).T143solves the max-Sharpe ratio of case A3 through the Charnes-Cooper transform (y = kappa*w,y'Sigma y <= 1,e'y = kappa), the ratio itself being non-convex.T144adds the bulk AFE surface (emptyafefrow/emptyafefcol,putafeglist/putafegslice/getafegslice,putafefentrylist,getafeftrip,getafefnumnz).T145adds the bulk ACC surface (appendaccs/getaccs/getaccntot/putaccb) and the sixth name table (putaccname/getaccname/getaccnamelen).getdimbarvarj/getlenbarvarj(lower triangled(d+1)/2) and thegetmaxnum*counters (T116; for var/con they are the current count — the arrays grow in place, a declared deviation — while cone/bar capacities are real).
- Information items (T160): the reference's information database is readable
end to end. PRIMALinftypee(DOU/INT/LINT) and the three contiguous enumsPRIMALdinfiteme(0..115),PRIMALiinfiteme(0..136),PRIMALliinfiteme(0..21) carry the reference's own indices (read fromconstants.html, MOSEK 11.2.4, 2026-09-19).PRIMAL_getdouinf/getintinf/getlintinfread by index,PRIMAL_getnadouinf/getnaintinfby name, andPRIMAL_getinfmax(max index + 1 =END),PRIMAL_getinfname(index -> name) andPRIMAL_getinfindex(name -> index) enumerate the database without knowing the symbols. The items this solver actually measures (objective values, primal/dual violations, solution norms, model sizes and bound-type breakdowns, problem/solution status,OPTIMIZER_TIME, matrix density,RD_*counts) are computed; the items of the phases it does not have (MIO cut separators, folding, remote, per-phase times) answer 0 — a declared deviation, not an invented number. The symbolic constants are also covered:PRIMAL_getsymbcondim/PRIMAL_getsymbcon(1438 entries, max name length 60),PRIMAL_symnamtovalue(name -> decimal value) andPRIMAL_iparvaltosymnam(parameter value -> symbolic name). That table was not copied from a doc page (the live docs do not expose it): it was extracted by calling the reference library's ownMSK_getsymbcondim/MSK_getsymbcon, so it is the reference's own data.
- Reference parity surface (category C, T159): the remaining implementable
entry points of the reference API are now present — getlintparam/putlintparam(64-bit parameter read/write,putparamby name,writeparamfile/readparamfile,printparam), the name generators (generate{varnames,connames,conenames,barvarnames,accnames,djcnames}),getconeinfo/readsummary,optimizetrm/optimizebatch(solve a batch of tasks over one env),primalrepair/dualsensitivity/primalsensitivity/toconic, the ACC/cone/bulk writersputacc/putacclist/putafefrowlist/putcone/putbararowlist, the "new" solution API (solutiondef/getsolutionnew/putsolutionnew/putsolution/getsolutioninfonewplus the per-indexputconsolutioni/putsolutionyi/putvarsolutionj), string/handle I/O (readlpstring;readopfstringreads OPF,readptfstringremains a declared deviation returningERR_ARG;readdatahandle/writedatahandle), basis solve (initbasissolve/solvewithbasis/basiscond),computesparsecholesky(sparse Cholesky viaspchol_factor_ord, real AMDperm;T264), the UTF-8 <-> wide-char conversions (utf8towchar/wchartoutf8,T265), andclonetask/getdualproblem/getinfeasiblesubproblem. Public surface 514 functions. The numerics of the new entry points are cross-checked against hand values inT159(Cholesky of a dense SPD,B x = bon the identity basis, the LP dual ofmin e'x, e'x = 2, a cloned task solving to the same optimum).
- KKT robustness: adaptive regularization for degenerate systems
make # build (gcc -std=c99 -Wall -Wextra -pedantic -O2, 0 warnings)
make test # reliability suite: 4706 checks (out/run_tests)
make run-samples # 171 runnable examples (out/samples/*)
make clean # remove the whole out/ directoryAll build artifacts go into out/ (gitignored); nothing is written into the
repository root or into samples/. No dependencies beyond libm. Clean under
AddressSanitizer + UndefinedBehaviorSanitizer across the whole suite and all
samples.
For the min-normalized problem (max → multiplied by s=−1):
- stationarity: c + Qx + A'y + z = 0, withz = slx + sux
- y_i > 0⇒ row i at its upper bound,- y_i < 0⇒ at its lower bound
- z_j > 0⇒ variable j at its upper bound,- z_j < 0⇒ at its lower bound
- Y/- Zabove are the min-normalized multipliers —- s·the reported- t->yand- t->slx + t->sux(the getters split them into- slc/sucand- slx/suxby sign) — and- b_i^act/- xb_j^actis whichever of- blc/buc(or- blx/bux) the multiplier's own sign selects at the point published.
- dual objective: dobj = cfix − Σ Y_i b_i^act − Σ Z_j xb_j^act − ½ x'Qx, published in the problem's own sense. The constant is a term of the objective as written, so it is added outside thesthat normalizes the sense — multiplied bysit misses by2·cfixon a max problem (measured:T92case A).
- strong duality |pobj − dobj| ~ 0verified byT17andT92(five routes: LP in max sense, QP, dense conic, bar/SDP, power cone). BeforeT92the two sums buildingdobjomittedcfixwhilepobjstarted from it: on the LP ofT17(cfix = 7.5) the published dual was2.0against a primal9.5— the dual of the model minus its constant — and the dense conic report had the same shape.
MOSEK's documentation is readable from this machine (11.2.4: capi/constants.html
for the enumerations, capi/parameters.html for the 277 identifiers with their
defaults and accepted ranges, capi/response-codes.html for the result codes,
capi/accessing-solution.html for which status each optimizer reports), and the
comparisons below were read from it rather than inferred. Where a number is
shared it is pinned by a test.
Solution status (PRIMAL_getsolsta) carries MOSEK's numbering, so a raw status
read against MSKsolsta means the same thing here:
Only the five statuses this solver produces are declared. Before T93 the
numbers were this solver's own and three of them read as a different verdict:
INTEGER_OPTIMAL was 2, which MSKsolsta calls PRIM_FEAS; PRIM_INFEAS_CER
was 4 (PRIM_AND_DUAL_FEAS); DUAL_INFEAS_CER was 5 (PRIM_INFEAS_CER).
Problem status is reported too, by PRIMAL_getprosta, with MOSEK's
MSKprostae numbering:
The reference fixes which pair of numbers may appear together (Tables 7.2 and
7.3 of accessing solution), so prosta is derived from solsta by one
function rather than stored beside it — the two cannot drift out of sync the way
the objective constant did:
The last-but-one row is not hypothetical: when the node cap stops the search with
work still on the stack, the incumbent is a feasible point, not a proof, and
sched does exactly that (100000 nodes, 27 open, incumbent 97.328509 — the true
makespan, never established). Before T94 that case printed
INTEGER_OPTIMAL with rc = OK. The pair is asserted on the raw numbers in
T94 H; the rows where the two halves disagree are the ones T98 added —
an infeasible integer problem used to be the only one, and now every verdict
that has no Farkas vector beside it is in the same situation. t->prosta is an
override field precisely for those: prosta is still derived by one function, so
a solve either follows the table or says, in one place, that the certificate is
missing.
The fourth column is T99. It is read from the getters and not from the flag:
six of those rows carry a verdict and have no answer to hand out, and two of
them — the certificated ones — hand out a ray while refusing the point. T96
had already made a not solved answer refuse; what T99 closes is the routes
that did reach a verdict, which until then raised the flag to be readable.
The same lesson applied to the input side. A bound key, a cone type or a solution
item passed in by its reference number used to select a different object here,
with no error to read: PRIMAL_BK_FX was 4, which is FR in
MSKboundkeye, so a caller fixing a variable asked this solver to free it. The
enumerations now carry the reference's numbers (T94):
Parameter defaults are the reference's where the reference states one and this
solver's behaviour matches (T84 reads each default back from a fresh task and
checks the declarative table against itself; the values below are the ones the
reference publishes):
PRIMAL_IPAR_OPTIMIZER used to default to 0, which after the numbering was
adopted is the reference's CONIC; the value is behaviourally the free choice
here, but a fresh task now reports what the reference reports.
PRIMAL_DPAR_MIP_TOL_REL_GAP used to default to 0 (off) as a declared
deviation, on the argument that pruning on the relative gap alone accepts a
deliberately suboptimal incumbent. With the reference's 1e-4 the whole corpus is
unchanged: the suite (1503 checks at the time of that measurement, now 2607) ran
green at -O0/-O1/-O2/-O3 and so did the 68 samples. The
deviation was therefore dropped and the default aligned.
Deviations that remain, each stated because the reference says otherwise:
- There is no third verdict, and there did not need to be one. MSKsolstahas noNEAR_member; MOSEK's categories are optimal, certificate, and unknown, and a stalled or numerically uncertain iterate is reportedMSK_SOL_STA_UNKNOWNwithMSK_PRO_STA_UNKNOWN. What MOSEK has instead isMSK_DPAR_INTPNT_CO_TOL_NEAR_REL(default 1000, range [1,+inf]): if the prescribed accuracy is not reached, the point is checked against the termination criteria with all tolerances multiplied by that factor, and if it passes it is declared optimal. This solver now has the same rule, asPRIMAL_DPAR_INTPNT_TOL_NEAR_REL(setting it to 1 disables it), and applies it where a measured triple decides: the unified conic IPM's gate. The policy is two-sided. Where an alternative algorithm exists — exp/power, the tangent cut outer approximation — the declared tolerances keep choosing the route, because a factor may not make the solver hand back a point the alternative beats. Where this IPM is the only conic route — bar and SOC models, no exp/power block — the factor is the verdict, somaxcut_sdp(rel_gap1.32e-8) and theT88model (rel_pri1.01e-8) are now publishedOPTIMALexactly as the reference would, instead of being reported "not solved". Nothing else had to be re-tuned: the iterate a path integrates does not depend on the scale of the three tolerances, because the snapshot selector ranks byviol = max(rel_i/tol_i)and its improvement testviol < 0.99·viol_bestis homogeneous of degree 0 — moving a tolerance moves the judgement, never the point (measured on the whole corpus, and locked by T96).
- Two integer tolerances, as the reference has them. Declaring an integer
constraint satisfied and declaring a point feasible are different questions,
and conflating them publishes an infeasible solution: with
PRIMAL_DPAR_MIP_TOL_INTHER(the reference'sMSK_DPAR_MIO_TOL_ABS_RELAX_INT, 1e-5) adopted and nothing else changed,samples/djc1answered the max case withx0 = 10— outside both disjuncts (x0 <= 2OR6 <= x0 <= 7) — atrc = OKandsolsta = INTEGER_OPTIMAL. A big-M disjunction amplifies the rounding: a selector at1 - 5e-6is integral within 1e-5, and rounding it moves its row by 5, not by 5e-6. The reference keeps a second tolerance for the second question,MSK_DPAR_MIO_TOL_FEAS(1e-6, accepted [1e-9;1e-3]), and so does this solver now:PRIMAL_DPAR_MIP_TOL_FEASis whatmip_point_measuresjudges a candidate incumbent against — bounds, rows, cones, semi domains and SOS sets — whileMIP_TOL_INTHERgoverns only the search. A candidate that does not measure is not discarded: the integers are pinned at their rounded values and the node's relaxation is re-solved (fix-and-optimize), which is still a relaxation of that node and loses no integer point. Measured:djc1answersx0 = 0(min) andx0 = 7(max) withgetprimalinfeas = 0at every threshold tried — default, 1e-9, 1e-5, 1e-3 (T95) — so the tolerance that used to be a deviation is the reference's own.
- Mixed-integer parameters live in MOSEK's MIOgroup, and there is noMSK_IPAR_MIP_MAX_NODES: the node budget isMSK_IPAR_MIO_MAX_NUM_BRANCHES(default -1, MOSEK chooses).PRIMAL_IPAR_MIP_MAX_NODESkeeps this solver's name and a concrete default of 100000.
- Three interior-point tolerance sets, as the reference has them. MOSEK splits
INTPNT_TOL_*(LP),INTPNT_CO_TOL_*(conic) andINTPNT_QO_TOL_*(quadratic); this solver now has all three, each read by its own route: the plain set by the LP route,CO_TOL_*by the conic route (sdp_ipm/socp_solve) andQO_TOL_*by the quadratic route (ipm_solve_qp_cscand the dense QPipm_solve_std). Defaults and accepted ranges are the reference's (parameters.html11.2.4: 1e-8,[0,1]). The previous single-set behaviour was a declared deviation, now closed.
- Accepted ranges now match the reference's. Read again from
parameters.html(11.2.4) and aligned 2026-09-22:TOL_PFEAS/TOL_DFEASare[0,1](the earlier table widened them to[DBL_MIN, DBL_MAX], a misread — onlyTOL_REL_GAPis[1e-14,+inf]),MIO_TOL_REL_GAPis[0,+inf](was[0,1]) andMIO_TOL_ABS_RELAX_INTis[1e-9,+inf](was[DBL_MIN,1]). The iteration limits used to require>= 1where MOSEK accepts 0 (its "no limit"): closed — the range is[0, INT_MAX], the stored value is 0 (what a getter reads back) and every use site maps 0 to the largest int (iter_cap), locked by T84 B2.T94F reads theTOL_PFEASrow.
- The optimizer time cap is enforced. PRIMAL_DPAR_OPTIMIZER_MAX_TIME(default -1 = no limit, range[-inf,+inf]) sets a wall-clock deadline thatPRIMAL_optimizepublishes to the IPM loops (ipm_set_deadline) and to the B&B loop; a fired cap returnsPRIMAL_RES_TRM_MAX_TIME(1008), and the conic verdicts are skipped (a timed-out solve has no answer to judge). MOSEK names the codeMSK_RES_TRM_MAX_TIME; the numbering is our own (as for everyrc).PRIMAL_DPAR_MIO_MAX_TIME(same default/range) caps the mixed-integer phase only: the B&B obeys the tighter of the two deadlines. Locked by T174/ T175.
- The lower objective cut is enforced. PRIMAL_DPAR_LOWER_OBJ_CUT(default-inf= no cut, range[-inf,+inf]) is checked in the IPM loops: a primal-feasible point whose objective is below the cut proves the optimum is below it, and the solve terminates withPRIMAL_RES_TRM_OBJECTIVE_RANGE(1009).PRIMAL_DPAR_UPPER_OBJ_CUT(default+inf) is its dual-side twin: a dual-feasible point whose dual objective (b'yfor an LP) is above the cut proves the optimum is above it. The upper cut is enforced on the LP route only (the QP/conic routes do not compute that bound), and both cuts are wired for minimization — a maximization would need the mirrored pair. Locked by T176/T177.
- The PSD tolerance is declared. PRIMAL_DPAR_SEMIDEFINITE_TOL_APPROX(default 1e-10, range[1e-15,+inf]) is read as the relative factor of the encoder's convexity threshold (bad_thr = tol*max(lmax,1)) and of the witness PSD-cone check (e >= -tol*(1+emax)). The reference's is an absolute tolerance — a declared deviation. Locked by T178.
- getprimalinfeasdoes not measure membership of the point in the cones. Both getters now read the row as written: scalar coefficients, the bar product- ⟨A_ij, B_j⟩(off-diagonal counted twice, the same pairing- sdp_bar_row()uses) and, for a- PRIMAL_putqconkrow,- quad_row_value()— the helper- [cones]already uses, so there is no second sign convention to keep in step.- PRIMAL_getdualinfeasadds- max(0, -λmin(barsj[j]))for every bar block that was actually published, and since T101 the dual cone of every readable conic block (- d = s·(c + A'y)in- K*, relative to the block, with a block left unmeasured when a member carries a finite bound or the model has a quadratic term). Before this, a row whose whole content was a bar lost its left side and kept its right side as a violation: the 2×2 AM-GM model of T96 A, solved to- rel_pri = 1.7e-9, published- getprimalinfeas = 1(T97). What the getters still do not do: primal membership in a QUAD/RQUAD/PEXP/PPOW cone is the fourth figure, printed as- [cones] rel_slack=and asserted by T90, and the signs of the row multipliers are not part of- getdualinfeas. The reference's own convention for either was not read, so neither was invented.
- Solution information is exposed per index and as a summary (T113).PRIMAL_getpviolcon/getpviolvar/getpviolbarvar/getpviolconesand the dualPRIMAL_getdviolcon/getdviolvar/getdviolbarvar/getdviolconeswrite the violation of the indices listed insub, andPRIMAL_getsolutioninforeports the maxima of those pluspobj/dobj. The primal definitions are the reference's exactly: rowmax(l − a'x, a'x − u)witha'xfromrow_activity(the three doors ofT97), variable the bound/domain ofvar_violation, barmax(−λmin(X), 0), conemax(0, −cone_signed_slack), integermin(x−⌊x⌋, ⌈x⌉−x). A per-index getter reads the current model against the published point, so perturbing a bound after the solve yields a known violation (on a row 2, on a variable 4, and — for the dual — a multiplier left on a side with no bound reads 1), andgetprimalinfeas == max(per-row, per-var). A refusal (bad index, no solution, null buffer) writes nothing. Dual convention and declared scope: thedviol*values follow this solver's dual convention (see «Dual conventions», whose sign is the reference's mirror),dviolconesis the reference's exact quadratic-cone formula but this solver's signed slack for the other cone types (the page does not spell out their generalization), and a variable that is a cone member is left unmeasured rather than measured wrong — the same boundarygetdualinfeasdraws.
- The bar route publishes no bound slacks on a cone's own variables. On a task
with bars and cones the conic build answers (the cut route is skipped by
numcones > 0), and it writesslx/sux = 0for the members of a cone block where the block's real dual is a nonzero vector — measured onT101I, where the dual of(x0,x1,x2)is(5/3,-4/3,1)and the publishedslx+suxreads0, while the same model without the bar (T101D) publishes-d. The consequence is not cosmetic: the dual feasibility of a cone cannot be read off those slots, socone_dual_worstrecomputes the block fromcandyinstead, which is also why the zeros do not turn into a reported violation.
- Result-code numbers are this solver's own. PRIMAL_RES_ERR_ARG = 1001,TRM_MAX_ITER = 1007and friends collide with MOSEK's license-error codes (1001LICENSE_EXPIRED, 1007FILE_LICENSE); MOSEK's iteration limit isMSK_RES_TRM_MAX_ITERATIONS = 100000. They are left as they are because MOSEK does not report infeasibility or unboundedness as a result code at all:MSK_optimizetrmreturnsMSK_RES_OKand the verdict is read fromMSK_getprosta/MSK_getsolsta. Matching the numbers would mean matching that policy, which is an API change rather than a correction.
- PRIMAL_SOL_ITGis an alias, not a third point. MOSEK keeps the integer relaxation under- ITR/- BASand the integral solution under- ITG; this solver's branch and bound stores one point, so all three keys read it.- T94C asserts the equality, so the deviation is visible rather than silent.
- One solution item is not produced. SNX("lagrange multipliers corresponding to the conic constraints on the variables") has no counterpart here: its conic duals live per cone block rather than per variable.PRIMAL_getsolutionsliceanswersPRIMAL_RES_ERR_ARGfor it.XC("solution for the constraints") is served now: it is the row activity, read from the model (the three doors ofT97), andT94D asserts it againstx0+x1.
- A name that is not there answers PRIMAL_RES_ERR_ARG. Here "no such name" is the same code with which every other getter of this API says no — the ray getters ofT85included — so one code means "no answer". The reference is recalled to have a dedicated code for a missing name, but its response-code page was unreachable this round (curlblocked by the permission classifier,WebFetchblocked on the retry), so that recollection is not the reason for the choice and is not stated as a fact: the deviation is declared against what this code measurably returns, whichgetidxvar/getidxcondocument.
- The whole-name getters hand out borrowed pointers, not characters. The
reference's getallvarname/getallconnamefill caller-owned buffers of fixed name length; here the array elements are the task's own strings — valid untilPRIMAL_deletetask, never to be freed — which is the same ownershipgetvarnameidxalready documents. An unnamed slot reads""in both readings: the empty string is not a borrowed object, so it is asserted by content and not by address.
- The readers and the writer name the task. This deviation is closed: the
writer emits each row and column under the name the task carries
(getconnameidx/getvarnameidx/getobjname), falling back toc%d/x%d/objonly for an empty name or one containing whitespace (MPS/LP names are whitespace-delimited), and it writes the reference'sOBJNAMEsection for the objective. A read → write → read round-trip now preserves the labels as well asnumcon(T123fixed the row count,T112G asserts the names). The two tables keep their own namespaces and cannot collide: a name is unique inside its table (T102/T105), and the file puts variables and constraints in different sections.
- The bar-name namespace is our decision, not a read rule. That a bar's name is
independent of the variable and constraint tables is the coherent extension of the
rule measured for the two scalar tables in T102, applied to a table that MPS-style naming never had: the reference's rule for naming bar variables was not read (the documentation fetch stayed blocked at this session's permission level) and was not invented. What is measured is the consequence, both halves of it — one string naming a variable, a constraint and a bar at once, and the same string refused on a second bar (T105C).
- The cone and objective namespaces are the same kind of decision, and one surface is
deliberately absent. The reference's rule for naming a cone block was not read
either (the documentation fetch stayed blocked at this session's permission level), so
the independence of the cone namespace is ours, stated as such and measured as a
consequence. The gap list also names a wider surface we do not implement: the
reference addresses some entities by a (type, index) pair (getidxfunoverMSKfunctiontypee). It is absent on purpose, not by oversight: the numbering of that enum was not readable in this round, and giving it our own numbers is precisely the defectT93andT94spent two rounds correcting. A fake enum is worse than a missing getter, so the gap is declared rather than filled.
- The two whole-Qreaders answer in the shape of our counter.getqobjenumerates the triplet store (write order, stored zeros, one entry per cross term) andgetqconkthe upper triangle of the dense symmetrized store — the very tablesgetnumqobjnz/getnumqconknzcount, sonumretcannot disagree with the number the user was able to ask for. The reference's shape for these two readings was not read (the documentation fetch stayed blocked at this session's permission level) and was not invented: the two*ijgetters answer the same value on both halves of a cross term either way, andx'Qxmultiplies it by 2, so which list is handed out does not change the operator that gets solved.
- Two contracts, one store: the store reads back entries, the scalar reads
back the operator. An (i,j)written twice is stored twice —putarowclears the row and appends, nothing deduplicates — and the store-reading surface says so:getnumanzscounts both entries andgetarow/getacolhand both back. The scalar read of one(i,j),getaij, answers the operator instead: it returns the sum of the stored entries, which is the coefficient every route that solves the model actually used. It used to read the first, so the number describing the answer contradicted the answer (measured:putarow(0,{0,0},{1,2})with3x <= 3solves atx=1on the coefficient 3 while the getter said 1). Chosen on parity, not preference:getqobjijalready sums andputqconkaccumulates into the dense matrix sogetqconkijalready summed —getaijwas the only exception. The reference's write-time policy was read (2026-09-25):MSK_putarow/putacolclear and then assign in sequence (a_{i,subi[k]}=vali[k]), so a repeated index keeps the last value (no summingappendrow); we keep the divergence — storage unchanged — because our standard form sums, sogetaijstays the operator.T107asserts the binding invariantgetaij(i,j) == Σ getarow(i) entries == Σ getacol(j) entriesover the whole matrix, plus that the summed number equals neither stored entry.
PrimalSolver is a research-grade convex solver, not a drop-in replacement for the industrial LP/MIP engines — but the combination of what it solves is rare among the permissively-licensed open solvers, and that combination is its niche. This is a comparison of capability, licence and accuracy class, not of speed.
What no other permissively-licensed open solver covers — said once: SDP + exponential/power cones + MIP + an interior-point (not first-order) method, in dependency-free C99.
- Clarabel (Apache-2.0) is the closest open conic IPM, but it has no PSD cone.
- SCS (MIT) has PSD and exp cones, but is first-order — accuracy around
1e-4/1e-5, not the1e-8an interior-point method reaches.
- CVXOPT / SDPA / Sedumi have SDP, but are GPL (viral for commercial use).
- HiGHS is fast on LP/MIP, but has no SDP and no non-symmetric cones.
The niche that follows: embedded/enterprise code that wants SDP, exponential cones and integer variables in C, with no dependencies and no licence server, and a migration path from a MOSEK-style API (the 48 C ports of MOSEK examples are reproduced with identical optimal values). Outside that niche — large LP/MIP at industrial scale — it is not (yet) a substitute.
PrimalSolver ports the published MOSEK 11.0 examples to its own API and reproduces
their optimal values — 48 C ports covering 47 distinct MOSEK examples (the
acc1 example is ported from both the Julia API and the C tutorial, with two
different models). Each port is independently verified: hand-derived values, KKT
via public getters, exhaustive enumeration, or a reference algorithm. A selection:
This is compatibility validation against the public examples: it demonstrates that PrimalSolver solves those problems and reproduces the expected optimal values. It is not a claim of equivalence with, or of industrial robustness comparable to, MOSEK on arbitrary large-scale problems. The example models are MOSEK's (Copyright (c) MOSEK ApS); PrimalSolver is an independent implementation that re-solves the same problems.
Beyond the MOSEK ports, samples/ (root) contains larger, deterministic
examples. Most stress the engine across every problem class: each one
checks its own result (strong duality, an exact dynamic program, a closed form, or
an independent numerical method) and reports its size and solve time, so it can
also be used as a scaling probe (./out/samples/<name> [args]). Others probe an
API shape or an algorithm rather than a size — barqcqp puts symmetric bars in
the same model as quadratic rows/objectives, lyapunov_roa and
sos_m1_certificate build SOS certificates, secure_ee_sdma runs the
Dinkelbach + SCA outer loop of rezarhp/Secure-EE-SDMA on a K=1/J=1 downlink
(rotated SOC, power SOC and exponential cone in one model), and
logcontrast_lasso casts audreyolmsted/admm-portfolio's log-contrast LASSO
logistic regression as exponential cones and compares the interior-point answer
with the repository's ADMM (Newton + joint proximal); quantum_separability
solves the PPT white-noise robustness SDP of y1-zhu/quantum-correlations over
a 9x9 semidefinite bar (checked against the closed form) and the CCNR
realignment criterion; binary_quadratic (Shor SDP), exact_cover (MIP), pwl_convex
(rotated cone regression), equilibrium (hanging masses, SOCP), dist_robust
(Wasserstein DRO, LP), facility_location (k-disk cover, MISOCP),
wasserstein (barycenter, LP), mle_density (log-concave MLE, exp cone), rank_one (best-subset MIP with
DJC), kmeans (MISOCP+DJC),
filter_design (trigonometric-polynomial SDP; its n=3 case needed the
multi-block PSD rescue, see below),
f_sparc (subcarrier/power allocation, MIP+exp),
hard_uncertain (robust uncertain inequality on box vertices, exp cones),
transformer_design (geometric program, exp cones), surface_cycles (shortest cycle
in a homology class), gp_toolbox (geometric program) and min_enclosing_ball
port MOSEK Tutorials; lovasz_theta solves the max-clique upper bound SDP
of Rudolfovoorg/Improving_Upper_Bounds_of_MCP_using_Reduction_Rules (Lovasz
theta of the complement) and checks it against sqrt(5) on C5:
The default sizes are chosen to complete in seconds; larger sizes can be passed
on the command line (e.g. ./out/samples/transport 150 150). Large LPs
(m·n > 1.5e6, m ≤ 2000) are routed automatically to the sparse
normal-equations interior point — e.g. transport 150 150 (22500 vars) solves
in ~0.07 s vs ~12.7 s for the dense simplex, and transport 300 300 (90000
vars) in ~0.61 s. lp_large (Ax=b) and lp_ineq_large (Ax≤b) are generic
sparse LPs with a planted KKT optimum, so their exact objective is known
independently of the solver; both route to the sparse IPM at default size
(~0.06 s for 6000 vars). Beyond the sparse-LP path, the dense QP/conic engine
and the outer-approximation round budgets remain the limiting factors — the
frontier described under Performance.
samples/finance/ collects portfolio-optimisation models, several following
Schmelzer, Hauser, Andersen and Dahl, "Regression techniques for Portfolio
Optimisation using MOSEK" (arXiv:1310.3397). Each is deterministic and checks
its own result:
All code is original, written from scratch in C99.
PrimalSolver solves LP, QP, SOCP, SDP, exponential/power-cone and mixed-integer problems. The algorithmic work still open (the maintained task list lives outside the repository):
- Basis warm start — no real crash basis, and solvebasis/solvewithbasisre-solve from scratch when the given basis does not measure instead of re-optimizing from it.
- I/O formats — TASK and PTF have no reader/writer (readdataformatrefuses both; MPS, LP, OPF and CBF are supported).
- Sparse augmented system in sdp.c— the unified conic KKT is assembled densely (Sys[Nsys²]), which caps the problem size.
- Symmetric KKT in socp.c— the SOCP KKT is assembled in an unsymmetric form, so the sparse symmetricLDLᵀcannot replace the sparseLUthere yet.
- Parameter coverage — sixteen parameters are wired against the reference's hundreds (declared deviation; the table is declarative, so adding rows is cheap).
- Exp/power polish floor — on logistic_largeno iterate enters the declared relative triple, so the tangent cuts answer, and the native route itself can move with the optimization level (the KKT is conditioned past1/eps).
- A conic extreme not attained inside a bar has no recession direction and is reported as "no verdict" rather than a value.