Implement threading for island constraint solving.

Humanoids22 No threads:
 Benchmark Execution Time: 9.475905812s

Humanoids22 with 10 Threads:
 Benchmark Execution Time: 4.871352s

PiperOrigin-RevId: 571214307
Change-Id: I1f4f2c761b4ae6bc8fac1f28c6c695f6c499f339

test/engine/engine_island_test.cc

@@ -14,15 +14,20 @@
 
 // Tests for engine/engine_island.c.
 
+#include <array>
+#include <cstring>
 #include <string>
 #include <vector>
 
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
 #include <mujoco/mjmodel.h>
+#include <mujoco/mjthread.h>
+#include <mujoco/mjxmacro.h>
 #include <mujoco/mujoco.h>
 #include "src/engine/engine_island.h"
 #include "src/engine/engine_util_sparse.h"
+#include "src/thread/thread_pool.h"
 #include "test/fixture.h"
 
 namespace mujoco {
@@ -382,5 +387,61 @@ TEST_F(IslandTest, IslandEfcElliptic) {
   mj_deleteModel(model);
 }
 
+TEST_F(IslandTest, IslandSingleAndMultiThreadedMatch) {
+  auto model_path = GetTestDataFilePath("benchmark/testdata/22_humanoids.xml");
+  std::array<char, 1024> error;
+  mjModel* model =
+      mj_loadXML(model_path.c_str(), nullptr, error.data(), error.size());
+  model->opt.solver = mjSOL_CG;             // use CG solver
+  model->opt.enableflags |= mjENBL_ISLAND;  // enable islands
+
+  mjModel* model_threaded =
+      mj_loadXML(model_path.c_str(), nullptr, error.data(), error.size());
+  model_threaded->opt.solver = mjSOL_CG;             // use CG solver
+  model_threaded->opt.enableflags |= mjENBL_ISLAND;  // enable islands
+
+  mjData* data = mj_makeData(model);
+  mjData* data_threaded = mj_makeData(model_threaded);
+
+  // warm-up rollout to get a steady state
+  static constexpr int kNumWarmupSteps = 100;
+  for (int i = 0; i < kNumWarmupSteps; i++) {
+    mj_step(model, data);
+  }
+
+  // sync the models and data
+  int spec = mjSTATE_INTEGRATION;
+  int size = mj_stateSize(model, spec);
+  std::vector<mjtNum> initial_state(size);
+  mj_getState(model, data, initial_state.data(), spec);
+  mj_setState(model_threaded, data_threaded, initial_state.data(), spec);
+
+  // bind a threadpool to the data_threaded
+  mjThreadPool* threadpool = mju_threadPoolCreate(10);
+  mju_bindThreadPool(data_threaded, threadpool);
+
+  for (int i = 0; i < 10; ++i) {
+    mj_step(model, data);
+    mj_step(model_threaded, data_threaded);
+  }
+
+  // compare the mjData's.
+  {
+    MJDATA_POINTERS_PREAMBLE((model))
+    #define X(type, name, nr, nc)                                     \
+        EXPECT_EQ(std::memcmp(data->name, data_threaded->name,        \
+                              sizeof(type)*(model->nr)*(nc)),         \
+                  0) << "mjData::" #name " differs";
+    MJDATA_POINTERS
+    #undef X
+  }
+
+  mj_deleteData(data);
+  mj_deleteModel(model);
+  mj_deleteData(data_threaded);
+  mj_deleteModel(model_threaded);
+  mju_threadPoolDestroy(threadpool);
+}
+
 }  // namespace
 }  // namespace mujoco