[FEATURE] Add incremental Hessian updates for Newton solver (#1172)

Skip the expensive J^T*D*J Hessian rebuild on subsequent solver
iterations by tracking constraint state changes and incrementally
updating H. Only constraints that changed QUADRATIC state contribute
a delta (±D * J[i] * J[j]) to the stored H, after which H is
re-factorized via the existing Cholesky kernels.

Implementation:
- Fuse state change detection into update_constraint_efc via
  @cache_kernel(track_changes) to avoid a separate kernel launch.
  Uses stream compaction (atomic add) to build a compact list of
  changed constraint IDs per world.
- Add incremental_h_update kernel: scalar per-thread kernel where
  each thread handles one (i,j) element of the lower triangle of H.
  Launches with dim=(nworld, nv*(nv+1)/2) using the standard
  triangular index formula. Skips zero J entries.
- On the first solver iteration per step, use the full JTDAJ path
  (unchanged). On subsequent iterations, use the incremental path:
  update H, then re-factorize and solve.
- Works for both small (nv <= 32, tile cholesky) and blocked
  (nv > 32, blocked cholesky) paths.
- Add correctness test comparing incremental vs full Hessian paths.

Author: adenzler-nvidia