I made my code 2.4x faster with this simple trick

27 October 2026

A copper bar hitting a wall took 13.85 s to simulate on one core. Now it takes 5.74 s, so 2.41× faster. The final length is 19.419 mm before and after.

Two changes did it. A closed-form solver for 3×3 eigenvalue problems, and a plasticity kernel that works on vertices in batches. Both live in the same build and I can switch each one on or off at runtime. So every comparison below is the same binary on the same inputs.

What BALEK does

BALEK is my explicit solid-dynamics code. Bars hitting walls, bodies colliding, that kind of thing. It’s written in C on top of PETSc and MPI. Explicit means there is no big system of equations to solve. The solver takes tiny time steps, thousands per run, and updates every point of the mesh at each one.

The hot loops are compute-bound. So any gain has to come from the maths done at each point, and the heaviest maths in a metal impact is plasticity.

Why a crash simulation needs eigenvalues

Metal behaves elastically (it springs back) until it’s stressed past its yield limit. Past that it flows, and the deformation stays. That’s plasticity. So at every time step, for every point, the solver roughly does this:

1. Assume the material just stretched elastically. That gives a trial state.
2. Check if that trial state is more than the metal can carry.
3. If yes, pull it back onto the yield limit and record how much permanent (plastic) deformation happened. This is the return mapping.

With large deformations, the elastic stretch at a point is a symmetric 3×3 matrix, call it b. Step 3 is much easier in the matrix’s own axes. Any symmetric 3×3 matrix can be written like this:

b = Q Λ Qᵀ

Λ holds three numbers, the eigenvalues: how much the material is stretched along three perpendicular directions. Q holds those directions, the eigenvectors. You rotate into those axes, the plasticity update becomes an easy problem on three numbers, and you rotate back.

Problem is, this happens at every point, at every stage of every time step. That’s millions of small decompositions in a short run.

The test case

A copper Taylor bar. G.I. Taylor came up with it in 1948 and it’s a very blunt test. You fire a short metal cylinder flat end first into a rigid wall, and you measure what’s left. The front squashes into a mushroom while the back is still travelling, and the final shape tells you how strong the metal is at high strain rate.

The strength depends on strain rate and temperature (a Johnson-Cook model). I run it on a coarse mesh of 475 vertices, in two ways:

• ALE. The mesh is allowed to slide relative to the metal so it doesn’t get mangled.
• Total Lagrangian. The mesh moves with the metal and everything is measured against the starting shape.

Plus a two-body version of the impact with 950 vertices and a simpler hardening law.

How I measured

All on a 16-thread desktop. I ran the jobs side by side, one per physical core, in six rounds. Between rounds I rotated which job got which core (a Latin square), so no level gets the lucky core. Each figure below is the median of the six.

The spread between rounds was 0.1 to 0.6 %. I only count a gain when it beats twice that spread, or 2 %, whichever is bigger. Every step below clears it easily.

Each run alone on an idle machine tells the same story. The total Lagrangian bar goes 13.72 s → 8.85 s → 5.61 s that way.

Level 0: where I started

A while back the profile showed 3×3 eigendecompositions eating most of the per-step work. I replaced a general-purpose library routine with a dedicated iterative 3×3 solver. That iterative solver is level 0 here, the oldest path I can still switch on.

Level 1: closed-form eigensolver, 22 to 35 % off

An iterative solver repeats rotations until the matrix is diagonal. A closed-form one computes the eigenvalues directly from a formula. At n = 3 that’s a short fixed sequence of operations, with no loop to converge.

The tricky part is accuracy. Impacts produce lots of states with two or three equal eigenvalues (uniaxial compression, near-hydrostatic), and that’s exactly where naive 3×3 formulas fall over. So the closed-form path takes extra care there.

The eigensolver is the only thing that changes between level 0 and level 1. So turns out on the total Lagrangian bar the iterative solver was at least 35 % of the whole run.

The closed-form solver is now the default.

Level 2: batches of vertices, another 6 to 36 %

The plasticity kernel used to handle one vertex at a time. One vertex per call leaves the compiler little room to vectorise. The new block kernel takes vertices in small batches. The same arithmetic then runs on several vertices at once with SIMD instructions.

The ALE number is smaller for a boring reason. On the ALE mesh only one of the two places that call the plasticity update uses the block kernel so far. On the total Lagrangian mesh both do.

I also kept a reference version of the block kernel. The fast one matches it bit for bit in every output frame.

All three levels side by side:

(Wall time per time step. The step counts agree within 6 steps at every level, so basically all of the speedup is work saved per step.)

A funny side effect is that the total Lagrangian bar used to be the slow one. Now it beats the ALE bar by 2.5 s, even though it takes almost twice as many time steps (17,186 against 9,717).

What the profile said

BALEK has its own phase timers. On the total Lagrangian bar, level 0 against level 2:

• The material update went from 7.36 s to 1.92 s, so 3.8× faster.
• The plasticity commit went from 3.30 s to 0.60 s, so 5.5× faster.
• The rest of the right-hand side barely moved, 2.50 s against 2.46 s.

So at level 0 the material update and the commit were 10.66 s of a 13.80 s run, about three quarters. At level 2 they’re 2.52 s of 5.61 s. The parts I didn’t touch moved by 0.04 s at most.

On the ALE bar the right-hand side went from 8.96 s to 7.14 s and the commit from 1.96 s to 0.38 s.

Checks after every level

Of course a speedup only counts if the answer stays put.

The final length agrees within 0.00024 mm and the radius within 0.00071 mm. The peak temperature agrees to two parts in 100,000.

Energy and mass are conserved at every level. Energy drift sits between 10⁻¹⁶ and 3 × 10⁻¹⁵, and mass drift is at most 4.4 × 10⁻¹⁶. That’s double-precision round-off. (On the total Lagrangian mesh mass is exact by construction, which is a bit of a cheat.)

The QA checkers pass at every level.

Same source, different bits

Fun one from when I first wrote the closed-form solver. It changed the output of the iterative path, the routine I hadn’t touched at all.

I swapped old and new object files one at a time until it came down to a single file, and the disassembly showed what happened. The new call sites made the compiler stop inlining the iterative routine. The out-of-line copy, under fast-math-style optimisation, rounds slightly differently. Floating-point addition isn’t associative, so when the compiler is allowed to reorder things, different machine code can mean different last bits. Same source, different bits 😅

I moved the new solver into its own file, the old machine code came back identical, and the iterative path reproduced its references exactly again.

Why the block kernel is opt-in for now

The closed-form solver is now the default. The block kernel stays opt-in until it has gone through the full QA bench. I also want it in the second ALE call site first, since that’s where the ALE bar is still leaving time on the table.

Until then it’s one switch, set at runtime.

Bref

• Same build, same inputs, three levels. The total Lagrangian Taylor bar went from 13.85 s to 5.74 s, so 2.41×. The ALE one is 1.41× and the two-body impact 1.37×.
• Closed-form eigensolver: 22 to 35 % off whole runs. Now the default.
• Batched plasticity kernel: another 6 to 36 %. Opt-in while it goes through the full QA bench.
• Final length and radius agree within a few ten-thousandths of a mm. Energy and mass stay at round-off.
• Next: the block kernel in the ALE right-hand side, then decide if it becomes the default.

References

• G.I. Taylor (1948), “The use of flat-ended projectiles for determining dynamic yield stress I. Theoretical considerations”, Proc. R. Soc. Lond. A 194, 289-299. doi:10.1098/rspa.1948.0081
• G.R. Johnson & W.H. Cook (1983), “A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures”, Proc. 7th International Symposium on Ballistics, The Hague, 541-547.
• J.C. Simo & T.J.R. Hughes (1998), Computational Inelasticity, Springer. doi:10.1007/b98904
• D. Goldberg (1991), “What every computer scientist should know about floating-point arithmetic”, ACM Computing Surveys 23(1), 5-48. doi:10.1145/103162.103163
• PETSc