diff --git a/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs b/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs index cee0d6ef8b..5b349abf05 100644 --- a/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs +++ b/ext/crates/algebra/src/module/homomorphism/free_module_homomorphism.rs @@ -59,31 +59,11 @@ where input_degree: i32, input_index: usize, ) { - assert!(input_degree >= self.source.min_degree()); - assert!(input_index < self.source.dimension(input_degree)); - let output_degree = input_degree - self.degree_shift; assert_eq!( - self.target.dimension(output_degree), + self.target.dimension(input_degree - self.degree_shift), result.as_slice().len() ); - let OperationGeneratorPair { - operation_degree, - generator_degree, - operation_index, - generator_index, - } = *self.source.index_to_op_gen(input_degree, input_index); - - if generator_degree >= self.min_degree() { - let output_on_generator = self.output(generator_degree, generator_index); - self.target.act( - result, - coeff, - operation_degree, - operation_index, - generator_degree - self.degree_shift, - output_on_generator.as_slice(), - ); - } + self.apply_to_basis_element_inner(result, coeff, input_degree, input_index); } fn quasi_inverse(&self, degree: i32) -> Option<&QuasiInverse> { @@ -162,6 +142,60 @@ where &self.outputs[generator_degree][generator_index] } + /// A truncating variant of [`ModuleHomomorphism::apply_to_basis_element`] that allows `result` + /// to span only a prefix of the target's degree-`(input_degree - degree_shift)` basis, rather + /// than the whole thing. + /// + /// The caller must guarantee that the image of the basis element is supported within that + /// prefix; otherwise `act` will panic on an out-of-bounds write. This is used by [Nassau's + /// algorithm](crate::module::homomorphism), which stores each differential truncated to the same + /// prefix (valid by minimality, since a generator's differential lands in the radical), so `act` + /// stops before writing past `result`. This lets a bidegree be computed while treating the + /// target module as if the generators added concurrently in the current internal degree do not + /// yet exist. + pub fn apply_to_basis_element_restricted( + &self, + result: FpSliceMut, + coeff: u32, + input_degree: i32, + input_index: usize, + ) { + assert!(result.as_slice().len() <= self.target.dimension(input_degree - self.degree_shift)); + self.apply_to_basis_element_inner(result, coeff, input_degree, input_index); + } + + /// The shared body of [`Self::apply_to_basis_element_restricted`] and + /// [`ModuleHomomorphism::apply_to_basis_element`], with no check on the length of `result` (the + /// two callers check it differently). + fn apply_to_basis_element_inner( + &self, + result: FpSliceMut, + coeff: u32, + input_degree: i32, + input_index: usize, + ) { + assert!(input_degree >= self.source.min_degree()); + assert!(input_index < self.source.dimension(input_degree)); + let OperationGeneratorPair { + operation_degree, + generator_degree, + operation_index, + generator_index, + } = *self.source.index_to_op_gen(input_degree, input_index); + + if generator_degree >= self.min_degree() { + let output_on_generator = self.output(generator_degree, generator_index); + self.target.act( + result, + coeff, + operation_degree, + operation_index, + generator_degree - self.degree_shift, + output_on_generator.as_slice(), + ); + } + } + pub fn differential_density(&self, degree: i32) -> f32 { let outputs = &self.outputs[degree]; if outputs.is_empty() { diff --git a/ext/src/nassau.rs b/ext/src/nassau.rs index 8d1919136a..15a9112bd1 100644 --- a/ext/src/nassau.rs +++ b/ext/src/nassau.rs @@ -34,6 +34,13 @@ use fp::{ vector::{FpSlice, FpSliceMut, FpVector}, }; use itertools::Itertools; +// If `concurrent` is enabled, the `enumerate`/`for_each` used in `restricted_partial_matrix` come +// from `rayon::prelude::IndexedParallelIterator`, loaded via the `maybe_rayon` prelude. If it is +// disabled, `MaybeIndexedParallelIterator` implements `Iterator`, and those methods come from +// `std::iter::Iterator` instead, leaving this import unused — the same single code path either way. +// Mirrors `algebra::module::homomorphism`. +#[allow(unused_imports)] +use maybe_rayon::prelude::*; use once::OnceBiVec; use sseq::coordinates::{Bidegree, BidegreeGenerator}; @@ -117,6 +124,12 @@ impl MilnorSubalgebra { /// Give a list of basis elements in degree `degree` that has signature `signature`. /// + /// Only basis elements coming from generators of degree strictly less than `max_gen_degree` are + /// considered. Passing [`i32::MAX`] recovers the full mask. Restricting the generator degree is + /// used by [`Resolution::compute_through_stem`] to read a module while ignoring the generators + /// of the current internal degree, which may be added concurrently by another thread. Because + /// generators are laid out in increasing degree, this is exactly a prefix of the full mask. + /// /// This requires passing the algebra for borrow checker reasons. fn signature_mask<'a>( &'a self, @@ -124,35 +137,62 @@ impl MilnorSubalgebra { module: &'a FreeModule, degree: i32, signature: &'a [PPartEntry], + max_gen_degree: i32, ) -> impl Iterator + 'a { - module.iter_gen_offsets([degree]).flat_map( - move |GeneratorData { - gen_deg, - start: [offset], - end: _, - }| { - algebra - .ppart_table(degree - gen_deg) - .iter() - .enumerate() - .filter_map(move |(n, op)| { - if self.has_signature(op, signature) { - Some(offset + n) - } else { - None - } - }) - }, - ) + module + .iter_gen_offsets([degree]) + .take_while(move |gen_data| gen_data.gen_deg < max_gen_degree) + .flat_map( + move |GeneratorData { + gen_deg, + start: [offset], + end: _, + }| { + algebra + .ppart_table(degree - gen_deg) + .iter() + .enumerate() + .filter_map(move |(n, op)| { + if self.has_signature(op, signature) { + Some(offset + n) + } else { + None + } + }) + }, + ) + } + + /// The number of basis elements in `degree` coming from generators of `module` of degree + /// strictly less than `max_gen_degree`. This is the dimension of `module` in `degree` when we + /// pretend the generators of degree `>= max_gen_degree` do not exist. Unlike + /// [`Module::dimension`], it only reads generator counts that are frozen once the previous + /// internal degrees have been committed, so it is safe to call while another thread is adding + /// generators of the current internal degree. + fn restricted_dimension( + module: &FreeModule, + degree: i32, + max_gen_degree: i32, + ) -> usize { + module + .iter_gen_offsets([degree]) + .take_while(|gen_data| gen_data.gen_deg < max_gen_degree) + .map(|gen_data| gen_data.end[0]) + .last() + .unwrap_or(0) } /// Get the matrix of a free module homomorphism when restricted to the subquotient given by /// the signature. + /// + /// Only generators of the target of degree strictly less than `target_max_gen_degree` are used + /// (see [`Self::signature_mask`]). fn signature_matrix( &self, hom: &FreeModuleHomomorphism>, degree: i32, signature: &[PPartEntry], + target_max_gen_degree: i32, ) -> Matrix { let p = hom.prime(); let source = hom.source(); @@ -161,19 +201,28 @@ impl MilnorSubalgebra { let target_degree = degree - hom.degree_shift(); let target_mask: Vec = self - .signature_mask(&algebra, &target, degree - hom.degree_shift(), signature) + .signature_mask( + &algebra, + &target, + target_degree, + signature, + target_max_gen_degree, + ) .collect(); let source_mask: Vec = self - .signature_mask(&algebra, &source, degree, signature) + .signature_mask(&algebra, &source, degree, signature, i32::MAX) .collect(); - let mut scratch = FpVector::new(p, target.dimension(target_degree)); + let mut scratch = FpVector::new( + p, + Self::restricted_dimension(&target, target_degree, target_max_gen_degree), + ); let mut result = Matrix::new(p, source_mask.len(), target_mask.len()); for (mut row, &masked_index) in std::iter::zip(result.iter_mut(), &source_mask) { scratch.set_to_zero(); - hom.apply_to_basis_element(scratch.as_slice_mut(), 1, degree, masked_index); + hom.apply_to_basis_element_restricted(scratch.as_slice_mut(), 1, degree, masked_index); row.add_masked(scratch.as_slice(), 1, &target_mask); } @@ -568,6 +617,28 @@ impl> Resolution { Ok(()) } + /// Build the matrix of `hom` on the basis elements `inputs`, using a target that has been + /// truncated to its first `target_dim` basis elements. This is a version of + /// [`ModuleHomomorphism::get_partial_matrix`] that does not read the (possibly concurrently + /// growing) target dimension. + fn restricted_partial_matrix( + hom: &FreeModuleHomomorphism>, + degree: i32, + inputs: &[usize], + target_dim: usize, + ) -> Matrix { + let mut matrix = Matrix::new(hom.prime(), inputs.len(), target_dim); + if target_dim > 0 { + matrix + .maybe_par_iter_mut() + .enumerate() + .for_each(|(i, row)| { + hom.apply_to_basis_element_restricted(row, 1, degree, inputs[i]) + }); + } + matrix + } + #[tracing::instrument(skip(self), fields(%b, %subalgebra, num_new_gens, density))] fn step_resolution_with_subalgebra( &self, @@ -588,21 +659,32 @@ impl> Resolution { let target = &*self.modules[b.s() - 1]; let algebra = target.algebra(); + // We compute this bidegree treating the target `C_{b.s() - 1}` as if it had no generators of + // degree `>= b.t()`, and `C_{b.s() - 2}` as if it had no generators of degree `>= b.t() - 1`. + // By minimality this loses no information (the differentials we care about land in the + // radical, hence in strictly lower-degree generators), and it makes the computation depend + // only on data that is frozen once `(b.s() - 1, b.t() - 1)` and `(b.s(), b.t() - 1)` have + // been committed. This is what lets [`Self::compute_through_stem`] compute `(b.s(), b.t())` + // concurrently with `(b.s() - 1, b.t())`, which is adding those degree-`b.t()` generators. + let target_bound = b.t(); + let next_bound = b.t() - 1; + let zero_sig = subalgebra.zero_signature(); - let target_dim = target.dimension(b.t()); + let target_dim = MilnorSubalgebra::restricted_dimension(target, b.t(), target_bound); let target_mask: Vec = subalgebra - .signature_mask(&algebra, target, b.t(), &zero_sig) + .signature_mask(&algebra, target, b.t(), &zero_sig, target_bound) .collect(); let target_masked_dim = target_mask.len(); let next = &self.modules[b.s() - 2]; next.compute_basis(b.t()); + let next_dim = MilnorSubalgebra::restricted_dimension(next, b.t(), next_bound); let mut f = if let Some(dir) = self.save_dir().write() { let mut f = self .save_file(SaveKind::NassauQi, b - Bidegree::s_t(1, 0)) .create_file(dir.to_owned(), true); - f.write_u64::(next.dimension(b.t()) as u64)?; + f.write_u64::(next_dim as u64)?; f.write_u64::(target_masked_dim as u64)?; subalgebra.to_bytes(&mut f)?; Some(f) @@ -612,13 +694,18 @@ impl> Resolution { let guard = tracing::info_span!("step", signature = ?zero_sig).entered(); let next_mask: Vec = subalgebra - .signature_mask(&algebra, &self.modules[b.s() - 2], b.t(), &zero_sig) + .signature_mask(&algebra, next, b.t(), &zero_sig, next_bound) .collect(); let next_masked_dim = next_mask.len(); let full_matrix = { let _guard = ParallelGuard::new(); - self.differentials[b.s() - 1].get_partial_matrix(b.t(), &target_mask) + Self::restricted_partial_matrix( + &self.differentials[b.s() - 1], + b.t(), + &target_mask, + next_dim, + ) }; let mut masked_matrix = AugmentedMatrix::new(p, target_masked_dim, [next_masked_dim, target_masked_dim]); @@ -639,14 +726,18 @@ impl> Resolution { &masked_matrix, )?; - if let Some(f) = &mut f - && target.max_computed_degree() < b.t() - { + // The quasi-inverse is always computed on the restricted (degree `< b.t()`) target basis, so + // from the point of view of a later `apply_quasi_inverse` it was computed with "incomplete + // information": the differentials on the degree-`b.t()` generators of the target were not + // available. We flag this unconditionally so the lift is corrected using those differentials + // once they are known. + if let Some(f) = &mut f { f.write_u64::(Magic::Fix as u64)?; } // Compute image - let mut n = subalgebra.signature_matrix(&self.differentials[b.s()], b.t(), &zero_sig); + let mut n = + subalgebra.signature_matrix(&self.differentials[b.s()], b.t(), &zero_sig, target_bound); n.row_reduce(); let next_row = n.rows(); @@ -659,7 +750,7 @@ impl> Resolution { self.add_generators(b, num_new_gens); let mut xs = vec![FpVector::new(p, target_dim); num_new_gens]; - let mut dxs = vec![FpVector::new(p, next.dimension(b.t())); num_new_gens]; + let mut dxs = vec![FpVector::new(p, next_dim); num_new_gens]; for ((x, x_masked), dx) in xs .iter_mut() @@ -682,13 +773,29 @@ impl> Resolution { let _guard = tracing::info_span!("step", ?signature).entered(); target_mask.clear(); next_mask.clear(); - target_mask.extend(subalgebra.signature_mask(&algebra, target, b.t(), &signature)); - next_mask.extend(subalgebra.signature_mask(&algebra, next, b.t(), &signature)); + target_mask.extend(subalgebra.signature_mask( + &algebra, + target, + b.t(), + &signature, + target_bound, + )); + next_mask.extend(subalgebra.signature_mask( + &algebra, + next, + b.t(), + &signature, + next_bound, + )); let full_matrix = { let _guard = ParallelGuard::new(); - self.differential(b.s() - 1) - .get_partial_matrix(b.t(), &target_mask) + Self::restricted_partial_matrix( + &self.differentials[b.s() - 1], + b.t(), + &target_mask, + next_dim, + ) }; let mut masked_matrix = @@ -925,6 +1032,18 @@ impl> Resolution { } /// This function resolves up till a fixed stem instead of a fixed t. + /// + /// The dependency graph we use is the relaxed one: computing `(s, t)` only requires `(s, t - 1)` + /// and `(s - 1, t - 1)` (for `s >= 2`), rather than `(s - 1, t)` and `(s, t - 1)`. The read-only + /// data `(s, t)` needs — the generators of the relevant modules of degree `< t` — is frozen once + /// those two bidegrees have been committed (`step_resolution_with_subalgebra` ignores the + /// degree-`t` generators of the target). This lets `(s, t)` run concurrently with + /// `(s - 1, t)`, the bidegree that produces those degree-`t` generators, and keeps many + /// `t`-diagonals (`n = t - s` fixed) in flight at once, which is where the parallelism comes + /// from. + /// + /// The rows `s = 0` and `s = 1` are kept on the strict schedule: they are cheap, and `step0` and + /// `step1` read their targets through full matrices, so they wait for `(s - 1, t)`. #[tracing::instrument(skip(self), fields(self = self.name, %max))] pub fn compute_through_stem(&self, max: Bidegree) { let _lock = self.lock.lock(); @@ -932,25 +1051,39 @@ impl> Resolution { self.extend_through_degree(max.s()); self.algebra().compute_basis(max.t()); + let min_degree = self.min_degree(); + let max_s = max.s(); + let max_n = max.n(); + + let in_region = |s: i32, t: i32| -> bool { + (0..=max_s).contains(&s) && t >= min_degree && t - s <= max_n + }; + + // `(s, t)` may be computed once its same-row predecessor `(s, t - 1)` and its diagonal + // predecessor are committed. For `s >= 2` the diagonal predecessor is `(s - 1, t - 1)` (the + // relaxed graph); for `s == 1` it is `(0, t)`; `s == 0` has none. `progress[s]` is the + // largest committed `t` in row `s`, so it doubles as a "predecessor committed" test. + let ready = |s: i32, t: i32, progress: &[i32]| -> bool { + in_region(s, t) + && progress[s as usize] >= t - 1 + && match s { + 0 => true, + // Row 1's diagonal predecessor is (0, t). At the stem edge (t = max_n + 1) that + // bidegree lies outside the computed region, so we treat it as satisfied. + 1 => t > max_n || progress[0] >= t, + _ => progress[(s - 1) as usize] >= t - 1, + } + }; + let tracing_span = tracing::Span::current(); maybe_rayon::in_place_scope(|scope| { let _tracing_guard = tracing_span.enter(); - // This algorithm is not optimal, as we compute (s, t) only after computing (s - 1, t) - // and (s, t - 1). In theory, it suffices to wait for (s, t - 1) and (s - 1, t - 1), - // but having the dimensions of the modules change halfway through the computation is - // annoying to do correctly. It seems more prudent to improve parallelism elsewhere. - - // Things that we have finished computing. - let mut progress: Vec = vec![-1; max.s() as usize + 1]; - // We will kickstart the process by pretending we have computed (0, - 1). So - // we must pretend we have only computed up to (0, - 2); - progress[0] = -2; + let mut progress: Vec = vec![min_degree - 1; max_s as usize + 1]; let (sender, receiver) = mpsc::channel(); - SenderData::send(Bidegree::s_t(0, -1), sender); - let f = |b, sender| { + let spawn_bidegree = |b: Bidegree, sender: mpsc::Sender| { if self.has_computed_bidegree(b) { SenderData::send(b, sender); } else { @@ -967,26 +1100,89 @@ impl> Resolution { } }; - while let Ok(SenderData { b, retry, sender }) = receiver.recv() { - if retry { - f(b, sender); - continue; + // Seed the base of every row. A bidegree `(s, min_degree)` has no in-region + // predecessors, so it is not spawned by the wavefront — except `(1, min_degree)`, whose + // diagonal predecessor `(0, min_degree)` is in region, so we let it be spawned instead. + for s in 0..=max_s { + if s != 1 { + spawn_bidegree(Bidegree::s_t(s, min_degree), sender.clone()); } - assert!(progress[b.s() as usize] == b.t() - 1); - progress[b.s() as usize] = b.t(); - - // How far we are from the last one for this s. - let distance = max.n() - b.n() + 1; + } + drop(sender); + + // Bidegrees whose spawned job was stolen onto a worker already inside a critical section + // (`is_in_parallel` set on that worker) and so bounced back a retry rather than causing + // a priority inversion. Because the check is per-thread, a job is only ever bounced when + // its worker is a blocked guard holder; a job picked up by a free worker just runs. Such + // bounces are therefore rare, but when the pool is momentarily saturated we still must + // avoid re-spawning immediately in a tight loop, so we park bounced bidegrees here. + // + // The scheduler thread never holds a guard, so it cannot itself observe when a worker + // frees; instead, while anything is parked we wait on the channel with a short timeout + // and retry the parked work whenever a completion arrives (a worker likely just freed) + // or the timeout elapses (periodic re-check). Incoming messages are still handled the + // instant they arrive; the timeout only governs how promptly we retry while otherwise + // idle. This cannot deadlock: parked entries keep their senders, so the channel stays + // open, and the timeout guarantees parked work is retried until a free worker takes it. + let mut deferred: Vec<(Bidegree, mpsc::Sender)> = Vec::new(); + // How long to wait for a message before retrying parked bidegrees. Small enough that a + // freed worker is used promptly, large enough that the poll is negligible; it only ticks + // while something is parked. + const RETRY_POLL_INTERVAL: std::time::Duration = std::time::Duration::from_micros(100); + + loop { + let event = if deferred.is_empty() { + // Nothing parked: block until a message arrives or all senders drop. + match receiver.recv() { + Ok(data) => Some(data), + Err(_) => break, + } + } else { + // Something parked: wake periodically to retry it. Parked entries hold senders, + // so the channel cannot be disconnected here. + match receiver.recv_timeout(RETRY_POLL_INTERVAL) { + Ok(data) => Some(data), + Err(mpsc::RecvTimeoutError::Timeout) => None, + Err(mpsc::RecvTimeoutError::Disconnected) => break, + } + }; - if b.s() < max.s() && progress[b.s() as usize + 1] == b.t() - 1 { - f(b + Bidegree::s_t(1, 0), sender.clone()); + if let Some(SenderData { b, retry, sender }) = event { + if retry { + // Park until a worker frees; retried below on a completion or timeout. + deferred.push((b, sender)); + continue; + } + assert!(progress[b.s() as usize] == b.t() - 1); + progress[b.s() as usize] = b.t(); + + // Completing `b` can only make ready its same-row successor `(s, t + 1)` and one + // diagonal successor. `ready` requires *both* predecessors, so of the two + // completions that could spawn a given bidegree, only the later one does. + let same_row = b + Bidegree::s_t(0, 1); + let diagonal = if b.s() == 0 { + Bidegree::s_t(1, b.t()) + } else { + b + Bidegree::s_t(1, 1) + }; + + for cand in [same_row, diagonal] { + if ready(cand.s(), cand.t(), &progress) { + spawn_bidegree(cand, sender.clone()); + } + } } - if distance > 1 && (b.s() == 0 || progress[b.s() as usize - 1] > b.t()) { - // We are computing a normal step - f(b + Bidegree::s_t(0, 1), sender); - } else if distance == 1 && b.s() < max.s() { - SenderData::send(b + Bidegree::s_t(0, 1), sender); + // Retry parked bidegrees — reached after a completion (a worker likely just freed) + // or a timeout (periodic re-check), but not after a retry (which `continue`s above, + // so a bounced job waits out the timeout before being retried). Each re-spawned job + // re-checks its own worker's flag: those on a free worker run, those stolen onto a + // blocked guard holder bounce and are re-parked. This stays cheap because per-thread + // bounces are rare, so `deferred` is normally empty. + if !deferred.is_empty() { + for (b, sender) in std::mem::take(&mut deferred) { + spawn_bidegree(b, sender); + } } } }); @@ -1083,6 +1279,7 @@ impl> ChainComplex for Resolution { source, b.t(), &subalgebra.zero_signature(), + i32::MAX, )); let mut scratch0 = FpVector::new(p, zero_mask_dim); @@ -1111,7 +1308,13 @@ impl> ChainComplex for Resolution { assert_eq!(mask.len(), zero_mask_dim + num_new_gens); let target_zero_mask: Vec = subalgebra - .signature_mask(&algebra, target, b.t(), &subalgebra.zero_signature()) + .signature_mask( + &algebra, + target, + b.t(), + &subalgebra.zero_signature(), + i32::MAX, + ) .collect(); let mut matrix = AugmentedMatrix::<3>::new( p, @@ -1143,7 +1346,13 @@ impl> ChainComplex for Resolution { let signature = subalgebra.signature_from_bytes(&mut f).unwrap(); mask.clear(); - mask.extend(subalgebra.signature_mask(&algebra, source, b.t(), &signature)); + mask.extend(subalgebra.signature_mask( + &algebra, + source, + b.t(), + &signature, + i32::MAX, + )); scratch0.set_scratch_vector_size(mask.len()); } else if col == Magic::Fix as usize { // We need to fix the differential problem @@ -1263,6 +1472,62 @@ mod tests { .assert_eq(&res.graded_dimension_string()); } + /// Cross-check the secondary (d2) computation on a *save-backed* Nassau resolution computed with + /// the relaxed [`Resolution::compute_through_stem`] against the standard resolution. This + /// exercises the quasi-inverse save files, which under the relaxed schedule are always written + /// using the "incomplete information" (`Magic::Fix`) path, since a bidegree is computed while + /// ignoring the same-degree generators of its target. + #[test] + fn test_stem_concurrent_secondary() { + use std::sync::Arc; + + use algebra::pair_algebra::PairAlgebra; + + use crate::{ + chain_complex::FreeChainComplex, resolution::secondary::SecondaryResolution, + secondary::SecondaryLift, utils::construct_standard, + }; + + fn d2_chart(lift: &SecondaryResolution) -> String + where + CC: FreeChainComplex, + CC::Algebra: PairAlgebra, + { + let underlying = lift.underlying(); + let mut out = String::new(); + // Mirror the guarded iteration in `SecondaryResolution::e3_page`. + for b in underlying.iter_stem() { + if b.t() > 0 && underlying.has_computed_bidegree(b + Bidegree::n_s(-1, 2)) { + let matrix = lift.homotopy(b.s() + 2).homotopies.hom_k(b.t()); + if matrix.iter().any(|row| !row.is_empty()) { + out.push_str(&format!("d2 {b}: {matrix:?}\n")); + } + } + } + out + } + + let max = Bidegree::n_s(20, 7); + + let dir = tempfile::TempDir::new().unwrap(); + let nassau = crate::utils::construct_nassau("S_2", Some(dir.path().to_owned())).unwrap(); + nassau.compute_through_stem(max); + let nassau_lift = SecondaryResolution::new(Arc::new(nassau)); + nassau_lift.extend_all(); + + let standard = construct_standard::("S_2", None).unwrap(); + standard.compute_through_stem(max); + let standard_lift = SecondaryResolution::new(Arc::new(standard)); + standard_lift.extend_all(); + + assert_eq!( + d2_chart(&nassau_lift), + d2_chart(&standard_lift), + "secondary d2 chart differs between Nassau (save-backed, relaxed schedule) and \ + standard" + ); + } + #[test] fn test_signature_iterator() { let subalgebra = MilnorSubalgebra::new(vec![2, 1]); diff --git a/ext/src/utils.rs b/ext/src/utils.rs index ea1c046eea..7225fff0fa 100644 --- a/ext/src/utils.rs +++ b/ext/src/utils.rs @@ -595,13 +595,31 @@ pub use logging::{LogWriter, ext_tracing_subscriber, init_logging}; pub(crate) mod parallel { - use std::sync::atomic::{AtomicUsize, Ordering}; - - static PARALLEL_DEPTH: AtomicUsize = AtomicUsize::new(0); + use std::cell::Cell; + + thread_local! { + /// Depth of `par_iter_mut` critical sections currently entered *on this thread*. + /// + /// The priority inversion we guard against is narrow: a `step_resolution` job initiates a + /// `par_iter` on some rayon worker, which then blocks in the join and work-steals to stay + /// busy. If that blocked worker steals another (heavy, itself nested-parallel) resolution + /// step, that step runs on it and stalls the critical section it is blocked on. A stolen + /// job runs on the *same* OS thread as the worker that stole it, so a per-thread depth is + /// exactly the right signal: [`is_in_parallel`] reports whether *this* worker is a blocked + /// guard holder. A job picked up by any other worker — idle, or busy on non-critical work — + /// reads zero and is free to run, which is what lets independent bidegrees resolve + /// concurrently. + /// + /// Deliberately per-thread rather than a global count of active critical sections: a global + /// flag blocks *all* new work whenever *any* thread is in a critical section, which under + /// the relaxed wavefront (many bidegrees in flight) is nearly always, producing a retry + /// storm that pegs every core doing no useful work. + static PARALLEL_DEPTH: Cell = const { Cell::new(0) }; + } - /// RAII guard that increments [`PARALLEL_DEPTH`] on creation and decrements on drop. Used to mark - /// regions where `par_iter_mut` work is active, so that `step_resolution` jobs can detect priority - /// inversion and retry. + /// RAII guard that increments this thread's [`PARALLEL_DEPTH`] on creation and decrements it on + /// drop. Used to mark regions where `par_iter_mut` work is active, so that a `step_resolution` + /// job stolen onto a blocked guard holder can detect the priority inversion and retry. pub(crate) struct ParallelGuard { #[allow(dead_code)] span: tracing::span::EnteredSpan, @@ -609,8 +627,11 @@ pub(crate) mod parallel { impl ParallelGuard { pub(crate) fn new() -> Self { - // We use Release to synchronize with `is_in_parallel` - let counter_start = PARALLEL_DEPTH.fetch_add(1, Ordering::Release); + let counter_start = PARALLEL_DEPTH.with(|d| { + let v = d.get(); + d.set(v + 1); + v + }); Self { span: tracing::info_span!( "parallel_guard", @@ -624,14 +645,57 @@ pub(crate) mod parallel { impl Drop for ParallelGuard { fn drop(&mut self) { - // We use Release to synchronize with `is_in_parallel` - let counter_end = PARALLEL_DEPTH.fetch_sub(1, Ordering::Release); - self.span.record("counter_end", counter_end - 1); + let counter_end = PARALLEL_DEPTH.with(|d| { + let v = d.get() - 1; + d.set(v); + v + }); + self.span.record("counter_end", counter_end); } } + /// Whether the *current* thread is inside a `par_iter_mut` critical section, i.e. whether it is + /// a blocked guard holder onto which stealing a resolution step would cause a priority + /// inversion. See [`PARALLEL_DEPTH`]. pub(crate) fn is_in_parallel() -> bool { - PARALLEL_DEPTH.load(Ordering::Acquire) > 0 + PARALLEL_DEPTH.with(|d| d.get() > 0) + } + + #[cfg(test)] + mod tests { + use std::sync::mpsc; + + use super::{ParallelGuard, is_in_parallel}; + + /// A [`ParallelGuard`] held on one thread must not be visible on another: the whole point of + /// making `PARALLEL_DEPTH` thread-local is that a resolution step stolen onto a free worker + /// reads zero. Guards against a regression to a shared counter. + #[test] + fn parallel_guard_is_thread_local() { + assert!(!is_in_parallel()); + + let (held_tx, held_rx) = mpsc::channel(); + let (release_tx, release_rx) = mpsc::channel(); + + let handle = std::thread::spawn(move || { + let guard = ParallelGuard::new(); + // Visible on the thread that holds it. + assert!(is_in_parallel()); + held_tx.send(()).unwrap(); + // Keep the guard alive until the main thread has checked. + release_rx.recv().unwrap(); + drop(guard); + // Cleared once dropped. + assert!(!is_in_parallel()); + }); + + // Once the other thread holds the guard, it must be invisible here. + held_rx.recv().unwrap(); + assert!(!is_in_parallel()); + release_tx.send(()).unwrap(); + + handle.join().unwrap(); + } } } diff --git a/ext/tests/milnor_vs_nassau.rs b/ext/tests/milnor_vs_nassau.rs index 2e6e18b81c..9ae45b85d4 100644 --- a/ext/tests/milnor_vs_nassau.rs +++ b/ext/tests/milnor_vs_nassau.rs @@ -21,3 +21,28 @@ fn compare(#[case] module_name: &str, #[case] max_degree: i32) { assert_eq!(a.graded_dimension_string(), b.graded_dimension_string()); } + +/// Cross-check the parallel [`compute_through_stem`](ext::nassau::Resolution::compute_through_stem) +/// against the standard resolution over a wide stem range. This exercises the relaxed +/// (column-parallel) dependency graph, in which a whole column of internal degree is computed +/// concurrently. +#[rstest] +#[trace] +#[case("S_2", 40, 30)] +#[case("C2", 40, 30)] +#[case("Joker", 30, 20)] +// `max_n = 2^i - 1` puts a nonzero class (`h_i`, i.e. `Sq^{2^i}`) exactly on the `s = 1` stem edge +// `(1, max_n + 1)`, exercising the row-1 boundary of the wavefront where the diagonal predecessor +// `(0, max_n + 1)` lies outside the computed region. +#[case("S_2", 7, 6)] +#[case("S_2", 15, 8)] +fn compare_stem(#[case] module_name: &str, #[case] n: i32, #[case] s: i32) { + let max = Bidegree::n_s(n, s); + let a = construct_standard::(module_name, None).unwrap(); + let b = construct_nassau(module_name, None).unwrap(); + + a.compute_through_stem(max); + b.compute_through_stem(max); + + assert_eq!(a.graded_dimension_string(), b.graded_dimension_string()); +}