Support Eigen 5: reimplement removed internal::make_coherent - #1
Closed
zfergus wants to merge 1 commit into
Closed
Conversation
Eigen 5.0 removed Eigen::internal::make_coherent from unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h, which AutomaticDifferentiation.hh references in the pow(AutoDiffScalar, AutoDiffScalar) overload. Because the call is qualified, name lookup happens at template definition time, so any TU including this header fails to compile against Eigen >= 5 even if pow is never instantiated (this also breaks all of MeshFEMSparse, whose SparseMatrices.hh includes this header). Reimplement it with the Eigen 3.4 semantics behind a version guard. Note Eigen 5 moved to semantic versioning: EIGEN_WORLD_VERSION remains 3 forever and the new major version lives in EIGEN_MAJOR_VERSION.
This was referenced Jul 31, 2026
Merged
Contributor
|
Thanks! Before noticing this PR, I fixed the issue in a different way via this commit. Rather than reimplementing the legacy Incidentally, I encountered this issue when integrating our code into PolySolve/PolyFEM (following your recent Eigen 5 update). The pull requests are here and here in case you're curious. |
zfergus
added a commit
to ipc-sim/ipc-toolkit
that referenced
this pull request
Sep 2, 2026
* Add contact assembly benchmarks (Phase 0 of block-assembly plan) Adds a reusable contact-scene fixture (8 scenes spanning 390 to 512k collisions, each padded with interior vertices so to_full_dof performs a genuine surface-to-volume scatter) and Catch2 benchmarks that isolate the three costs of contact Hessian/gradient assembly: 1. per-collision (local) derivative evaluation, 2. global assembly (triplets + setFromTriplets), 3. the reduced-DOF map (CollisionMesh::to_full_dof). Baseline findings: local derivative evaluation is only 1.8-7.6% of Hessian cost; the rest is assembly bookkeeping (42-62%) and to_full_dof SpGEMMs (30-56%). On the largest scene (puffer-ball, 512k collisions) bookkeeping costs ~560 ms per Newton iteration vs 21 ms of derivative math. Also adds a memory-guarded scene probe ([assembly-probe], hidden) that counts broad-phase candidates before building the collision set, since an oversized dhat can exhaust host memory. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Assemble potential derivatives directly in full-mesh DOF (Phase 1) Adds an in_full_dof parameter to Potential<T>::gradient/hessian. When the mesh's DOF map is a pure selection matrix (the default; tracked by the new CollisionMesh::is_selection_dof_map()), stencil vertex IDs are remapped to full-mesh IDs during triplet generation, producing the full-DOF result directly instead of applying to_full_dof afterwards. This eliminates the two serial SpGEMMs (S^T H S), which were 30-56% of end-to-end Hessian cost. With a user-provided displacement map, in_full_dof falls back to to_full_dof internally, so the flag is always safe. Measured end-to-end Hessian speedups: 1.29-1.83x across 8 scenes (390-512k collisions). Gradient folding is not beneficial on large scenes (the thread-local accumulators grow to full_ndof while the SpMV saved is cheap) and is left off by default; documented in the benchmark. Note: the defensive storage-empty path in hessian() now returns a correctly-sized (ndof x ndof) empty matrix instead of 0x0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Introduce HessianAssembler seam for pluggable assembly backends (Phase 2) Extracts the global-matrix construction out of Potential<T>::hessian into an abstract HessianAssembler interface (begin / thread-safe add_local_hessian / end). The historical triplet + setFromTriplets path moves verbatim into TripletHessianAssembler, and hessian() becomes a thin wrapper over the new public Potential<T>::assemble_hessian driver, which also owns the Phase 1 full-DOF stencil remap so every future backend gets it for free. No behavior change; benchmarks confirm collision-DOF assembly times are within run-to-run noise of the previous implementation on all 8 scenes. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add MeshFEMSparse block-CSC assembly backend (Phase 3) Adds MeshFEMHessianAssembler, a HessianAssembler backed by MeshFEMSparse's block-CSC data structures (Mohammadian et al., SIGGRAPH 2026): begin() builds a block sparsity pattern from the collision stencils and add_local_hessian() scatters each local Hessian directly into the value array via MeshFEM's sorted column-merge with per-column spin locks — no triplets, no setFromTriplets. Guarded by IPC_TOOLKIT_WITH_MESHFEM_SPARSE (default OFF). The HessianAssembler seam gains a StencilGetter argument to begin() so pattern-based backends can see stencils up front. The dependency is fetched with CPM DOWNLOAD_ONLY (pinned SHAs + SHA256 archive hashes) and compiled into a minimal static target (matrix data structures and assembly only, no sparse direct solvers), avoiding upstream's PUBLIC -fvisibility=hidden, its solver sources (which clash with Eigen 5's BLAS declarations), and its transitive dependency fetching. Compatibility notes: - MeshFEM targets Eigen 3.4; Eigen 5 removed internal::make_coherent, which MeshFEMCore/AutomaticDifferentiation.hh references (included by SparseMatrices.hh at the root of the header chain). A force-included shim (meshfem_eigen_compat.hpp) reimplements the Eigen 3.4 semantics. - BlockCSCHessian::toEigen/toScalar read out of bounds on empty block columns (impossible for FE Hessians, ubiquitous for contact Hessians: most vertices are collision-free), causing intermittent segfaults. Replaced with a custom direct block-CSC -> symmetric Eigen conversion, which is also ~2x faster than upstream's two-step expansion. Measured on 8 scenes (390-512k collisions), full-DOF Hessian, pattern rebuilt every call: 2.5-11x end-to-end vs the triplet path to an Eigen matrix, 3-15x to the block-CSC format. Matches the triplet assembler to <= 1e-13 relative across scenes x PSD projection x DOF space. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Reuse the block sparsity pattern across assemblies (Phase 4) MeshFEMHessianAssembler is now designed to live across assemblies (e.g., one instance per Newton solve). begin() compares the stencils against the cached block pattern via MeshFEM's detectChangedEntries and reuses it (values-only reset + scatter) unless the contact set gained a new vertex pair or lost more than stale_block_tolerance() blocks; stale blocks assemble to explicit zeros. The Eigen conversion structure (symmetrized pattern + index arrays) is cached the same way, so get_matrix() — now returning a const reference valid until the next begin() — reduces to a parallel value refill while the pattern holds. For callers that know the collision set is identical to the previous assembly (change detection costs a sizable fraction of a rebuild on large scenes), set_assume_unchanged_stencils(true) skips detection entirely; a differing stencil count falls back to detection automatically and debug builds verify the assumption. Amortization is automatic through the existing assemble_hessian seam — no API changes beyond the new accessors. Steady-state contact Hessians (Eigen output included) reach 3.3-29x over the triplet baseline across the 8 benchmark scenes (e.g., cloth-ball 14 ms -> 0.48 ms, puffer-ball ~600 ms -> 32 ms), with reuse semantics covered by new tests (identical/shrunken/grown sets, tolerance behavior, assume-unchanged fallback). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Hybrid gather/scatter gradient assembly (Phase 5) Potential<T>::gradient now selects between two assembly strategies based on problem shape (no API change, no new dependency): - gather (new): local gradients are written to a flat per-slot buffer, a vertex->slot adjacency is built with a parallel counting sort, and each vertex sums its contributions independently. Cost scales with the number of contributions rather than ndof. - scatter+reduce (previous behavior): thread-local dense accumulators whose zero+combine cost scales with ndof. Gather is selected when out_ndof > 4 * num_collisions, the empirical crossover on the benchmark scenes: contact-sparse large meshes get gather (cloth-ball 512-612 -> 381 us, n-body 917 -> 695 us), while collision-dense scenes (rod-twist: 1.3M contributions on 120k DOF, where gather's buffer + adjacency traffic measured 1.6x worse) keep the scatter path. This also removes the Phase 1 caveat that in_full_dof gradients could be slower: with gather the accumulators no longer grow with full_ndof (cloth-ball 595 -> 444 us, puffer-ball 13.7 -> 11.3 ms folded). Summation order remains floating-point nondeterministic on both paths; sorting each gather bucket would make that path reproducible if ever needed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Make MeshFEMSparse block assembly the default (Phase 6) IPC_TOOLKIT_WITH_MESHFEM_SPARSE now defaults to ON (auto-disabled for IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=ColMajor, which the block layout does not support), and Potential<T>::hessian() assembles through the block-CSC backend when compiled in, via a new zero-copy MeshFEMHessianAssembler::take_matrix(). The triplet path remains as the fallback when the option is off. Every existing hessian() caller gets the speedup with no code change: cloth-ball 5-6.7 -> 1.5 ms, armadillo-rollers 11-18 -> 2.2 ms, rod-twist 165-212 -> 29.5 ms, puffer-ball 375-1020 -> 48.9 ms (identical results up to floating-point summation order; full 286-test suite passes in both configurations). The dependency is now pinned to fork commits carrying the two fixes submitted upstream (MeshFEM/MeshFEMCore#1 for Eigen 5 support, MeshFEM/MeshFEMSparse#1 for an out-of-bounds read on empty block columns) -- marked TEMPORARY in the recipe; repoint to upstream SHAs once merged. This allowed deleting the force-included make_coherent compatibility shim entirely. Also: document the HessianAssembler classes in the C++ API docs (with IPC_TOOLKIT_WITH_MESHFEM_SPARSE added to Doxygen's PREDEFINED so the guarded class renders) and add MeshFEMSparse to the optional-dependency docs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Suppress clang-tidy identifier-naming on m_H/m_M matrix members Single-capital-letter names for matrices (H = Hessian, M = matrix) are the codebase's mathematical convention; NOLINT the readability-identifier-naming check on them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Expose the block-CSC matrix and bind the assemblers to Python Adds MeshFEMHessianAssembler::block_matrix(), which returns the assembled matrix in MeshFEMSparse's native block-CSC form so a downstream user can feed it to MeshFEM's block SpMV or Cholesky factorizers instead of paying for the Eigen conversion (0.11 vs 0.30 ms on bunny, 42.8 vs 51.9 ms on puffer-ball). MeshFEM::BlockCSCHessianBase is forward declared, so our header still does not pull in MeshFEMSparse's; callers that want the block matrix include <MeshFEMSparse/BlockCSCHessian.hh> themselves and everyone else pays nothing. Binds assemble_hessian, HessianAssembler, TripletHessianAssembler, and MeshFEMHessianAssembler to Python, so Python callers can now hold an assembler across iterations and get pattern reuse (previously they were limited to the cold path inside hessian()). All three classes are py::is_final(): a Python-defined assembler would take the GIL once per collision, which is hundreds of thousands of times per assembly on the larger scenes. Exercising block_matrix() turned up a third instance of the empty-block- column assumption upstream, in visitDiagonalScalarEntries, which made trace() read the preceding column's storage (1951.93 against a dense trace of 447.82) and addDiag()/setDiag() write to the wrong entries. Fixed in the pinned fork commit alongside the other two (MeshFEM/MeshFEMSparse#1); the tests now cover trace() agreement and that addDiag() rejects a pattern with missing diagonal blocks. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Document in_full_dof and assembler reuse in the simulation tutorial The tutorial still showed to_full_dof as the only way to get full-mesh derivatives, and said nothing about holding an assembler across a Newton solve, which is where most of the speedup lives. Adds an in_full_dof example next to the existing to_full_dof one (with a note on the pure-selection requirement and the silent fallback when a displacement map is present), and a section on reusing a MeshFEMHessianAssembler: what the cached pattern covers, when it is rebuilt, block_matrix() for solvers that speak block CSC, and the assume_unchanged_stencils escape hatch and its caveat. Also drops the now-wrong "two fixes" count for the pinned forks; the MeshFEMSparse PR carries two empty-block-column fixes of its own. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Use a more formal register in the new tutorial prose * Update dependencies.rst * Measure the triplet baseline explicitly in the breakdown table hessian() now routes through whichever backend is compiled in, so the table's local%/asm%/full% columns were comparing the MeshFEM path against itself. Time the triplet assembler directly instead, and rename the column to match. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Harden the assembly entry points Short-circuit an empty collision set in hessian(): building a sparsity pattern to produce an all-zero matrix costs O(ndof) for nothing. Validate dim in MeshFEMHessianAssembler::begin() before dividing by it, so an unsupported dimension throws instead of trapping. Make EIGEN_DONT_VECTORIZE PUBLIC on the MeshFEMSparse target: the setting has to travel with the target, since anything including its headers must agree with how its own translation units were compiled. Return by value rather than through the ternary so the returns are implicitly moved, and include <memory> where unique_ptr is used. * Mark MeshFEMSparse as PRIVATE * Add MeshFEM::Sparse alias and use CPM gh: shorthand for recipe URIs Guards against being included twice (via MeshFEM::Sparse or MeshFEMSparse targets) and aliases MeshFEMSparse as MeshFEM::Sparse for consistent namespaced usage. Also switches CPMAddPackage calls to the gh: URI shorthand instead of manual URL/URL_HASH pairs. * Switch back to canonical MeshFEM source after PR fixes merged * Share gradient assembly across potentials, choosing its strategy per call Five call sites each carried their own copy of the parallel scatter that sums per-stencil local gradients into a dense global vector, and all five used a dense per-thread accumulator whose zero-and-reduce cost scales with the DOF count rather than the contact set. At 5k collisions over 1M DOF that reduction was ~94% of the call. assemble_gradient (utils/gradient_assembler.hpp) now owns that work and picks between two strategies: - a serial scatter over buffered locals: num_slots adds plus zeroing the output once, and bitwise deterministic since the order is fixed. Wins when contact is sparse relative to the DOF count, by up to ~75x. - per-thread accumulators reduced in parallel over DOF blocks. Wins in the dense self-contact regime, by up to ~2x. Their costs scale complementarily and the crossover (out_ndof > num_slots) sits inside the range real scenes span, so it is evaluated per call from the problem shape -- no switch, no global state. Callers whose stencils are unbounded (smooth contact) take the thread-local path only, gated at compile time on the local gradient's storage. Converted: Potential<T>::gradient, NormalPotential::gauss_newton_hessian_diagonal, TangentialPotential::force, TangentialPotential::smooth_contact_force, and SmoothContactPotential::gradient. Summation order changes for the latter four, so their results shift at the ULP level. Tests exercise both strategies by sizing scenes on either side of the crossover, and fail if a later edit to those sizes strands one untested. The sweep that calibrated the crossover is retained behind [.][grad-assembly-sweep] so the condition can be re-derived; the one-off dhat probe test it superseded is gone. Full release suite green (290 cases, 4112585 assertions). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Deprecate CollisionMesh::to_full_dof with m_is_selection_dof_map * Add warn_to_full_dof_is_selection_dof_map helper * Harden the assembler accessors and fix assembly doc/CMake papercuts Code-review follow-ups on the assembly work: - TripletHessianAssembler::get_matrix() and MeshFEMHessianAssembler:: get_matrix()/take_matrix() guarded their storage pointers with only an assert, so calling them before the first assembly was a null dereference in release builds (reachable straight from Python). Throw instead, as block_matrix() already did. - Assert that a stencil's vertex ids are pairwise distinct in MeshFEMHessianAssembler: we scatter with ElementHessianContribAssembler <true>, whose sorted column-merge assumes each block variable occurs once. The broad phase never produces a repeated id, so this only pins the precondition for hand-built collision sets. - Drop the hard-coded caller from the to_full_dof deprecation warning: it claimed Potential::gradient/hessian even when the mapped quantity was unrelated to a potential. Name the alternative instead. - Document that is_selection_dof_map()'s DOF-index formula is layout dependent (C++ and Python docstrings). - Stop FORCE-writing IPC_TOOLKIT_WITH_MESHFEM_SPARSE=OFF into the cache when the layout is ColMajor; shadow it with a normal variable so that switching back to RowMajor re-enables the backend. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Key the CI build cache on the runner CPU FindSIMD compiles with -march=native, so every cached object file carries the building runner's ISA. The cache key was runner.os + config with no CPU component, so objects built on one runner model were restored onto another and the test step died with SIGILL across unrelated suites (friction, candidates, CFL, plane-vertex collisions). The run that populated the cache compiled from scratch in 4m58s and passed; the next run restored 313 MB, built in 1m22s, and failed in 29s. (cherry picked from commit 4769e4b) * Handle degenerate stencils in the MeshFEM Hessian assembler The block-merge scatter (ElementHessianContribAssembler<true>) scans each block column once, so it can only place a stencil whose block variables are pairwise distinct. A stencil that repeats a vertex maps several local blocks onto one global block: the merge re-requests a row it has already passed, writing the wrong entries. Such stencils are reachable. The friction fixtures build an edge-vertex collision whose point is an endpoint of the edge (ids {1, 1, 2}), and the assembled block (1,1) came out 39.4 off the triplet result. Dispatch on the stencil: keep the merge for distinct ids (every collision the broad phase produces, since it rejects candidates that share a vertex) and use MeshFEMSparse's ...SupportingStencilDuplicates otherwise, which accumulates the repeated blocks instead. Covered by a unit test that drives both assemblers with a duplicated id and compares against the triplet result. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Eigen 5.0 removed
Eigen::internal::make_coherentfromunsupported/Eigen/src/AutoDiff/AutoDiffScalar.h, whichAutomaticDifferentiation.hhreferences in thepow(AutoDiffScalar, AutoDiffScalar)overload:Because the call is qualified, name lookup happens at template definition time, so any translation unit including this header fails to compile against Eigen ≥ 5 — even if
powis never instantiated. SinceMeshFEMSparse/SparseMatrices.hhincludes this header, this also makes all of MeshFEMSparse unusable with Eigen 5.This PR reimplements
make_coherentwith the Eigen 3.4 semantics (if exactly one of the two derivative vectors is empty, resize it to match the other and zero it) behind a version guard.Version-check note: Eigen 5 moved to semantic versioning —
EIGEN_WORLD_VERSIONremains 3 forever and the new major version lives inEIGEN_MAJOR_VERSION, so the guard testsEIGEN_MAJOR_VERSION >= 5(Eigen 3.4 hasEIGEN_MAJOR_VERSION == 4, so it keeps using Eigen's own implementation there).Validation: with this change (and no other workarounds), MeshFEMCore headers and all of MeshFEMSparse's matrix/assembly code (
BlockCSCHessian.cc,BorderedSparseHessian.cc,SystemAssembler.hh) compile and pass tests against Eigen 5.0.1, exercised through the IPC Toolkit's new MeshFEMSparse-backed contact Hessian assembly (which pins Eigen 5.0.1).One known remaining Eigen-5 issue, out of scope here:
MeshFEMSparse/Solvers/AccelerateFactorizer.ccconflicts with Eigen 5's BLAS prototype declarations on macOS (misc/blas.hvs the Accelerate framework's).🤖 Generated with Claude Code