Parallel kernels, one accounting layer, one reverse index

CHANGELOG.md

@@ -184,6 +184,23 @@ reconstructed from Git history, release commits, and the shipped docs.
   requires every kernel's Arrow columns to equal its declared outputs.
   `projectionStats` and `estimateCsr` are refused there: they are not kernels
   over a projection.
+### Added
+
+- `ExecutionContext::with_concurrency` permits a kernel to use threads, and
+  `ExecutionContext::work_meter` gives one worker a batched share of the work
+  budget. An execution that does not ask for threads behaves exactly as before.
+- `grust-algorithms` runs `degree`, `pagerank`, `bfs`, `multiSourceBfs`, `wcc`
+  and `projectionStats` in parallel when an execution asks for threads, behind
+  the default `parallel` feature. Results do not depend on the worker count.
+
+### Changed
+
+- PageRank computes an unweighted projection by pulling into each target rather
+  than pushing out of each source when concurrency is requested. The
+  distribution is the same and the two paths agree to within a relative 1e-9,
+  but the low digits of a score differ from the sequential kernel's, at any
+  worker count including one. Evidence that pins PageRank scores to more than
+  nine digits should be regenerated rather than compared across this change.
 ## 0.21.0 — Tadpole — 2026-09-18
 
 ### Cypher language

crates/grust-algorithms/Cargo.toml

@@ -12,7 +12,9 @@ categories.workspace = true
 [features]
 default = ["parallel"]
 arrow = ["dep:arrow-array", "dep:grust-arrow"]
-# Run kernels on the caller's rayon pool. Results are identical at any width.
+# Threads for the kernels, bounded by what each execution asks for. Results are
+# identical at any width. Off compiles only the sequential code, which is also
+# the differential oracle the parallel paths are tested against.
 parallel = ["dep:rayon"]
 
 [dependencies]

crates/grust-algorithms/examples/scaling.rs

@@ -0,0 +1,212 @@
+//! Time the parallel kernels against their sequential paths on a real graph.
+//!
+//!     scaling <edge-list> [--workers 0,1,2,4,8,16] [--iterations 10]
+//!             [--undirected] [--weighted] [--repeat 2]
+//!             [--kernels degree,pagerank,bfs,wcc]
+//!
+//! The edge list is one `source target` pair per line, `#` comments ignored,
+//! as SNAP and the strain harness datasets ship them. The same projection is
+//! built once per worker count, because concurrency belongs to the execution
+//! that owns the projection; build time is reported separately from kernel
+//! time, so a kernel's scaling is not diluted by a sequential load.
+//!
+//! Every row prints the result the kernel produced (a component count, a
+//! reachable count, a score checksum) so that a speedup cannot come from doing
+//! less work: the values must match across worker counts.
+use std::collections::HashMap;
+use std::time::Instant;
+
+use grust_algorithms::{
+    ExecutionContext, ExecutionLimits, GraphProjection, Orientation, PageRankOptions,
+    ProjectionEdge, SnapshotIdentity, bfs, degree, pagerank, weakly_connected_components,
+};
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut args = std::env::args().skip(1);
+    let path = args.next().ok_or("usage: scaling <edge-list> [options]")?;
+    // Zero means the sequential kernels: an execution that never asks for
+    // threads, which is what an embedder gets by default.
+    let mut workers = vec![0usize, 1, 2, 4, 8, 16];
+    let mut iterations = 10usize;
+    // A kernel that builds a cached structure on first use (PageRank builds the
+    // reverse topology) pays for it once. Timing a second call separates that
+    // one-off cost from the iteration it amortises over.
+    let mut repeat = 1usize;
+    let mut orientation = Orientation::Outgoing;
+    let mut weighted = false;
+    let mut kernels = vec![
+        "degree".to_string(),
+        "pagerank".to_string(),
+        "bfs".to_string(),
+        "wcc".to_string(),
+    ];
+    while let Some(flag) = args.next() {
+        let mut value = || args.next().ok_or_else(|| format!("{flag} needs a value"));
+        match flag.as_str() {
+            "--workers" => {
+                workers = value()?
+                    .split(',')
+                    .map(|part| part.parse::<usize>())
+                    .collect::<Result<_, _>>()?;
+            }
+            "--iterations" => iterations = value()?.parse()?,
+            "--repeat" => repeat = value()?.parse()?,
+            "--undirected" => orientation = Orientation::Undirected,
+            // Weights are synthesised from the edge position, so a run is
+            // reproducible and every row has a different scale.
+            "--weighted" => weighted = true,
+            "--kernels" => kernels = value()?.split(',').map(String::from).collect(),
+            other => return Err(format!("unknown argument `{other}`").into()),
+        }
+    }
+
+    let started = Instant::now();
+    let text = std::fs::read_to_string(&path)?;
+    // Rows are numbered on first appearance, as the strain harness does, so the
+    // node table is the graph's own order rather than a sorted one.
+    let mut rows: HashMap<&str, usize> = HashMap::new();
+    let mut names: Vec<&str> = Vec::new();
+    let mut pairs: Vec<(usize, usize)> = Vec::new();
+    for line in text.lines() {
+        if line.starts_with('#') || line.trim().is_empty() {
+            continue;
+        }
+        let mut fields = line.split_whitespace();
+        let (Some(source), Some(target)) = (fields.next(), fields.next()) else {
+            continue;
+        };
+        let source = match rows.get(source) {
+            Some(&row) => row,
+            None => {
+                let row = names.len();
+                names.push(source);
+                rows.insert(source, row);
+                row
+            }
+        };
+        let target = match rows.get(target) {
+            Some(&row) => row,
+            None => {
+                let row = names.len();
+                names.push(target);
+                rows.insert(target, row);
+                row
+            }
+        };
+        pairs.push((source, target));
+    }
+    println!(
+        "{path}: {} nodes, {} edges, parsed in {:.1}s",
+        names.len(),
+        pairs.len(),
+        started.elapsed().as_secs_f64()
+    );
+
+    println!(
+        "{:<10} {:>10} {:>10} {:>10} {:>9}  result",
+        "kernel", "workers", "build s", "kernel s", "vs seq"
+    );
+    let mut baseline: HashMap<String, f64> = HashMap::new();
+    for &threads in &workers {
+        let context = ExecutionContext::new(ExecutionLimits {
+            memory_bytes: 18 << 30,
+            work_units: usize::MAX,
+            batch_rows: 8192,
+            deadline: None,
+        })?;
+        let context = if threads == 0 {
+            context
+        } else {
+            context.with_concurrency(threads)?
+        };
+        let build = Instant::now();
+        let graph = GraphProjection::from_topology(
+            SnapshotIdentity::new("scaling".into(), "r1".into(), "bench".into())?,
+            names.iter().map(|name| (*name).into()).collect(),
+            pairs
+                .iter()
+                .enumerate()
+                .map(|(ordinal, &(source, target))| ProjectionEdge {
+                    source,
+                    target,
+                    ordinal,
+                    id: None,
+                })
+                .collect(),
+            weighted.then(|| {
+                pairs
+                    .iter()
+                    .enumerate()
+                    .map(|(ordinal, _)| 1.0 + (ordinal % 9) as f64)
+                    .collect()
+            }),
+            orientation,
+            &context,
+        )?;
+        let build = build.elapsed().as_secs_f64();
+        for kernel in &kernels {
+            for pass in 1..=repeat {
+                let started = Instant::now();
+                let result = match kernel.as_str() {
+                    "degree" => {
+                        let degrees = degree(&graph)?;
+                        format!("arcs {}", degrees.counts().iter().sum::<usize>())
+                    }
+                    "pagerank" => {
+                        let ranks = pagerank(
+                            &graph,
+                            PageRankOptions {
+                                damping: 0.85,
+                                tolerance: 0.0,
+                                max_iterations: iterations,
+                                personalization: None,
+                            },
+                        )?;
+                        format!(
+                            "iterations {} top {:.9}",
+                            ranks.iterations(),
+                            ranks.values().iter().copied().fold(0.0f64, f64::max)
+                        )
+                    }
+                    "bfs" => {
+                        let distances = bfs(&graph, names[0])?;
+                        let reached = distances.values().iter().filter(|d| d.is_finite()).count();
+                        let sum: f64 = distances.values().iter().filter(|d| d.is_finite()).sum();
+                        format!("reached {reached} sum {sum}")
+                    }
+                    "wcc" => {
+                        let components = weakly_connected_components(&graph)?;
+                        let mut labels: Vec<usize> = components.values().to_vec();
+                        labels.sort_unstable();
+                        labels.dedup();
+                        format!("components {}", labels.len())
+                    }
+                    other => return Err(format!("unknown kernel `{other}`").into()),
+                };
+                let seconds = started.elapsed().as_secs_f64();
+                let key = format!("{kernel}/{pass}");
+                let speedup = match baseline.get(&key) {
+                    Some(&first) => format!("{:.2}x", first / seconds),
+                    None => {
+                        baseline.insert(key, seconds);
+                        "-".to_string()
+                    }
+                };
+                let label = if threads == 0 {
+                    "sequential".to_string()
+                } else {
+                    threads.to_string()
+                };
+                let name = if repeat > 1 {
+                    format!("{kernel}#{pass}")
+                } else {
+                    kernel.clone()
+                };
+                println!(
+                    "{name:<10} {label:>10} {build:>10.1} {seconds:>10.3} {speedup:>9}  {result}"
+                );
+            }
+        }
+    }
+    Ok(())
+}

crates/grust-algorithms/src/biconnected.rs

@@ -153,7 +153,7 @@ pub fn biconnectivity(graph: &GraphProjection) -> Result<Biconnectivity> {
                 let arc = *next;
                 *next += 1;
                 meter.tick(1)?;
-                let slot = adjacency.edge_slots.values[arc];
+                let slot = adjacency.edge_slot(arc);
                 let other = adjacency.targets.values[arc];
                 if slot == arrived_by || other == node {
                     continue;

crates/grust-algorithms/src/buffer.rs

@@ -39,6 +39,36 @@ impl<T> Buffer<T> {
         })
     }
 
+    /// A buffer the caller will fill by index, initialized without charging
+    /// work, for kernels that charge the same elements in their own pass.
+    ///
+    /// [`Self::filled`] charges initialization because its caller does not
+    /// otherwise visit the elements. A kernel that writes every element and
+    /// charges as it goes would charge twice, which would change what a work
+    /// budget means for no change in the work performed.
+    pub(crate) fn indexed(count: usize, value: T, context: &ExecutionContext) -> Result<Self>
+    where
+        T: Clone,
+    {
+        let mut buffer = Self::capacity(count, context)?;
+        context.checkpoint()?;
+        buffer.values.resize(count, value);
+        Ok(buffer)
+    }
+
+    /// As [`Self::indexed`], for elements that cannot be cloned, such as the
+    /// atomics a parallel kernel claims entries in.
+    pub(crate) fn indexed_with(
+        count: usize,
+        make: impl FnMut() -> T,
+        context: &ExecutionContext,
+    ) -> Result<Self> {
+        let mut buffer = Self::capacity(count, context)?;
+        context.checkpoint()?;
+        buffer.values.resize_with(count, make);
+        Ok(buffer)
+    }
+
     pub(crate) fn filled(count: usize, value: T, context: &ExecutionContext) -> Result<Self>
     where
         T: Clone,

crates/grust-algorithms/src/degree.rs

@@ -36,51 +36,69 @@ impl Degrees {
 /// execution takes O(V + A) work and sums weights in stable CSR order. Projection
 /// weights must already be finite and nonnegative; a nonfinite sum is a numerical
 /// error. No reverse index, graph copy, or implicit normalization is built.
+/// Each node's count and strength depend only on its own arc range, so workers
+/// write disjoint slices of the output and need no synchronization. Weighted
+/// sums stay in CSR order within a node, which is why the result does not depend
+/// on how the nodes were divided.
 pub fn degree(graph: &GraphProjection) -> Result<Degrees> {
     let context = graph.execution();
     context.checkpoint()?;
     let adjacency = graph.outgoing();
-    let mut counts = Buffer::capacity(graph.node_count(), context)?;
+    let n = graph.node_count();
+    let mut counts = Buffer::indexed(n, 0usize, context)?;
     let mut strengths = adjacency
         .weights
         .as_ref()
-        .map(|_| Buffer::capacity(graph.node_count(), context))
+        .map(|_| Buffer::indexed(n, 0.0f64, context))
         .transpose()?;
-    // Precharge bounded chunks, as Buffer initialization does, instead of
-    // locking the execution context for every scalar. Arc charges span node
-    // boundaries so low-degree graphs also amortize the accounting lock.
-    let mut arcs_charged_until = 0;
-    for first in (0..graph.node_count()).step_by(1024) {
-        let last = first.saturating_add(1024).min(graph.node_count());
-        context.charge_work(last - first)?;
-        for node in first..last {
-            let start = adjacency.offsets.values[node];
-            let end = adjacency.offsets.values[node + 1];
-            counts.values.push(end - start);
-            if let (Some(weights), Some(strengths)) = (&adjacency.weights, &mut strengths) {
-                let mut sum = 0.0;
-                let mut position = start;
-                while position < end {
-                    if position == arcs_charged_until {
-                        arcs_charged_until =
-                            position.saturating_add(1024).min(weights.values.len());
-                        context.charge_work(arcs_charged_until - position)?;
-                    }
-                    let stop = end.min(arcs_charged_until);
-                    for &weight in &weights.values[position..stop] {
+    let arcs = adjacency.targets.values.len();
+    let workers = crate::parallel::workers(context, n.saturating_add(arcs)).unwrap_or(1);
+    let offsets = &adjacency.offsets.values;
+    let weights = adjacency.weights.as_ref().map(|weights| &weights.values);
+
+    // Counts first: O(V) from the CSR offsets alone.
+    crate::parallel::for_each_chunk(
+        context,
+        workers,
+        &mut counts.values,
+        |first, slice, meter| {
+            meter.charge(slice.len())?;
+            for (index, count) in slice.iter_mut().enumerate() {
+                let node = first + index;
+                *count = offsets[node + 1] - offsets[node];
+            }
+            Ok(())
+        },
+    )?;
+
+    if let (Some(weights), Some(strengths)) = (weights, &mut strengths) {
+        crate::parallel::for_each_chunk(
+            context,
+            workers,
+            &mut strengths.values,
+            |first, slice, meter| {
+                for (index, strength) in slice.iter_mut().enumerate() {
+                    let node = first + index;
+                    let (start, end) = (offsets[node], offsets[node + 1]);
+                    // The counts pass already charged this node; charge its arcs
+                    // before summing them, so a high-degree node cannot run
+                    // unmetered and cancellation stays prompt.
+                    meter.charge(end - start)?;
+                    let mut sum = 0.0;
+                    for &weight in &weights[start..end] {
                         sum += weight;
                     }
-                    position = stop;
-                }
-                // Nonnegative finite inputs cannot recover after overflowing.
-                if !sum.is_finite() {
-                    return Err(AlgorithmError::Numerical(
-                        "weighted degree exceeds finite Float64 domain".into(),
-                    ));
+                    // Nonnegative finite inputs cannot recover after overflowing.
+                    if !sum.is_finite() {
+                        return Err(AlgorithmError::Numerical(
+                            "weighted degree exceeds finite Float64 domain".into(),
+                        ));
+                    }
+                    *strength = sum;
                 }
-                strengths.values.push(sum);
-            }
-        }
+                Ok(())
+            },
+        )?;
     }
     context.checkpoint()?;
     Ok(Degrees {