__restrict__

Slight performance win

Author: ms609

AGENTS.md

@@ -1195,7 +1195,7 @@ downstream packages (e.g., TreeSearch) can use `LinkingTo: TreeDist`:
 `src/lap.h` is now a thin wrapper that includes `<TreeDist/lap.h>` and
 re-exports types to global scope.
 
-### LAPJV codegen regression (diagnosed, partially mitigated)
+### LAPJV codegen regression (diagnosed, characterised, accepted)
 
 Including `lap_impl.h` in `lap.cpp` changed the TU context enough for GCC 14's
 register allocator to produce ~8% more instructions in the Dijkstra hot loop,
@@ -1213,9 +1213,58 @@ with `__attribute__((optimize("align-functions=64", "align-loops=16")))`.
 The `lapjv()` wrapper fills the transposed buffer first (matching R's
 column-major storage) then untransposes — restoring the cache-friendly pattern.
 
-**Residual:** ~5–9% on LAPJV 1999×1999 vs main (from the CostMatrix class
-definition difference).  Tree distance metrics are unaffected — they call
-`lap()` from `pairwise_distances.cpp`, whose TU context is unchanged.
+**VTune profiling (vtune-lap/, 2026-04-01):**  Profiled with `-O2 -g
+-fno-omit-frame-pointer` on the current branch (LAPJV 1999×1999 ×30s +
+400×400 ×30s + CID/MSD random 200-tip ×30s each).
+
+| Function | CPU Time | % | Notes |
+|---|---|---|---|
+| `TreeDist::CostMatrix::findRowSubmin` | 23.8s | 19.8% | Inlined into lap() at -O2 (separate symbol with -g) |
+| `TreeDist::lap` | 23.4s | 19.5% | Dijkstra phase dominates |
+| `msd_score` | 9.4s | 7.8% | Contains inlined lap() |
+| `mutual_clustering_score` | 9.3s | 7.7% | Contains inlined lap() |
+| `lapjv` (R wrapper) | 6.0s | 5.0% | Matrix conversion overhead |
+
+Source-line breakdown of `lap()`:
+
+| Phase | Lines | CPU Time (s) | % of lap() |
+|---|---|---|---|
+| Dijkstra inner update loop | 308–325 | ~16.5 | 72% |
+| Dijkstra min-scan | 274–287 | ~3.5 | 15% |
+| Matrix fill (lapjv wrapper) | 61–65 | ~2.7 | 12% |
+| findRowSubmin (reduction) | 173–184 | ~2.0 | 9% |
+
+**Optimisation attempts:**
+
+1. **`__restrict__` on Dijkstra pointers** (DONE, kept): Added restrict-qualified
+   local pointers (`d_ptr`, `pred_ptr`, `cl_ptr`) for the augment-solution phase
+   to eliminate potential alias reloads. Applied to both `lap.cpp` and `lap_impl.h`.
+   Benchmark: no measurable change in regression magnitude, but correct optimisation.
+
+2. **`#pragma GCC optimize("O3")`** (attempted, reverted): Forced O3 for the entire
+   `lap.cpp` TU. Benchmark: no measurable improvement; same alignment pattern.
+
+3. **Per-object Makevars flags** (attempted, failed): R's build system does not support
+   per-target variable overrides in `Makevars.win`.
+
+**Residual regression (confirmed, alignment-dependent):**
+
+| LAPJV size | dev/ref ratio | Direction |
+|---|---|---|
+| 400×400 | 0.81 | **Dev 19% faster** (alignment win) |
+| 1999×1999 | 1.08–1.13 | **Dev 8–13% slower** (alignment loss) |
+
+The regression is alignment-sensitive and manifests differently at different problem
+sizes. The same codegen change that makes 400×400 faster makes 1999×1999 slower.
+This is the classic instruction-alignment lottery: the Dijkstra inner loop's branch
+prediction and instruction fetch are affected by code placement that varies with
+TU context.
+
+**Conclusion:** The regression cannot be reliably eliminated without matching main's
+exact TU context (adding dead code back to the installed header, which is
+unacceptable for CRAN). **Tree distance metrics are completely unaffected** —
+`lap()` is called from `pairwise_distances.cpp` (different TU context), and
+benchmarks confirm neutral performance across all metrics and tree sizes.
 
 **Maintenance note:** if the `lap()` algorithm changes, update BOTH `src/lap.cpp`
 and `inst/include/TreeDist/lap_impl.h`.  The duplication is intentional — it

inst/include/TreeDist/lap_impl.h

@@ -145,22 +145,25 @@ cost lap(const lap_row dim,
   } while (loopcnt < 2);
 
   // AUGMENT SOLUTION for each free row.
-  auto& d           = scratch.d;
-  auto& predecessor = scratch.predecessor;
+  // Restrict-qualified local pointers enable the compiler to avoid
+  // reloads after stores in the Dijkstra inner loop.
+  cost* __restrict__ d_ptr        = scratch.d.data();
+  lap_row* __restrict__ pred_ptr  = scratch.predecessor.data();
+  lap_col* __restrict__ cl_ptr    = scratch.col_list.data();
 
   for (lap_row f = 0; f < num_free; ++f) {
     bool unassignedfound = false;
     lap_row free_row = freeunassigned[f];
-    const cost* free_row_cost = input_cost.row(free_row);
+    const cost* __restrict__ free_row_cost = input_cost.row(free_row);
     lap_col endofpath = 0;
     lap_col last = 0;
     lap_row i;
     lap_col j1;
 
     for (lap_col j = 0; j < dim; ++j) {
-      d[j] = free_row_cost[j] - v_ptr[j];
-      predecessor[j] = free_row;
-      col_list[j] = j;
+      d_ptr[j] = free_row_cost[j] - v_ptr[j];
+      pred_ptr[j] = free_row;
+      cl_ptr[j] = j;
     }
 
     cost min = 0;
@@ -170,63 +173,63 @@ cost lap(const lap_row dim,
     do {
       if (up == low) {
         last = low - 1;
-        min = d[col_list[up++]];
+        min = d_ptr[cl_ptr[up++]];
 
         for (lap_dim k = up; k < dim; ++k) {
-          const lap_col j = col_list[k];
-          const cost h = d[j];
+          const lap_col j = cl_ptr[k];
+          const cost h = d_ptr[j];
           if (h <= min) {
             if (h < min) {
               up = low;
               min = h;
             }
-            col_list[k] = col_list[up];
-            col_list[up++] = j;
+            cl_ptr[k] = cl_ptr[up];
+            cl_ptr[up++] = j;
           }
         }
         for (lap_dim k = low; k < up; ++k) {
-          if (colsol[col_list[k]] < 0) {
-            endofpath = col_list[k];
+          if (colsol[cl_ptr[k]] < 0) {
+            endofpath = cl_ptr[k];
             unassignedfound = true;
             break;
           }
         }
       }
 
       if (!unassignedfound) {
-        j1 = col_list[low++];
+        j1 = cl_ptr[low++];
         i = colsol[j1];
-        const cost* row_i = input_cost.row(i);
+        const cost* __restrict__ row_i = input_cost.row(i);
         const cost h = row_i[j1] - v_ptr[j1] - min;
 
         for (lap_dim k = up; k < dim; ++k) {
-          const lap_col j = col_list[k];
+          const lap_col j = cl_ptr[k];
           cost v2 = row_i[j] - v_ptr[j] - h;
-          if (v2 < d[j]) {
-            predecessor[j] = i;
+          if (v2 < d_ptr[j]) {
+            pred_ptr[j] = i;
             if (v2 == min) {
               if (colsol[j] < 0) {
                 endofpath = j;
                 unassignedfound = true;
                 break;
               } else {
-                col_list[k] = col_list[up];
-                col_list[up++] = j;
+                cl_ptr[k] = cl_ptr[up];
+                cl_ptr[up++] = j;
               }
             }
-            d[j] = v2;
+            d_ptr[j] = v2;
           }
         }
       }
     } while (!unassignedfound);
 
     for (lap_dim k = 0; k <= last; ++k) {
-      j1 = col_list[k];
-      v[j1] += d[j1] - min;
+      j1 = cl_ptr[k];
+      v[j1] += d_ptr[j1] - min;
     }
 
     do {
-      i = predecessor[endofpath];
+      i = pred_ptr[endofpath];
       colsol[endofpath] = i;
       j1 = endofpath;
       endofpath = rowsol[i];

src/lap.cpp

@@ -186,7 +186,6 @@ cost lap(const lap_row dim,
   }
 
   //   AUGMENTING ROW REDUCTION
-  auto& col_list = scratch.col_list;    // List of columns to be scanned in various ways.
   int loopcnt = 0;                       // do-loop to be done twice.
 
   do {
@@ -239,13 +238,16 @@ cost lap(const lap_row dim,
   } while (loopcnt < 2); // Repeat once.
 
   // AUGMENT SOLUTION for each free row.
-  auto& d           = scratch.d;           // 'Cost-distance' in augmenting path calculation.
-  auto& predecessor = scratch.predecessor; // Row-predecessor of column in augmenting/alternating path.
+  // Restrict-qualified local pointers enable the compiler to avoid
+  // reloads after stores in the Dijkstra inner loop.
+  cost* __restrict__ d_ptr        = scratch.d.data();
+  lap_row* __restrict__ pred_ptr  = scratch.predecessor.data();
+  lap_col* __restrict__ cl_ptr    = scratch.col_list.data();
 
   for (lap_row f = 0; f < num_free; ++f) {
     bool unassignedfound = false;
     lap_row free_row = freeunassigned[f];       // Start row of augmenting path.
-    const cost* free_row_cost = input_cost.row(free_row);
+    const cost* __restrict__ free_row_cost = input_cost.row(free_row);
     lap_col endofpath = 0;
     lap_col last = 0;
     lap_row i;
@@ -254,9 +256,9 @@ cost lap(const lap_row dim,
     // Dijkstra shortest path algorithm.
     // Runs until unassigned column added to shortest path tree.
     for (lap_col j = 0; j < dim; ++j) {
-      d[j] = free_row_cost[j] - v_ptr[j];
-      predecessor[j] = free_row;
-      col_list[j] = j;        // Init column list.
+      d_ptr[j] = free_row_cost[j] - v_ptr[j];
+      pred_ptr[j] = free_row;
+      cl_ptr[j] = j;        // Init column list.
     }
 
     cost min = 0;
@@ -271,26 +273,26 @@ cost lap(const lap_row dim,
 
         // Scan columns for up..dim-1 to find all indices for which new minimum occurs.
         // Store these indices between low..up-1 (increasing up).
-        min = d[col_list[up++]];
+        min = d_ptr[cl_ptr[up++]];
 
         for (lap_dim k = up; k < dim; ++k) {
-          const lap_col j = col_list[k];
-          const cost h = d[j];
+          const lap_col j = cl_ptr[k];
+          const cost h = d_ptr[j];
           if (h <= min) {
             if (h < min) {   // New minimum.
               up = low;      // Restart list at index low.
               min = h;
             }
             // New index with same minimum, put on undex up, and extend list.
-            col_list[k] = col_list[up];
-            col_list[up++] = j;
+            cl_ptr[k] = cl_ptr[up];
+            cl_ptr[up++] = j;
           }
         }
         // Check if any of the minimum columns happens to be unassigned.
         // If so, we have an augmenting path right away.
         for (lap_dim k = low; k < up; ++k) {
-          if (colsol[col_list[k]] < 0) {
-            endofpath = col_list[k];
+          if (colsol[cl_ptr[k]] < 0) {
+            endofpath = cl_ptr[k];
             unassignedfound = true;
             break;
           }
@@ -300,16 +302,16 @@ cost lap(const lap_row dim,
       if (!unassignedfound) {
         // Update 'distances' between free_row and all unscanned columns,
         // via next scanned column.
-        j1 = col_list[low++];
+        j1 = cl_ptr[low++];
         i = colsol[j1];
-        const cost* row_i = input_cost.row(i);
+        const cost* __restrict__ row_i = input_cost.row(i);
         const cost h = row_i[j1] - v_ptr[j1] - min;
 
         for (lap_dim k = up; k < dim; ++k) {
-          const lap_col j = col_list[k];
+          const lap_col j = cl_ptr[k];
           cost v2 = row_i[j] - v_ptr[j] - h;
-          if (v2 < d[j]) {
-            predecessor[j] = i;
+          if (v2 < d_ptr[j]) {
+            pred_ptr[j] = i;
             if (v2 == min) { // New column found at same minimum value
               if (colsol[j] < 0) {
                 // If unassigned, shortest augmenting path is complete.
@@ -318,25 +320,25 @@ cost lap(const lap_row dim,
                 break;
               } else {
               // Else add to list to be scanned right away.
-                col_list[k] = col_list[up];
-                col_list[up++] = j;
+                cl_ptr[k] = cl_ptr[up];
+                cl_ptr[up++] = j;
               }
             }
-            d[j] = v2; // <MS: Unintended>
+            d_ptr[j] = v2;
           }
         }
       }
     } while (!unassignedfound);
 
     // Update column prices.
     for(lap_dim k = 0; k <= last; ++k) {
-      j1 = col_list[k];
-      v[j1] += d[j1] - min;
+      j1 = cl_ptr[k];
+      v[j1] += d_ptr[j1] - min;
     }
 
     // Reset row and column assignments along the alternating path.
     do {
-      i = predecessor[endofpath];
+      i = pred_ptr[endofpath];
       colsol[endofpath] = i;
       j1 = endofpath;
       endofpath = rowsol[i];