diff --git a/doc/includes/references.h b/doc/includes/references.h
index 4e4f282a..8eb37d4d 100644
--- a/doc/includes/references.h
+++ b/doc/includes/references.h
@@ -272,7 +272,7 @@ struct mjData_ {
   mjtNum* qM;                // total inertia (sparse)                           (nM x 1)
 
   // computed by mj_fwdPosition/mj_factorM
-  mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nC x 1)
+  mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nM x 1)
   mjtNum* qLDiagInv;         // 1/diag(D)                                        (nv x 1)
 
   // computed by mj_collisionTree
@@ -305,7 +305,7 @@ struct mjData_ {
   mjtNum* subtree_angmom;    // angular momentum about subtree com               (nbody x 3)
 
   // computed by mj_Euler or mj_implicit
-  mjtNum* qH;                // L'*D*L factorization of modified M               (nC x 1)
+  mjtNum* qH;                // L'*D*L factorization of modified M               (nM x 1)
   mjtNum* qHDiagInv;         // 1/diag(D) of modified M                          (nv x 1)
 
   // computed by mj_resetData
diff --git a/include/mujoco/mjdata.h b/include/mujoco/mjdata.h
index d69a1d4e..ddb42215 100644
--- a/include/mujoco/mjdata.h
+++ b/include/mujoco/mjdata.h
@@ -300,7 +300,7 @@ struct mjData_ {
   mjtNum* qM;                // total inertia (sparse)                           (nM x 1)
 
   // computed by mj_fwdPosition/mj_factorM
-  mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nC x 1)
+  mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nM x 1)
   mjtNum* qLDiagInv;         // 1/diag(D)                                        (nv x 1)
 
   // computed by mj_collisionTree
@@ -333,7 +333,7 @@ struct mjData_ {
   mjtNum* subtree_angmom;    // angular momentum about subtree com               (nbody x 3)
 
   // computed by mj_Euler or mj_implicit
-  mjtNum* qH;                // L'*D*L factorization of modified M               (nC x 1)
+  mjtNum* qH;                // L'*D*L factorization of modified M               (nM x 1)
   mjtNum* qHDiagInv;         // 1/diag(D) of modified M                          (nv x 1)
 
   // computed by mj_resetData
diff --git a/include/mujoco/mjxmacro.h b/include/mujoco/mjxmacro.h
index 2c21aaea..34a6dee4 100644
--- a/include/mujoco/mjxmacro.h
+++ b/include/mujoco/mjxmacro.h
@@ -650,7 +650,7 @@
     X   ( mjtNum,    actuator_moment,   nJmom,       1           ) \
     X   ( mjtNum,    crb,               nbody,       10          ) \
     X   ( mjtNum,    qM,                nM,          1           ) \
-    X   ( mjtNum,    qLD,               nC,          1           ) \
+    X   ( mjtNum,    qLD,               nM,          1           ) \
     X   ( mjtNum,    qLDiagInv,         nv,          1           ) \
     XMJV( mjtNum,    bvh_aabb_dyn,      nbvhdynamic, 6           ) \
     XMJV( mjtByte,   bvh_active,        nbvh,        1           ) \
@@ -667,7 +667,7 @@
     X   ( mjtNum,    qfrc_passive,      nv,          1           ) \
     X   ( mjtNum,    subtree_linvel,    nbody,       3           ) \
     X   ( mjtNum,    subtree_angmom,    nbody,       3           ) \
-    X   ( mjtNum,    qH,                nC,          1           ) \
+    X   ( mjtNum,    qH,                nM,          1           ) \
     X   ( mjtNum,    qHDiagInv,         nv,          1           ) \
     X   ( int,       B_rownnz,          nbody,       1           ) \
     X   ( int,       B_rowadr,          nbody,       1           ) \
diff --git a/mjx/mujoco/mjx/_src/io.py b/mjx/mujoco/mjx/_src/io.py
index a054af78..408647d0 100644
--- a/mjx/mujoco/mjx/_src/io.py
+++ b/mjx/mujoco/mjx/_src/io.py
@@ -372,7 +372,7 @@ def make_data(
         'efc_aref': (nefc, float),
         'efc_force': (nefc, float),
         '_qM_sparse': (m.nM, float),
-        '_qLD_sparse': (m.nC, float),
+        '_qLD_sparse': (m.nM, float),
         '_qLDiagInv_sparse': (m.nv, float),
     }
 
@@ -511,7 +511,7 @@ def get_data_into(
         value = value[dof_i, dof_j]
       elif field.name == 'qLD' and not support.is_sparse(m):
         # TODO(erikfrey): provide correct qLDs
-        value = np.zeros(m.nC)
+        value = np.zeros(m.nM)
       elif field.name == 'qLDiagInv' and not support.is_sparse(m):
         value = np.ones(m.nv)
 
diff --git a/mjx/mujoco/mjx/_src/smooth_test.py b/mjx/mujoco/mjx/_src/smooth_test.py
index bf9176ba..01ce8abb 100644
--- a/mjx/mujoco/mjx/_src/smooth_test.py
+++ b/mjx/mujoco/mjx/_src/smooth_test.py
@@ -92,8 +92,8 @@ class SmoothTest(absltest.TestCase):
     # factor_m
     dx = jax.jit(mjx.factor_m)(mx, mjx.put_data(m, d))
     qLDLegacy = np.zeros(mx.nM)  # pylint:disable=invalid-name
-    for i in range(m.nC):
-      qLDLegacy[d.mapM2C[i]] = d.qLD[i]
+    for i in range(m.nM):
+      qLDLegacy[d.mapM2M[i]] = d.qLD[i]
     _assert_eq(qLDLegacy, dx.qLD, 'qLD')
     _assert_attr_eq(d, dx, 'qLDiagInv')
     _assert_eq(dx._qLD_sparse, np.zeros(0), '_qLD_sparse')
diff --git a/mjx/mujoco/mjx/_src/types.py b/mjx/mujoco/mjx/_src/types.py
index b8945fff..228f3f3d 100644
--- a/mjx/mujoco/mjx/_src/types.py
+++ b/mjx/mujoco/mjx/_src/types.py
@@ -1352,7 +1352,7 @@ class Data(PyTreeNode):
     efc_aref: reference pseudo-acceleration                     (nefc,)
     efc_force: constraint force in constraint space             (nefc,)
     _qM_sparse: qM in sparse representation                     (nM,)
-    _qLD_sparse: qLD in sparse representation                   (nC,)
+    _qLD_sparse: qLD in sparse representation                   (nM,)
     _qLDiagInv_sparse: qLDiagInv in sparse representation       (nv,)
   """  # fmt: skip
   # constant sizes:
diff --git a/python/mujoco/introspect/structs.py b/python/mujoco/introspect/structs.py
index 9a98b07c..c0f605b6 100644
--- a/python/mujoco/introspect/structs.py
+++ b/python/mujoco/introspect/structs.py
@@ -5356,7 +5356,7 @@ STRUCTS: Mapping[str, StructDecl] = dict([
                      inner_type=ValueType(name='mjtNum'),
                  ),
                  doc="L'*D*L factorization of M (sparse)",
-                 array_extent=('nC',),
+                 array_extent=('nM',),
              ),
              StructFieldDecl(
                  name='qLDiagInv',
@@ -5492,7 +5492,7 @@ STRUCTS: Mapping[str, StructDecl] = dict([
                      inner_type=ValueType(name='mjtNum'),
                  ),
                  doc="L'*D*L factorization of modified M",
-                 array_extent=('nC',),
+                 array_extent=('nM',),
              ),
              StructFieldDecl(
                  name='qHDiagInv',
diff --git a/src/engine/engine_core_constraint.c b/src/engine/engine_core_constraint.c
index af773edd..47918d48 100644
--- a/src/engine/engine_core_constraint.c
+++ b/src/engine/engine_core_constraint.c
@@ -2131,7 +2131,7 @@ void mj_projectConstraint(const mjModel* m, mjData* d) {
   // inverse square root of D from inertia LDL decomposition
   mjtNum* sqrtInvD = mjSTACKALLOC(d, nv, mjtNum);
   for (int i=0; i < nv; i++) {
-    int diag = d->C_rowadr[i] + d->C_rownnz[i] - 1;
+    int diag = d->M_rowadr[i] + d->M_rownnz[i] - 1;
     sqrtInvD[i] = 1 / mju_sqrt(d->qLD[diag]);
   }
 
@@ -2167,11 +2167,11 @@ void mj_projectConstraint(const mjModel* m, mjData* d) {
           continue;
         }
 
-        // traverse row j of C, marking new unique nonzeros
-        int nnzC = d->C_rownnz[j];
-        int adrC = d->C_rowadr[j];
-        for (int k=0; k < nnzC; k++) {
-          int c = d->C_colind[adrC + k];
+        // traverse row j of M, marking new unique nonzeros
+        int nnzM = d->M_rownnz[j];
+        int adrM = d->M_rowadr[j];
+        for (int k=0; k < nnzM; k++) {
+          int c = d->M_colind[adrM + k];
           if (marker[c] != r) {
             marker[c] = r;
             nnz++;
@@ -2251,10 +2251,10 @@ void mj_projectConstraint(const mjModel* m, mjData* d) {
           continue;
         }
         int j = B_colind[i];
-        int adrC = d->C_rowadr[j];
-        mju_addToSclSparseInc(B + adrB, d->qLD + adrC,
+        int adrM = d->M_rowadr[j];
+        mju_addToSclSparseInc(B + adrB, d->qLD + adrM,
                               nnzB, B_colind + adrB,
-                              d->C_rownnz[j]-1, d->C_colind + adrC, -b);
+                              d->M_rownnz[j]-1, d->M_colind + adrM, -b);
       }
 
       // B(r,:) <- sqrt(inv(D)) * B(r,:)
diff --git a/src/engine/engine_core_smooth.c b/src/engine/engine_core_smooth.c
index da7fd2c7..2e8cf8bb 100644
--- a/src/engine/engine_core_smooth.c
+++ b/src/engine/engine_core_smooth.c
@@ -1469,11 +1469,11 @@ void mj_factorI_legacy(const mjModel* m, mjData* d, const mjtNum* M, mjtNum* qLD
 // sparse L'*D*L factorizaton of the inertia matrix M, assumed spd
 void mj_factorM(const mjModel* m, mjData* d) {
   TM_START;
-  int nC = m->nC;
-  for (int i=0; i < nC; i++) {
-    d->qLD[i] = d->qM[d->mapM2C[i]];
+  int nM = m->nM;
+  for (int i=0; i < nM; i++) {
+    d->qLD[i] = d->qM[d->mapM2M[i]];
   }
-  mj_factorI(d->qLD, d->qLDiagInv, m->nv, d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+  mj_factorI(d->qLD, d->qLDiagInv, m->nv, d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
   TM_ADD(mjTIMER_POS_INERTIA);
 }
 
@@ -1715,7 +1715,7 @@ void mj_solveM(const mjModel* m, mjData* d, mjtNum* x, const mjtNum* y, int n) {
     mju_copy(x, y, n*m->nv);
   }
   mj_solveLD(x, d->qLD, d->qLDiagInv, m->nv, n,
-             d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+             d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 }
 
 
@@ -1727,14 +1727,14 @@ void mj_solveM_island(const mjModel* m, const mjData* d, mjtNum* restrict x, int
   const mjtNum* qLDiagInv = d->qLDiagInv;
   if (island < 0) {
     mj_solveLD(x, qLD, qLDiagInv, m->nv, 1,
-               d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+               d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
     return;
   }
 
   // local copies of key variables
-  const int* rownnz = d->C_rownnz;
-  const int* rowadr = d->C_rowadr;
-  const int* colind = d->C_colind;
+  const int* rownnz = d->M_rownnz;
+  const int* rowadr = d->M_rowadr;
+  const int* colind = d->M_colind;
   const int* diagnum = m->dof_simplenum;
 
   // local constants: island specific
@@ -1785,9 +1785,9 @@ void mj_solveM2(const mjModel* m, mjData* d, mjtNum* x, const mjtNum* y,
   int nv = m->nv;
 
   // local copies of key variables
-  const int* rownnz = d->C_rownnz;
-  const int* rowadr = d->C_rowadr;
-  const int* colind = d->C_colind;
+  const int* rownnz = d->M_rownnz;
+  const int* rowadr = d->M_rowadr;
+  const int* colind = d->M_colind;
   const int* diagnum = m->dof_simplenum;
   const mjtNum* qLD = d->qLD;
 
diff --git a/src/engine/engine_forward.c b/src/engine/engine_forward.c
index 2bb06790..a61d48a0 100644
--- a/src/engine/engine_forward.c
+++ b/src/engine/engine_forward.c
@@ -770,7 +770,7 @@ static void mj_advance(const mjModel* m, mjData* d,
 // Euler integrator, semi-implicit in velocity, possibly skipping factorisation
 void mj_EulerSkip(const mjModel* m, mjData* d, int skipfactor) {
   TM_START;
-  int nv = m->nv, nC = m->nC;
+  int nv = m->nv, nM = m->nM;
   mj_markStack(d);
   mjtNum* qfrc = mjSTACKALLOC(d, nv, mjtNum);
   mjtNum* qacc = mjSTACKALLOC(d, nv, mjtNum);
@@ -795,22 +795,22 @@ void mj_EulerSkip(const mjModel* m, mjData* d, int skipfactor) {
   else {
     if (!skipfactor) {
       // qH = M + h*diag(B)
-      for (int i=0; i < nC; i++) {
-        d->qH[i] = d->qM[d->mapM2C[i]];
+      for (int i=0; i < nM; i++) {
+        d->qH[i] = d->qM[d->mapM2M[i]];
       }
       for (int i=0; i < nv; i++) {
-        d->qH[d->C_rowadr[i] + d->C_rownnz[i] - 1] += m->opt.timestep * m->dof_damping[i];
+        d->qH[d->M_rowadr[i] + d->M_rownnz[i] - 1] += m->opt.timestep * m->dof_damping[i];
       }
 
       // factorize in-place
-      mj_factorI(d->qH, d->qHDiagInv, nv, d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+      mj_factorI(d->qH, d->qHDiagInv, nv, d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
     }
 
     // solve
     mju_add(qfrc, d->qfrc_smooth, d->qfrc_constraint, nv);
     mju_copy(qacc, qfrc, m->nv);
     mj_solveLD(qacc, d->qH, d->qHDiagInv, nv, 1,
-               d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+               d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
   }
 
   // advance state and time
@@ -939,7 +939,7 @@ void mj_RungeKutta(const mjModel* m, mjData* d, int N) {
 // fully implicit in velocity, possibly skipping factorization
 void mj_implicitSkip(const mjModel* m, mjData* d, int skipfactor) {
   TM_START;
-  int nv = m->nv, nM = m->nM, nD = m->nD, nC = m->nC;
+  int nv = m->nv, nM = m->nM, nD = m->nD;
 
   mj_markStack(d);
   mjtNum* qfrc = mjSTACKALLOC(d, nv, mjtNum);
@@ -987,18 +987,18 @@ void mj_implicitSkip(const mjModel* m, mjData* d, int skipfactor) {
       mju_addScl(MhB, d->qM, MhB, -m->opt.timestep, nM);
 
       // copy into qH
-      for (int i=0; i < nC; i++) {
-        d->qH[i] = MhB[d->mapM2C[i]];
+      for (int i=0; i < nM; i++) {
+        d->qH[i] = MhB[d->mapM2M[i]];
       }
 
       // factorize in-place
-      mj_factorI(d->qH, d->qHDiagInv, nv, d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+      mj_factorI(d->qH, d->qHDiagInv, nv, d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
     }
 
     // solve for qacc: (qM - dt*qDeriv) * qacc = qfrc
     mju_copy(qacc, qfrc, nv);
     mj_solveLD(qacc, d->qH, d->qHDiagInv, nv, 1,
-               d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+               d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
   } else {
     mjERROR("integrator must be implicit or implicitfast");
diff --git a/src/engine/engine_print.c b/src/engine/engine_print.c
index c165ae7a..fcea96b2 100644
--- a/src/engine/engine_print.c
+++ b/src/engine/engine_print.c
@@ -1127,12 +1127,12 @@ void mj_printFormattedData(const mjModel* m, const mjData* d, const char* filena
 
   printInertia("QM", d->qM, m, fp, float_format);
 
-  printSparse("QLD", d->qLD, m->nv, d->C_rownnz,
-              d->C_rowadr, d->C_colind, fp, float_format);
+  printSparse("QLD", d->qLD, m->nv, d->M_rownnz,
+              d->M_rowadr, d->M_colind, fp, float_format);
   printArray("QLDIAGINV", m->nv, 1, d->qLDiagInv, fp, float_format);
 
   if (!mju_isZero(d->qHDiagInv, m->nv)) {
-    printSparse("QH", d->qH, m->nv, d->C_rownnz, d->C_rowadr, d->C_colind, fp, float_format);
+    printSparse("QH", d->qH, m->nv, d->M_rownnz, d->M_rowadr, d->M_colind, fp, float_format);
     printArray("QHDIAGINV", m->nv, 1, d->qHDiagInv, fp, float_format);
   }
 
diff --git a/src/engine/engine_support.c b/src/engine/engine_support.c
index 36cdb6b7..60e05973 100644
--- a/src/engine/engine_support.c
+++ b/src/engine/engine_support.c
@@ -1097,14 +1097,14 @@ void mj_mulM2(const mjModel* m, const mjData* d, mjtNum* res, const mjtNum* vec)
 
     // non-simple: add off-diagonals
     if (!m->dof_simplenum[i]) {
-      int adr = d->C_rowadr[i];
-      res[i] += mju_dotSparse(qLD+adr, vec, d->C_rownnz[i] - 1, d->C_colind+adr, /*flg_unc1=*/0);
+      int adr = d->M_rowadr[i];
+      res[i] += mju_dotSparse(qLD+adr, vec, d->M_rownnz[i] - 1, d->M_colind+adr, /*flg_unc1=*/0);
     }
   }
 
   // res *= sqrt(D)
   for (int i=0; i < nv; i++) {
-    int diag = d->C_rowadr[i] + d->C_rownnz[i] - 1;
+    int diag = d->M_rowadr[i] + d->M_rownnz[i] - 1;
     res[i] *= mju_sqrt(qLD[diag]);
   }
 }
diff --git a/src/engine/engine_vis_interact.c b/src/engine/engine_vis_interact.c
index 3366e731..2e400e56 100644
--- a/src/engine/engine_vis_interact.c
+++ b/src/engine/engine_vis_interact.c
@@ -556,7 +556,7 @@ void mjv_initPerturb(const mjModel* m, mjData* d, const mjvScene* scn, mjvPertur
 
   // compute average spatial inertia at selection point
   for (int i=0; i < nv; i++) {
-    int diag = d->C_rowadr[i] + d->C_rownnz[i] - 1;
+    int diag = d->M_rowadr[i] + d->M_rownnz[i] - 1;
     sqrtInvD[i] = 1 / mju_sqrt(d->qLD[diag]);
   }
   mj_jac(m, d, jac, NULL, selpos, sel);
diff --git a/test/benchmark/factorI_benchmark_test.cc b/test/benchmark/factorI_benchmark_test.cc
index e434df5d..dd73d3da 100644
--- a/test/benchmark/factorI_benchmark_test.cc
+++ b/test/benchmark/factorI_benchmark_test.cc
@@ -44,9 +44,9 @@ static void BM_factorI(benchmark::State& state, bool legacy, bool coil) {
   mj_markStack(d);
 
   // M: mass matrix in CSR format
-  mjtNum* M = mj_stackAllocNum(d, m->nC);
-  for (int i=0; i < m->nC; i++) {
-    M[i] = d->qM[d->mapM2C[i]];
+  mjtNum* M = mj_stackAllocNum(d, m->nM);
+  for (int i=0; i < m->nM; i++) {
+    M[i] = d->qM[d->mapM2M[i]];
   }
 
   // LDlegacy: legacy LD matrix (size nM)
@@ -60,7 +60,7 @@ static void BM_factorI(benchmark::State& state, bool legacy, bool coil) {
       } else {
         mju_copy(d->qLD, M, m->nC);
         mj_factorI(d->qLD, d->qLDiagInv, m->nv,
-                   d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+                   d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
       }
     }
   }
diff --git a/test/benchmark/inertia_benchmark_test.cc b/test/benchmark/inertia_benchmark_test.cc
index 76cd73a4..e787fa24 100644
--- a/test/benchmark/inertia_benchmark_test.cc
+++ b/test/benchmark/inertia_benchmark_test.cc
@@ -46,9 +46,9 @@ static void BM_solve(benchmark::State& state, SolveType type) {
   mj_markStack(d);
 
   // M: mass matrix in CSR format
-  mjtNum* M = mj_stackAllocNum(d, m->nC);
-  for (int i=0; i < m->nC; i++) {
-    M[i] = d->qM[d->mapM2C[i]];
+  mjtNum* M = mj_stackAllocNum(d, m->nM);
+  for (int i=0; i < m->nM; i++) {
+    M[i] = d->qM[d->mapM2M[i]];
   }
 
   // LDlegacy: legacy LD matrix (size nM)
@@ -74,9 +74,9 @@ static void BM_solve(benchmark::State& state, SolveType type) {
         case SolveType::kCsr:
           mju_copy(d->qLD, M, m->nC);
           mj_factorI(d->qLD, d->qLDiagInv, m->nv,
-                     d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+                     d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
           mj_solveLD(res, d->qLD, d->qLDiagInv, m->nv, 1,
-                     d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+                     d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
       }
     }
   }
diff --git a/test/benchmark/solveLD_benchmark_test.cc b/test/benchmark/solveLD_benchmark_test.cc
index 436e39ca..64276ba9 100644
--- a/test/benchmark/solveLD_benchmark_test.cc
+++ b/test/benchmark/solveLD_benchmark_test.cc
@@ -52,9 +52,8 @@ static void BM_solveLD(benchmark::State& state, bool featherstone, bool coil) {
 
   // make legacy matrix
   mjtNum* LDlegacy = mj_stackAllocNum(d, m->nM);
-  mju_zero(LDlegacy, m->nM);
-  for (int i=0; i < m->nC; i++) {
-    LDlegacy[d->mapM2C[i]] = d->qLD[i];
+  for (int i=0; i < m->nM; i++) {
+    LDlegacy[d->mapM2M[i]] = d->qLD[i];
   }
 
   // benchmark
@@ -65,7 +64,7 @@ static void BM_solveLD(benchmark::State& state, bool featherstone, bool coil) {
         mj_solveLD_legacy(m, res, 1, LDlegacy, d->qLDiagInv);
       } else {
         mj_solveLD(res, d->qLD, d->qLDiagInv, m->nv, 1,
-                   d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+                   d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
       }
     }
   }
diff --git a/test/engine/engine_core_smooth_test.cc b/test/engine/engine_core_smooth_test.cc
index d0c5fa00..3263acc5 100644
--- a/test/engine/engine_core_smooth_test.cc
+++ b/test/engine/engine_core_smooth_test.cc
@@ -476,7 +476,7 @@ TEST_F(CoreSmoothTest, FactorI) {
   int nv = model->nv;
   vector<mjtNum> Ldense(nv*nv, 0);
   mju_sparse2dense(Ldense.data(), data->qLD, nv, nv,
-                   data->C_rownnz, data->C_rowadr, data->C_colind);
+                   data->M_rownnz, data->M_rowadr, data->M_colind);
   for (int i=0; i < nv; i++) {
     // set diagonal to 1
     Ldense[i*nv+i] = 1;
@@ -485,7 +485,7 @@ TEST_F(CoreSmoothTest, FactorI) {
   // dense D matrix
   vector<mjtNum> Ddense(nv*nv);
   mju_sparse2dense(Ddense.data(), data->qLD, nv, nv,
-                   data->C_rownnz, data->C_rowadr, data->C_colind);
+                   data->M_rownnz, data->M_rowadr, data->M_colind);
   for (int i=0; i < nv; i++) {
     for (int j=0; j < nv; j++) {
       // zero everything except the diagonal
@@ -521,18 +521,17 @@ TEST_F(CoreSmoothTest, SolveLDs) {
 
   int nv = m->nv;
   int nM = m->nM;
-  int nC = m->nC;
 
   // copy M into LD: Legacy format
-  vector<mjtNum> LDlegacy(nM, 0);
-  for (int i=0; i < nC; i++) {
-    LDlegacy[d->mapM2C[i]] = d->qLD[i];
+  vector<mjtNum> LDlegacy(nM);
+  for (int i=0; i < nM; i++) {
+    LDlegacy[d->mapM2M[i]] = d->qLD[i];
   }
 
   // compare LD and LDs densified matrices
   vector<mjtNum> LDdense(nv*nv);
   mju_sparse2dense(LDdense.data(), d->qLD, nv, nv,
-                   d->C_rownnz, d->C_rowadr, d->C_colind);
+                   d->M_rownnz, d->M_rowadr, d->M_colind);
   vector<mjtNum> LDdense2(nv*nv);
   mj_fullM(m, LDdense2.data(), LDlegacy.data());
 
@@ -551,7 +550,7 @@ TEST_F(CoreSmoothTest, SolveLDs) {
 
   mj_solveLD_legacy(m, vec.data(), 1, LDlegacy.data(), d->qLDiagInv);
   mj_solveLD(vec2.data(), d->qLD, d->qLDiagInv, nv, 1,
-             d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+             d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
   // expect vectors to match up to floating point precision
   for (int i=0; i < nv; i++) {
@@ -573,12 +572,11 @@ TEST_F(CoreSmoothTest, SolveLDmultipleVectors) {
 
   int nv = m->nv;
   int nM = m->nM;
-  int nC = m->nC;
 
   // copy LD into LDlegacy: Legacy format
-  vector<mjtNum> LDlegacy(nM, 0);
-  for (int i=0; i < nC; i++) {
-    LDlegacy[d->mapM2C[i]] = d->qLD[i];
+  vector<mjtNum> LDlegacy(nM);
+  for (int i=0; i < nM; i++) {
+    LDlegacy[d->mapM2M[i]] = d->qLD[i];
   }
 
   // compare n LD and LDs vector solve
@@ -590,7 +588,7 @@ TEST_F(CoreSmoothTest, SolveLDmultipleVectors) {
 
   mj_solveLD_legacy(m, vec.data(), n, LDlegacy.data(), d->qLDiagInv);
   mj_solveLD(vec2.data(), d->qLD, d->qLDiagInv, nv, n,
-             d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+             d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
   // expect vectors to match up to floating point precision
   for (int i=0; i < nv*n; i++) {
@@ -614,7 +612,7 @@ TEST_F(CoreSmoothTest, SolveM2) {
   int nv = m->nv;
   vector<mjtNum> sqrtInvD(nv);
   for (int i=0; i < nv; i++) {
-    int diag = d->C_rowadr[i] + d->C_rownnz[i] - 1;
+    int diag = d->M_rowadr[i] + d->M_rownnz[i] - 1;
     sqrtInvD[i] = 1 / mju_sqrt(d->qLD[diag]);
   }
 
@@ -628,7 +626,7 @@ TEST_F(CoreSmoothTest, SolveM2) {
 
   mj_solveM2(m, d, res.data(), vec.data(), sqrtInvD.data(), n);
   mj_solveLD(vec2.data(), d->qLD, d->qLDiagInv, nv, n,
-             d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+             d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
   // expect equality of dot(v, M^-1 * v) and dot(M^-1/2 * v, M^-1/2 * v)
   for (int i=0; i < n; i++) {
@@ -649,29 +647,29 @@ TEST_F(CoreSmoothTest, FactorIs) {
   mjData* d = mj_makeData(m);
   mj_forward(m, d);
 
-  int nC = m->nC, nM = m->nM, nv = m->nv;
+  int nM = m->nM, nv = m->nv;
 
   // copy qM into into qLDlegacy and factorize
   vector<mjtNum> qLDlegacy(nM);
   mj_factorI_legacy(m, d, d->qM, qLDlegacy.data(), d->qLDiagInv);
 
   // copy qLDlegacy into qLDexpected: CSR format
-  vector<mjtNum> qLDexpected(nC);
-  for (int i=0; i < nC; i++) {
-    qLDexpected[i] = qLDlegacy[d->mapM2C[i]];  // mj_factorIs is in-place
+  vector<mjtNum> qLDexpected(nM);
+  for (int i=0; i < nM; i++) {
+    qLDexpected[i] = qLDlegacy[d->mapM2M[i]];
   }
 
   // copy qM into qLD: CSR format
-  vector<mjtNum> qLD(nC);
-  for (int i=0; i < nC; i++) {
-    qLD[i] = d->qM[d->mapM2C[i]];  // mj_factorIs is in-place
+  vector<mjtNum> qLD(nM);
+  for (int i=0; i < nM; i++) {
+    qLD[i] = d->qM[d->mapM2M[i]];  // mj_factorI is in-place
   }
 
   vector<mjtNum> qLDiagInvExpected(d->qLDiagInv, d->qLDiagInv + nv);
   vector<mjtNum> qLDiagInv(nv, 0);
 
   mj_factorI(qLD.data(), qLDiagInv.data(), nv,
-             d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+             d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
   // expect outputs to match to floating point precision
   EXPECT_THAT(qLD, Pointwise(DoubleNear(1e-12), qLDexpected));
@@ -681,12 +679,12 @@ TEST_F(CoreSmoothTest, FactorIs) {
   vector<mjtNum> LDdense(nv*nv);
 
   mju_sparse2dense(LDdense.data(), qLDexpected.data(), nv, nv,
-                   d->C_rownnz, d->C_rowadr, d->C_colind);
-  PrintMatrix(LDdense.data(), nv, nv, 2);
+                   d->M_rownnz, d->M_rowadr, d->M_colind);
+  PrintMatrix(LDdense.data(), nv, nv, 2, "qLDexpected");
 
-  mju_sparse2dense(LDdense.data(), qLDs.data(), nv, nv,
-                   d->C_rownnz, d->C_rowadr, d->C_colind);
-  PrintMatrix(LDdense.data(), nv, nv, 2);
+  mju_sparse2dense(LDdense.data(), qLD.data(), nv, nv,
+                   d->M_rownnz, d->M_rowadr, d->M_colind);
+  PrintMatrix(LDdense.data(), nv, nv, 2, "qLD");
   */
 
   mj_deleteData(d);
diff --git a/test/engine/engine_derivative_test.cc b/test/engine/engine_derivative_test.cc
index 12b9ca88..027578e7 100644
--- a/test/engine/engine_derivative_test.cc
+++ b/test/engine/engine_derivative_test.cc
@@ -436,7 +436,7 @@ static void LinearSystem(const mjModel* m, mjData* d, mjtNum* A, mjtNum* B) {
       Ac[nv*nv + i*nv + i] = -m->dof_damping[i];
     }
     mj_solveLD(Ac, d->qH, d->qHDiagInv, nv, 2*nv,
-               d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+               d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
 
     // A = [dt*Ac; Ac]
     mju_transpose(A, Ac, 2*nv, nv);
@@ -464,7 +464,7 @@ static void LinearSystem(const mjModel* m, mjData* d, mjtNum* A, mjtNum* B) {
     mju_sparse2dense(Bc, d->actuator_moment, nu, nv, d->moment_rownnz,
                      d->moment_rowadr, d->moment_colind);
     mj_solveLD(Bc, d->qH, d->qHDiagInv, nv, nu,
-               d->C_rownnz, d->C_rowadr, m->dof_simplenum, d->C_colind);
+               d->M_rownnz, d->M_rowadr, m->dof_simplenum, d->M_colind);
     mju_transpose(BcT, Bc, nu, nv);
     mju_scl(B, BcT, dt*dt, nu*nv);
     mju_scl(B+nu*nv, BcT, dt, nu*nv);
diff --git a/test/engine/engine_support_test.cc b/test/engine/engine_support_test.cc
index 8e915328..20dbfa29 100644
--- a/test/engine/engine_support_test.cc
+++ b/test/engine/engine_support_test.cc
@@ -32,7 +32,6 @@ namespace {
 using ::std::vector;
 using ::testing::ContainsRegex;  // NOLINT
 using ::testing::DoubleNear;
-using ::testing::ElementsAreArray;
 using ::testing::Eq;
 using ::testing::MatchesRegex;
 using ::testing::NotNull;
diff --git a/test/fixture.h b/test/fixture.h
index d2860d7c..df568943 100644
--- a/test/fixture.h
+++ b/test/fixture.h
@@ -113,9 +113,10 @@ inline std::vector<mjtNum> AsVector(const mjtNum* array, int n) {
 }
 
 // Prints a matrix to stderr, useful for debugging.
-inline void PrintMatrix(const mjtNum* mat, int nrow, int ncol, int p = 5) {
+inline void PrintMatrix(const mjtNum* mat, int nrow, int ncol, int p = 5,
+                        std::string_view name = "") {
   std::cerr.precision(p);
-  std::cerr << "\n";
+  std::cerr << name << "\n";
   for (int r = 0; r < nrow; r++) {
     for (int c = 0; c < ncol; c++) {
       mjtNum val = mat[c + r*ncol];