Parallel addJMInvJtLocal overflows columnResult buffer when constraint rows > DOFs (parallelInverseProduct)

[AdoHaha]

Problem

Description

MatrixLinearSolver<Matrix,Vector>::addJMInvJtLocal computes the compliance product J M⁻¹ Jᵀ in parallel and stores each result column in a per-row scratch vector that is sized with the wrong dimension:

sofa::type::vector<Vector> columnResult(J->rowSize());
...
columnResult[row].resize(J->colSize());   // <-- must be J->rowSize()
...
columnResult[row][row2] += acc;           // row2 iterates over constraint rows
result->add(row2, row, columnResult[row][row2]);

J is created by buildComplianceMatrix (same file) as

j_local->resize(result->rowSize(), l_linearSystem->getSystemMatrix()->colSize());

so J->colSize() is the number of mechanical DOFs and J->rowSize() is the number of constraint rows (= result->rowSize()). columnResult[row] is column row of W, therefore it must have J->rowSize() entries. Whenever the constraint row count exceeds the DOF count (J->rowSize() > J->colSize()) — the normal situation for contact-rich scenes, and already when a small mechanical model carries contact constraints — every column computation writes out of bounds.

The bug is deterministic (not a race): it occurs on the first compliance build. In release builds (NDEBUG, i.e. shipped installs) the out-of-bounds write is silent and is usually reported later as malloc(): corrupted top size / munmap_chunk(): invalid pointer on an unrelated allocation, or as wrong W entries. In debug builds FullVector::checkIndex asserts with invalid index.

The serial path (singleThreadAddJMInvJtLocal) is correct because it reuses the system RHS/solution vectors sized by the system matrix.

This is still present on master (checked at d43dc1927fb25d7ad6cca00703b2e4b6cd44d466, line 308) as well as in the v26.06.00 and v25.12.00 releases.

Steps to reproduce

A self-contained script and logs are available in this gist:
https://gist.github.com/AdoHaha/6ba87f1b18e194ad75722b1cce7ab199

The scene is deliberately tiny in DOFs and rich in constraint rows:

• FreeMotionAnimationLoop + BlockGaussSeidelConstraintSolver + FrictionContactConstraint(mu=0.5)
• /rigid: one Rigid3d MechanicalObject → 6 mechanical DOFs (J->colSize() == 6), EulerImplicitSolver, native EigenSimplicialLDLT(template=CompressedRowSparseMatrixd, parallelInverseProduct=true), GenericConstraintCorrection, mapped SphereCollisionModel
• /probes: static node (no solver) with N points and a SphereCollisionModel, placed just inside first contact with the rigid sphere. With mu=0.5 each contact pair yields 3 constraint rows, so J->rowSize() == 3N.

1. Download sofa_parallel_compliance_overflow_repro.py and run_repro.sh from the gist into one directory.
2. Set up a v26.06 (or v25.12) runtime, e.g.:
   SOFA_ROOT=<install>, PYTHONPATH=<install>/plugins/SofaPython3/lib/python3/site-packages, LD_LIBRARY_PATH=<install>/lib.
3. Run (see run_repro.sh for the exact env):

   python3.12 sofa_parallel_compliance_overflow_repro.py --spheres 3 --steps 1 --threads 2

Observed results:

--spheres	constraint rows	mechanical DOFs	outcome
1	3	6	OK
2	6	6	OK (boundary; indices 0…5 still in range)
3	9	6	malloc(): corrupted top size → SIGABRT
10	30	6	malloc(): corrupted top size → SIGABRT
10, --serial	30	6	OK (serial path)

Crash backtrace with 3 spheres (9 rows vs 6 DOFs):

[repro] spheres(contact rows)=3 mechanical DOFs=6 parallel=1 threads=2
malloc(): corrupted top size

munmap_chunk(): invalid pointer
########## SIG 6 - SIGABRT ##########
  sofa::helper::BackTrace::sig(int)
  pthread_kill
  gsignal
  abort
  __libc_free
  void Eigen::internal::permutation_matrix_product<...>::run<...>(...)
  sofa::component::linearsolver::direct::EigenDirectSparseSolver<double, ...>::solve(..., FullVector<double>&, FullVector<double>&)
  sofa::component::linearsolver::MatrixLinearSolver<...>::addJMInvJtLocal(...)::{lambda(int)#1}::operator()(int) const
  sofa::simulation::WorkerThread::runTask(sofa::simulation::Task*)
  sofa::simulation::WorkerThread::doWork(sofa::simulation::Task::Status*)
  sofa::simulation::WorkerThread::run()

The corruption is detected inside the worker-thread lambda that performs the parallel compliance product.

For a real-scale instance of the same crash: a scene with 82 wall vertices (246 mechanical coordinates) and 940 contact rows aborted with the same malloc(): corrupted top size when parallelInverseProduct=true was enabled on EigenSimplicialLDLT.

Expected behavior

parallelInverseProduct=true should compute the same compliance product as the serial path, without memory corruption, for any J->rowSize() > J->colSize(). Fix:

-            columnResult[row].resize(J->colSize());
+            columnResult[row].resize(J->rowSize());
+            columnResult[row].clear();

(The explicit clear() keeps the += zero-initialization explicit; currently it only works because resize is called on a freshly constructed vector.)

Note: SparseLDLSolver is not affected because it overrides the product with its own buffers. A regression test can be added cheaply at the unit level by calling addJMInvJt with a small system matrix and a SparseMatrix<Real> where rowSize() > colSize(), run under ASan/valgrind. A short matched simulation run is not a sufficient acceptance test for this bug (a 0.05 s run passed with the buggy native path because the corruption happened not to be fatal in that window).

---

Environment

Context

• System: Ubuntu, kernel 6.8.0-134-generic, x86_64
• Version of SOFA: v26.06.00 (7c18e95d5c5f2839079892c69e7d89a313c79603, origin/v26.06); also verified in v25.12.00 (4d6c2b8e8cd976836c5e071e9db14f1963d84cc9) and still present on master (d43dc1927fb25d7ad6cca00703b2e4b6cd44d466)
• State: install directory (release build, NDEBUG)

Command called

python3.12 sofa_parallel_compliance_overflow_repro.py --spheres 3 --steps 1 --threads 2

Env vars

--- sys.version ---
3.12.3 (main, Aug 31 2026, 10:18:26) [GCC 13.3.0]
--- SOFA_ROOT ---
<install>/SOFA_v26.06.00_Linux
--- uname ---
Linux igor-z 6.8.0-134-generic #134-Ubuntu SMP PREEMPT_DYNAMIC x86_64
--- relevant settings ---
EigenSimplicialLDLT.parallelInverseProduct = true
Sofa.Core.TaskScheduler.init(2)

---

Logs

Full output

[repro] spheres(contact rows)=3 mechanical DOFs=6 parallel=1 threads=2
malloc(): corrupted top size

munmap_chunk(): invalid pointer
########## SIG 6 -
########## SIG 6 - SIGABRT: usually caused by an abort() or assert() ##########
  sofa::helper::BackTrace::sig(int)
  pthread_kill
  gsignal
  abort
  __libc_free
  void Eigen::internal::permutation_matrix_product<Eigen::Map<Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0, Eigen::Stride<0, 0> >, 1, false, Eigen::DenseShape>::run<...>(...)
  sofa::component::linearsolver::direct::EigenDirectSparseSolver<double, ...>::solve(..., FullVector<double>&, FullVector<double>&)
  sofa::component::linearsolver::MatrixLinearSolver<...>::addJMInvJtLocal(...)::{lambda(int)#1}::operator()(int) const
  sofa::simulation::WorkerThread::runTask(sofa::simulation::Task*)
  sofa::simulation::WorkerThread::doWork(sofa::simulation::Task::Status*)
  sofa::simulation::WorkerThread::run()
timeout: the monitored command dumped core

Complete logs (3spheres crash, serial 10spheres pass) and the full script are attached in the gist linked above. The serial control with 10 spheres (30 rows vs 6 DOFs) completes and produces a 30-entry multiplier vector.

Content of build_dir/CMakeCache.txt

N/A — this report uses a pre-built install directory, not a local build.

---

Thank you for your report.