feat: add MATLAB LDL linear system solver backend (#29)

* feat: add MATLAB LDL linear system solver backend

Add a new linear system solver backend that uses MATLAB's built-in
sparse ldl() factorization (MA57) instead of the bundled QDLDL solver.
This may be faster for large problems due to MA57's supernodal
factorization and advanced fill-reducing orderings.

Usage: set pars.matlab_ldl = true before calling scs() or scs_init().

The backend implements the standard linsys.h interface via mexCallMATLAB:
- Init/update: forms KKT matrix in C, factorizes via ldl(K,'vector')
- Solve: calls MATLAB helper for triangular solves with stored factors
- Factors are kept as persistent mxArrays across MEX calls

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: do LDL solve inline in C instead of via MATLAB helper

Eliminate the private/scs_matlab_ldl_solve.m helper function to avoid
potential mexCallMATLAB scoping issues with the private/ directory.

The solve now operates entirely in C with three mexCallMATLAB calls
(L\bp, D\y, Lt\z) and does the permutation/unpermutation in C.
L' is pre-computed and stored at factorization time.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: simplify MATLAB LDL backend to use direct backslash solve

Replace the manual LDL factorization + 3-step mldivide solve with
a direct K\b approach using the full symmetric KKT matrix. This
eliminates several potential issues:
- mexCallMATLAB potentially modifying persistent factor arrays
- Manual permutation handling
- Separate L\, D\, L'\ steps

The KKT matrix (stored as upper triangle in CSC) is now symmetrized
via K_sym = K + triu(K,1)' and stored as a persistent mxArray.
Each solve simply calls mldivide(K_sym, b).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* test: run all tests with matlab_ldl solver alongside direct and indirect

Replace the boolean `use_indirect` test parameter with a three-way
`solver` parameter ('direct', 'indirect', 'matlab_ldl') across all
test files. This ensures the MATLAB LDL backend is exercised by every
test in the suite: LP, QP, SOC, SDP, exponential, power, mixed cones,
infeasible, unbounded, workspace, info fields, defaults, and string
params.

Remove the standalone matlab_direct.m test file since its coverage
is now subsumed by the parameterized tests.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* feat: make MATLAB LDL the default linear system solver

The MATLAB built-in LDL backend (via backslash) is now the default
solver, replacing QDLDL. Users can explicitly select QDLDL with
pars.use_qdldl = true. Test parameters updated accordingly: 'default'
(matlab_ldl), 'qdldl', and 'indirect'.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* perf: cache LDL factors in C, eliminate per-iteration refactorization

Instead of calling MATLAB's mldivide (backslash) every ADMM iteration
(which refactorizes from scratch each time), we now:

1. Call ldl(K, 0, 'vector') once at init and on diag_r updates only
2. Extract L, D, permutation into C arrays
3. Do pure C forward/backward triangular solves each iteration

The thresh=0 argument disables Bunch-Kaufman pivoting so D is purely
diagonal (no 2x2 blocks), matching QDLDL's structure. The upper-
triangular KKT is passed directly to ldl (which only reads the upper
triangle), eliminating the symmetrization step.

Also adds memory_leak test that runs scs 100 times in a loop to verify
no persistent mxArray leaks.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: use ldl(K,'vector') without threshold, handle 2x2 D blocks

The threshold parameter in ldl(K, thresh, 'vector') is not supported
for sparse matrices (only dense). Removed it and properly handle the
tridiagonal D matrix from Bunch-Kaufman pivoting, which can have 2x2
blocks for indefinite KKT systems.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* test: add diagnostic test to verify LDL extraction and solve logic

Pure MATLAB test that replicates the C-level LDL extraction (strip unit
diagonal from L, extract D diagonal/sub-diagonal) and solve (forward,
block-diagonal, backward, with permutation). Compares against backslash.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: symmetrize KKT before calling ldl (sparse ldl needs full symmetric)

MATLAB's sparse ldl() documentation says "only the upper triangle is
accessed" but this means for efficiency on symmetric input — it does NOT
reconstruct a symmetric matrix from upper-triangular-only input. The
diagnostic test (ldl_diag) confirms: ldl on triu(K) gives wrong factors,
while ldl on the full symmetric K works correctly.

Added symmetrize_upper back to the factorization path. This runs only
at init and on diag_r updates (not every iteration), so the overhead
is minimal. The per-iteration solve remains pure C.

Updated ldl_diag test to verify:
1. C-level solve algorithm is correct on symmetric input
2. ldl on upper-tri-only gives wrong results (documents the limitation)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* test: improve ldl diagnostic to show exactly what ldl does with triu input

Prints whether L*D*L' matches K_upper(p,p) (factoring asymmetric) or
K_full(p,p) (correctly reading upper triangle), plus any warnings.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* refactor: build symmetric sparse matrix in C instead of MATLAB ops

Replace the 3 mexCallMATLAB calls (triu, ctranspose, plus) with a pure C
function that constructs the full symmetric MATLAB sparse matrix directly
from the upper-triangular CSC. For each off-diagonal entry (i,j) with
i < j, both (i,j) and (j,i) are written. Row indices are sorted per
column via insertion sort (columns are small in practice).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* test: add benchmark and cross-validation tests, document LDL limitations

Add comments in private.c explaining:
- Why full symmetric matrix is required for MATLAB's sparse ldl
  (structurally asymmetric input silently produces wrong factors)
- Why AMD permutation is recomputed on every refactor (MATLAB's ldl
  wraps MA57 but doesn't expose separate symbolic/numeric phases)

New test files:
- test/benchmark.m: head-to-head matlab_ldl vs qdldl across problem
  sizes, sparsity patterns, and LP/QP types with summary table
- test/matlab_ldl.m: cross-validation ensuring both solvers agree
- test/ldl_diag.m: extended with 2x2 block, QP, and multi-RHS tests

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: construct feasible problems in cross-validation and benchmark tests

Random b vectors don't guarantee feasible/bounded problems. Use the
standard pattern b = A * x_feas + s_feas with s_feas in the cone
interior (LP: s > 0, SOC: s(1) > ||s(2:end)||).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: benchmark should verify solvers agree, not require all problems solved

Random LPs can be genuinely unbounded (especially with m=2n). The
benchmark compares solvers against each other, so verify they return
the same status and only compare solutions when both solve.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: ensure test problems are both primal and dual feasible

Construct c = A'*y_feas (y_feas in cone) to guarantee boundedness,
not just primal feasibility. Use AbsTol for cross-solver comparisons
since RelTol fails on near-zero values.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: correct SCS dual sign convention (c = -A'*y for dual feasibility)

SCS dual constraint is A'y + c = 0, so bounded LP requires
c = -A'*y_feas, not c = A'*y_feas. Also relax benchmark AbsTol
to 1e-2 for sparse ill-conditioned problems.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* revert: restore benchmark AbsTol to 1e-3

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* docs: document matlab_ldl as default solver backend

Update README and SCS docs to reflect that the default linear system
solver is now MATLAB's built-in sparse LDL (MA57), and document
use_qdldl as an explicit alternative for the bundled QDLDL solver.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* revert: restore scs submodule pointer to upstream master

The docs change is in a separate PR (cvxgrp/scs#382). Keep the
submodule pointer at upstream master so this PR can merge independently.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* test: increase benchmark problem sizes for meaningful comparison

Previous sizes (n up to 200) were too small — overhead dominated and
factorization cost was negligible. Increase to n up to 2000 with
sparse densities (0.01-0.1) which better represents real problems.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* feat: make QDLDL the default solver, matlab_ldl opt-in via use_matlab_ldl

Benchmarks show QDLDL is generally faster due to the mexCallMATLAB
overhead in the matlab_ldl pipeline. Switch default back to QDLDL
and rename the setting from use_qdldl to use_matlab_ldl.

Update all test parameterizations from 'qdldl' to 'matlab_ldl' and
update README/benchmark labels accordingly.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* fix: address code review — memory safety, error handling, sort removal

1. Update scs_matlab_direct.m docs: per-iteration solves run in C,
   no MATLAB callbacks during iterations
2. Use mexCallMATLABWithTrap instead of mexCallMATLAB to prevent
   hard-abort on ldl() errors (e.g., out of memory) which would
   leak scs_calloc'd memory
3. Check extract_L return value in matlab_ldl_factor and propagate
   errors instead of silently continuing with a null/partial L
4. Add allocation checks for L->p, L->i, L->x in extract_L
5. Remove redundant insertion sort in scs_to_mxsparse_symmetric —
   row indices are already sorted by construction (upper-tri entries
   have rows <= col, mirror entries have rows > col)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* feat: make matlab_ldl the default solver, qdldl opt-in via use_qdldl

Switch the default MATLAB solver backend from QDLDL to MATLAB's
built-in sparse LDL factorization (MA57). Users can opt into QDLDL
with settings.use_qdldl = true. Updates README, SCS docs, all tests,
and benchmark to reflect the new default.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com>

Author: bodono

README.md

@@ -59,12 +59,14 @@ scs_finish(work);                        % free workspace
 
 ### Solver backends
 
-By default SCS uses the sparse direct (LDL) solver. Alternatives:
+By default SCS uses MATLAB's built-in sparse LDL factorization (MA57 under
+the hood). Alternatives:
 
 ```matlab
-settings.use_indirect = true;  % conjugate gradient
-settings.dense = true;         % dense Cholesky (for dense A)
-settings.gpu = true;           % GPU solver
+settings.use_qdldl = true;      % bundled QDLDL sparse direct solver
+settings.use_indirect = true;    % conjugate gradient (iterative)
+settings.dense = true;           % dense Cholesky (for dense A)
+settings.gpu = true;             % GPU solver
 ```
 
 ### Cones

make_scs.m

@@ -71,6 +71,7 @@
 compile_direct(flags, common_scs);
 compile_indirect(flags, common_scs);
 compile_dense(flags, common_scs);
+compile_matlab_direct(flags, common_scs);
 if (gpu)
     compile_gpu(flags, common_scs);
 end

private/compile_matlab_direct.m

@@ -0,0 +1,4 @@
+function compile_matlab_direct(flags, common_scs)
+% compile MATLAB LDL direct solver (uses MATLAB's built-in ldl())
+cmd = sprintf('mex -O -v %s %s %s %s COMPFLAGS="$COMPFLAGS %s" CFLAGS="$CFLAGS %s" -Iscs -Iscs/linsys -Iscs/include -Iprivate/matlab_linsys private/matlab_linsys/private.c %s %s %s %s -output scs_matlab_direct', flags.arr, flags.LCFLAG, flags.INCS, flags.INT, flags.COMPFLAGS, flags.CFLAGS, common_scs, flags.link, flags.LOCS, flags.BLASLIB);
+eval(cmd);