Direct solver on Apple Accelerate: 120,000 DOF in 3.8 s instead of 87

A direct solver is the honest way to solve a contact problem: it gives an exact answer, does not stall on ill-conditioned constraint equations, and never quietly returns a half-converged field. Its problem has always been cost — factorization time and memory grow far faster than the model does. On the desktop build we now hand that factorization to Apple's Accelerate sparse Cholesky, and the cost dropped by a factor of 20 or more.

What changed

The native C++ core extracts the lower triangle of the stiffness matrix in CSC form, runs a symbolic analysis with AMD reordering, and factorizes through Accelerate. On Apple Silicon that spreads the dense supernodes across the performance cores and the matrix accelerators, which is exactly the part of a Cholesky factorization that dominates the clock.

Measured against the previous direct solver (Eigen SimplicialLDLT) on 3D hex models of beams and assemblies, same machine, same matrices:

2,940 DOF

41.2 ms → 2.9 ms

9,720 DOF

337 ms → 18.6 ms

54,000 DOF

24.8 s → 0.95 s

120,000 DOF

87.2 s → 3.8 s

Speed-up

14× – 26×

Agreement with the reference

1e-12

The displacement fields match the old solver to 1e-12 — machine zero for this problem. Nothing about the formulation changed; only the factorization backend did.

A higher ceiling

Because the factorization got cheap, the threshold at which the product still chooses the direct path moved from 120,000 to 300,000 degrees of freedom on macOS. Below that line you get an exact solution with no convergence monitor to watch. Above it, the solver switches to the iterative path as before.

The cascade, and when we refuse to use AMG

The solver is a cascade rather than a single method: a direct sparse factorization for small and medium systems, an algebraic multigrid (AMG) preconditioner for large clean meshes, an incomplete-Cholesky conjugate gradient for ill-conditioned contact systems, and a basic CG at the bottom. What we learned this month is that the choice matters more than the tuning.

On an assembly with exact multipoint constraints — 74,262 reduced degrees of freedom — AMG did not merely converge slowly. The residual grew: 4.06 at iteration zero, 34.4 by iteration 480, with an empty convergence chart in the UI and no answer at the end. The same payload through the direct solver: about 3 s, converged. The cause is structural — after the constraint reduction the near-nullspace the multigrid relies on belongs to the unreduced solid, not to the reduced matrix the smoother actually sees.

This is not an indictment of multigrid. On a clean first-order mesh with no constraint equations, the same AMG converges in 8 iterations at 6 million degrees of freedom. So the fix was not a better preconditioner but an explicit rule: when the model uses exact MPC reduction, the cascade skips AMG outright, and a guard test fails the build if that rule is ever removed. There is also a hard stop — if the residual exceeds twice its starting value after 15 multigrid cycles, the solver gives up instead of burning minutes on a diverging run.

What this buys on a real model

A frame of 88 bodies tied through 1,106 contact pairs — 3,318 constraint equations — now solves as one monolithic structure, with a maximum deflection under self-weight of 0.126–0.143 mm. Before the constraint assembly was fixed, the same model produced deflections of tens of metres and 87 of the 88 bodies flew apart.

Where the limits still are

Honest numbers, because the shape of the curve matters more than any single measurement: factorization work grows roughly with the square of the problem size, and the fill-in of the factor grows faster than the matrix itself. Measured on a Linux benchmark VM, a 121,680-DOF model takes 87.2 s and 1.17 GB; at 242,760 DOF it is about 14 minutes and 3.8 GB. Past that the direct path is no longer interactive, and beyond a million degrees of freedom on assemblies with exact constraints we currently have no good route at all — the iterative prototypes we tested (deflated PCG, a Schur complement formulation) all failed to converge on a 1.1M-DOF system. That gap is open, and we would rather say so than quote a number we cannot reproduce.

The figures above were measured on Apple Silicon for the Accelerate path and on an 8-vCPU Linux VM for the scaling study. They are not interchangeable, and we do not mix them in one table.