RustedSciThe 0.4.19

#![cfg(test)]

//! Performance and equivalence diagnostics for symbolic BVP code generation.
//!
//! This module compares several execution paths:
//! - the legacy `parallel2` lambdified pipeline,
//! - sequential AOT execution,
//! - parallel AOT execution,
//! - and, for some sweeps, IR-level AOT execution without checked-in Rust fixtures.
//!
//! Baselines and runtime paths used across the file:
//! - legacy BVP residual/Jacobian generation and execution,
//! - IR-level AOT lowering/evaluation,
//! - checked-in compiled fixture modules,
//! - sequential and parallel orchestrator execution over the same chunk plans.
//!
//! The file is intentionally test-heavy. It acts as a lab notebook for several
//! hypotheses about when AOT becomes worthwhile and when intra-call parallelism
//! starts to pay off.
//!
//! # What The Main Test Groups Cover
//!
//! ## 1. Real Small BVP Checks
//! These use `build_real_bvp_damp1_case(...)`.
//! - Problem shape: `2` equations on a modest grid.
//! - Typical unknown count: tens of variables.
//! - Purpose: semantic sanity checks, early perf smoke tests, and regression
//!   tests for the basic codegen/orchestrator path.
//! - Relevant tests include the `real_bvp_*` and small compiled fixture tests.
//!
//! ## 2. Large Single-Problem Chunking Experiments
//! These still use the `2`-equation BVP family, but with larger `n_steps` and
//! multiple emitted fixture variants.
//! - Problem shape: still `2` equations, but roughly hundreds of unknowns.
//! - Purpose: test chunk-size hypotheses such as `O4/O8/O16`, `R4/R8/R16`,
//!   super-block chunking, and orchestration overhead.
//! - Relevant tests include:
//!   - `benchmark_large_real_bvp_*`
//!   - `benchmark_large_real_bvp_compiled_chunk_variants_runtime`
//!   - `benchmark_large_real_bvp_superblock_variants_runtime`
//!
//! ## 3. Stress Multifield BVP Family
//! These use `build_multifield_stress_bvp_case(field_count, n_steps)`.
//! - Problem shape: `field_count` coupled equations over `n_steps` grid points.
//! - Unknown count: approximately `field_count * n_steps`.
//! - Jacobian shape: approximately `(field_count * n_steps) x (field_count * n_steps)`.
//! - Purpose: probe realistic scaling toward production cases with many coupled
//!   right-hand sides.
//! - These are the most important tests for finding where parallel AOT begins
//!   to outperform sequential AOT.
//!
//! ## 4. Compiled Stress Fixtures
//! These are real checked-in Rust modules generated by the AOT pipeline.
//! They are the strongest "production-like" evidence in this file because the
//! measured code is emitted Rust, compiled by Rust, and then executed through
//! the orchestrator layer.
//!
//! Current compiled stress fixtures:
//! - `test_codegen_generated_stress_bvp_fixtures`
//!   - baseline stress case
//!   - currently `field_count = 8`, `n_steps = 96`
//!   - roughly `768` unknowns
//! - `test_codegen_generated_stress_large_bvp_fixtures`
//!   - near-threshold compiled case
//!   - currently `field_count = 16`, `n_steps = 96`
//!   - roughly `1536` unknowns
//! - `test_codegen_generated_stress_above_threshold_bvp_fixtures`
//!   - above-threshold compiled case
//!   - currently `field_count = 32`, `n_steps = 96`
//!   - roughly `3072` unknowns
//!
//! In current measurements:
//! - the `16 x 96` compiled case behaves like a near-threshold case,
//! - the `32 x 96` compiled case behaves like an above-threshold case,
//! - and the IR-level crossover search suggests the practical boundary is
//!   around `field_count = 8` for residuals and around `field_count = 16` for
//!   sparse Jacobian values when `n_steps` is around `96..128`.
//!
//! # Important Named Benchmarks
//!
//! - `benchmark_stress_multifield_bvp_pipeline`
//!   Measures symbolic build, lambdify compilation, AOT planning, AOT assembly,
//!   and emitted-source cost on a stress family.
//!
//! - `benchmark_stress_multifield_scaling_sweep`
//!   IR-level scaling benchmark. This is not a compiled fixture benchmark.
//!   It lowers runtime plans into `LinearBlock`s and compares sequential AOT IR
//!   vs parallel AOT IR over a size sweep.
//!
//! - `benchmark_stress_multifield_parallel_crossover`
//!   Searches for the first stress size where parallel AOT IR becomes faster
//!   than sequential AOT IR for:
//!   - residual evaluation
//!   - sparse Jacobian value evaluation
//!
//! - `benchmark_compiled_stress_multifield_aot_runtime`
//!   Compiled Rust benchmark for the baseline stress fixture.
//!
//! - `benchmark_compiled_large_stress_multifield_aot_runtime`
//!   Compiled Rust benchmark for the near-threshold `16 x 96` stress fixture.
//!
//! - `benchmark_compiled_above_threshold_stress_multifield_aot_runtime`
//!   Compiled Rust benchmark for the above-threshold `32 x 96` stress fixture.
//!
//! # Reading The Results
//!
//! These tests are diagnostics, not strict performance assertions:
//! - they print stage timings and hot-path timings,
//! - they retain correctness checks where practical,
//! - but they do not fail based on performance ratios.
//!
//! The main questions they help answer are:
//! - Does AOT match `parallel2` semantically?
//! - Does sequential AOT beat the legacy lambdified path?
//! - At what problem size does parallel AOT start beating sequential AOT?

use crate::numerical::BVP_Damp::BVP_traits::{MatrixType, VectorType};
use crate::symbolic::codegen::CodegenIR::{CodegenModule, LinearBlock};
use crate::symbolic::codegen::codegen_adapters::{
    banded_ir_blocks, banded_module, residual_and_sparse_module, residual_ir_blocks,
    residual_module, sparse_ir_blocks, sparse_module,
};
use crate::symbolic::codegen::codegen_orchestrator::{
    GeneratedChunkFn, ParallelExecutorConfig, ParallelFallbackPolicy, ParallelResidualExecutor,
    ParallelSparseJacobianExecutor, ResidualChunkBinding, SequentialResidualExecutor,
    SequentialSparseJacobianExecutor, SparseJacobianChunkBinding,
};
use crate::symbolic::codegen::codegen_runtime_api::{
    ResidualChunkingStrategy, RuntimeArguments, recommended_residual_chunking_for_parallelism,
    recommended_row_chunking_for_parallelism,
};
use crate::symbolic::codegen::codegen_tasks::{
    ResidualTask, SparseChunkingStrategy, SparseExprEntry, SparseJacobianTask,
};
use crate::symbolic::codegen::tests::codegen_test_support::{
    benchmark_parallel_config, benchmark_sparse_parallel_config_with_jobs,
    build_combustion_bvp_case, build_combustion_bvp_case_with_backend,
    build_multifield_stress_bvp_case, env_usize, median_duration, per_iter_ns,
    stress_env_list_usize, stress_jacobian_iters, stress_residual_iters, stress_samples,
};
use crate::somelinalg::banded::banded_assembly::BandedAssembly;
use crate::somelinalg::banded::block_tridiagonal_lu::BlockTridiagonalLu;
use crate::somelinalg::banded::block_tridiagonal_lu_consistent::BlockTridiagonalLuConsistent;
use crate::somelinalg::banded::linear_solver::{build_linear_solver, build_solver_for_system};
use crate::somelinalg::banded::node_major_layout::NodeMajorLayout;
use crate::somelinalg::banded::solver_policy::{
    FallbackPolicy, LinearSolverConfig, LinearSolverPolicy,
};
use crate::somelinalg::banded::superblock_layout::SuperBlockLayout;
use crate::somelinalg::banded::{LinearSystemRef, solver_traits::DirectLinearSolver};
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_bvp_fixtures::generated_bvp_fixture;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_bvp_huge_fixtures::generated_bvp_huge_fixture_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_bvp_large_chunk_variants::generated_bvp_large_chunk_variant_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_bvp_large_fixtures::generated_bvp_large_fixture_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_bvp_large_superblock_variants::generated_bvp_large_superblock_variant_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_stress_above_threshold_bvp_fixtures::generated_stress_above_threshold_bvp_fixture_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_stress_bvp_fixtures::generated_stress_bvp_fixture_bindings;
use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_stress_large_bvp_fixtures::generated_stress_large_bvp_fixture_bindings;
use crate::symbolic::symbolic_engine::Expr;
use crate::symbolic::symbolic_functions_BVP::{BvpSymbolicAssemblyBackend, Jacobian};
use faer::col::{Col, ColRef};
use faer::linalg::solvers::Solve;
use nalgebra::DMatrix;
use rayon::prelude::*;
use std::collections::HashSet;
use std::collections::HashMap;
use std::hint::black_box;
use std::time::{Duration, Instant};
use tabled::{builder::Builder, settings::Style};

type EvalFn<'a> = Box<dyn Fn(&[f64]) -> f64 + Send + Sync + 'a>;

/// Builds a simple two-equation BVP case for performance testing.
///
/// The system consists of equations `y - z` and `-z^3` with border conditions
/// `z(0) = 1` and `y(1) = 1`. The problem is discretized on a uniform grid
/// with `n_steps` intervals using forward differences.
///
/// # Arguments
/// * `n_steps` - Number of grid steps (intervals) for discretization.
///
/// # Returns
/// A [`Jacobian`] containing the symbolic residual functions and sparse Jacobian
/// entries for the discretized system.
fn build_real_bvp_damp1_case(n_steps: usize) -> Jacobian {
    let eq1 = Expr::parse_expression("y-z");
    let eq2 = Expr::parse_expression("-z^3");
    let eq_system = vec![eq1, eq2];
    let values = vec!["z".to_string(), "y".to_string()];

    let mut border_conditions = HashMap::new();
    border_conditions.insert("z".to_string(), vec![(0usize, 1.0f64)]);
    border_conditions.insert("y".to_string(), vec![(1usize, 1.0f64)]);

    let mut jac = Jacobian::new();
    jac.discretization_system_BVP_par(
        eq_system,
        values,
        "x".to_string(),
        0.0,
        Some(n_steps),
        None,
        None,
        border_conditions,
        None,
        None,
        "forward".to_string(),
    );
    jac.calc_jacobian_parallel_smart_optimized();
    jac
}

#[test]
#[ignore = "explicitly regenerate checked-in BVP codegen fixtures"]
fn regenerate_real_bvp_codegen_fixtures() {
    use std::fmt::Write as _;

    fn append_bindings(
        source: &mut String,
        bindings_module: &str,
        generated_module: &str,
        n_steps: usize,
        arrays: &[(&str, &[String])],
    ) {
        writeln!(source, "\npub mod {bindings_module} {{").unwrap();
        writeln!(source, "    use super::{generated_module};").unwrap();
        writeln!(source, "    pub const N_STEPS: usize = {n_steps};").unwrap();
        for (name, functions) in arrays {
            writeln!(
                source,
                "    pub const {name}: [fn(&[f64], &mut [f64]); {}] = [",
                functions.len()
            )
            .unwrap();
            for function in *functions {
                writeln!(source, "        {generated_module}::{function},").unwrap();
            }
            source.push_str("    ];\n\n");
        }
        source.push_str("}\n");
    }

    fn names(chunks: impl Iterator<Item = String>) -> Vec<String> {
        chunks.collect()
    }

    fn tasks(
        n_steps: usize,
        residual_chunk: usize,
        sparse_rows: usize,
        module_name: &str,
        residual_name: &str,
        sparse_name: &str,
    ) -> (CodegenModule, Vec<String>, Vec<String>) {
        let jac = build_real_bvp_damp1_case(n_steps);
        let variables_owned = jac.variable_string.clone();
        let variables: Vec<&str> = variables_owned.iter().map(String::as_str).collect();
        let entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let entries = borrowed_sparse_entries(&entries_owned);
        let residual_task = ResidualTask {
            fn_name: residual_name,
            residuals: &jac.vector_of_functions,
            variables: &variables,
            params: None,
        };
        let residual = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: residual_chunk,
        });
        let sparse_task = SparseJacobianTask {
            fn_name: sparse_name,
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &entries,
            variables: &variables,
            params: None,
        };
        let sparse = sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount {
            rows_per_chunk: sparse_rows,
        });
        let residual_names = names(residual.chunks.iter().map(|chunk| chunk.fn_name.clone()));
        let sparse_names = names(sparse.chunks.iter().map(|chunk| chunk.fn_name.clone()));
        let mut module = CodegenModule::new(module_name);
        for chunk in &residual.chunks {
            module.push_residual_block_plan(&chunk.plan);
        }
        for chunk in &sparse.chunks {
            module.push_sparse_values_plan(&chunk.plan);
        }
        (module, residual_names, sparse_names)
    }

    fn write_fixture(file: &str, source: String) {
        let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("src/symbolic/codegen/testing_fixtures")
            .join(file);
        std::fs::write(path, source).expect("generated fixture should be writable");
    }

    let (module, _, _) = tasks(
        8,
        8,
        1,
        "generated_bvp_fixture",
        "fixture_bvp_residual",
        "fixture_bvp_sparse_values",
    );
    write_fixture(
        "test_codegen_generated_bvp_fixtures.rs",
        module.emit_source(),
    );

    let (module, residual, sparse) = tasks(
        32,
        8,
        1,
        "generated_bvp_large_fixture",
        "fixture_bvp32_residual",
        "fixture_bvp32_sparse_values",
    );
    let mut source = module.emit_source();
    append_bindings(
        &mut source,
        "generated_bvp_large_fixture_bindings",
        "generated_bvp_large_fixture",
        32,
        &[
            ("RESIDUAL_CHUNKS", &residual),
            ("SPARSE_VALUE_CHUNKS", &sparse),
        ],
    );
    write_fixture("test_codegen_generated_bvp_large_fixtures.rs", source);

    let (module, residual, sparse) = tasks(
        128,
        64,
        64,
        "generated_bvp_huge_fixture",
        "fixture_bvp128_residual",
        "fixture_bvp128_sparse_values",
    );
    let mut source = module.emit_source();
    append_bindings(
        &mut source,
        "generated_bvp_huge_fixture_bindings",
        "generated_bvp_huge_fixture",
        128,
        &[
            ("RESIDUAL_CHUNKS", &residual),
            ("SPARSE_VALUE_CHUNKS", &sparse),
        ],
    );
    write_fixture("test_codegen_generated_bvp_huge_fixtures.rs", source);

    let jac = build_real_bvp_damp1_case(32);
    let variables_owned = jac.variable_string.clone();
    let variables: Vec<&str> = variables_owned.iter().map(String::as_str).collect();
    let entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let entries = borrowed_sparse_entries(&entries_owned);
    let mut module = CodegenModule::new("generated_bvp_large_chunk_variants");
    let mut arrays = Vec::new();
    for outputs in [4usize, 8, 16] {
        let task = ResidualTask {
            fn_name: match outputs {
                4 => "fixture_bvp32_residual_o4",
                8 => "fixture_bvp32_residual_o8",
                _ => "fixture_bvp32_residual_o16",
            },
            residuals: &jac.vector_of_functions,
            variables: &variables,
            params: None,
        };
        let plan = task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: outputs,
        });
        let residual = names(plan.chunks.iter().map(|chunk| chunk.fn_name.clone()));
        for chunk in &plan.chunks {
            module.push_residual_block_plan(&chunk.plan);
        }
        arrays.push((
            match outputs {
                4 => "RESIDUAL_O4",
                8 => "RESIDUAL_O8",
                _ => "RESIDUAL_O16",
            },
            residual,
        ));
    }
    for rows in [4usize, 8, 16] {
        let task = SparseJacobianTask {
            fn_name: match rows {
                4 => "fixture_bvp32_sparse_r4",
                8 => "fixture_bvp32_sparse_r8",
                _ => "fixture_bvp32_sparse_r16",
            },
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &entries,
            variables: &variables,
            params: None,
        };
        let plan = task.runtime_plan(SparseChunkingStrategy::ByRowCount {
            rows_per_chunk: rows,
        });
        let sparse = names(plan.chunks.iter().map(|chunk| chunk.fn_name.clone()));
        for chunk in &plan.chunks {
            module.push_sparse_values_plan(&chunk.plan);
        }
        arrays.push((
            match rows {
                4 => "SPARSE_R4",
                8 => "SPARSE_R8",
                _ => "SPARSE_R16",
            },
            sparse,
        ));
    }
    let mut source = module.emit_source();
    let borrowed_arrays: Vec<(&str, &[String])> = arrays
        .iter()
        .map(|(name, functions)| (*name, functions.as_slice()))
        .collect();
    append_bindings(
        &mut source,
        "generated_bvp_large_chunk_variant_bindings",
        "generated_bvp_large_chunk_variants",
        32,
        &borrowed_arrays,
    );
    write_fixture("test_codegen_generated_bvp_large_chunk_variants.rs", source);

    let mut module = CodegenModule::new("generated_bvp_large_superblock_variants");
    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual_t2",
        residuals: &jac.vector_of_functions,
        variables: &variables,
        params: None,
    };
    let residual_plan = residual_task
        .runtime_plan(ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let residual = names(
        residual_plan
            .chunks
            .iter()
            .map(|chunk| chunk.fn_name.clone()),
    );
    for chunk in &residual_plan.chunks {
        module.push_residual_block_plan(&chunk.plan);
    }
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_t2",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &entries,
        variables: &variables,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let sparse = names(sparse_plan.chunks.iter().map(|chunk| chunk.fn_name.clone()));
    for chunk in &sparse_plan.chunks {
        module.push_sparse_values_plan(&chunk.plan);
    }
    let mut source = module.emit_source();
    append_bindings(
        &mut source,
        "generated_bvp_large_superblock_variant_bindings",
        "generated_bvp_large_superblock_variants",
        32,
        &[("RESIDUAL_T2", &residual), ("SPARSE_T2", &sparse)],
    );
    write_fixture(
        "test_codegen_generated_bvp_large_superblock_variants.rs",
        source,
    );
}

fn borrowed_sparse_entries(owned_entries: &[(usize, usize, Expr)]) -> Vec<SparseExprEntry<'_>> {
    owned_entries
        .iter()
        .map(|(row, col, expr)| SparseExprEntry {
            row: *row,
            col: *col,
            expr,
        })
        .collect()
}

fn compile_sparse_lambdify<'a>(
    input_names: &[&str],
    owned_entries: &'a [(usize, usize, Expr)],
) -> Vec<EvalFn<'a>> {
    owned_entries
        .iter()
        .map(|(_, _, expr)| {
            let fun = expr.lambdify_borrowed_thread_safe(input_names);
            Box::new(move |args: &[f64]| fun(args)) as EvalFn<'a>
        })
        .collect()
}

fn compile_banded_lambdify<'a>(
    prepared: &'a crate::symbolic::symbolic_functions_BVP::BvpPreparedSparseAotProblem,
) -> Vec<EvalFn<'a>> {
    let banded = prepared.as_prepared_banded_problem();
    let input_names = banded.input_names().to_vec();
    banded
        .jacobian_plan
        .chunks
        .iter()
        .flat_map(|chunk| chunk.entries.iter())
        .map(|entry| {
            let fun = entry
                .expr
                .lambdify_borrowed_thread_safe(input_names.as_slice());
            Box::new(move |args: &[f64]| fun(args)) as EvalFn<'a>
        })
        .collect()
}

fn eval_functions(functions: &[EvalFn<'_>], args: &[f64]) -> Vec<f64> {
    functions.iter().map(|fun| fun(args)).collect()
}

fn banded_assembly_to_dense(assembly: &BandedAssembly) -> DMatrix<f64> {
    let mut dense = DMatrix::zeros(assembly.n(), assembly.n());
    for offset in assembly.min_offset()..=assembly.max_offset() {
        if let Some(diag) = assembly.diag(offset) {
            for (pos, value) in diag.iter().enumerate() {
                let (row, col) = assembly
                    .diag_pos_to_ij(offset, pos)
                    .expect("banded diagonal position should map to dense coordinates");
                dense[(row, col)] = *value;
            }
        }
    }
    dense
}

fn max_abs_diff(lhs: &[f64], rhs: &[f64]) -> f64 {
    lhs.iter()
        .zip(rhs.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max)
}

fn sample_rhs(n: usize) -> Vec<f64> {
    (0..n).map(|index| 1.0 + index as f64 * 1.0e-3).collect()
}

fn base_variable_name(name: &str) -> &str {
    name.rsplit_once('_').map_or(name, |(base, _)| base)
}

fn infer_node_major_layout(jac: &Jacobian, n_total: usize) -> NodeMajorLayout {
    let mut seen = HashSet::new();
    let mut vars_per_node = 0usize;
    for name in &jac.variable_string {
        if seen.insert(base_variable_name(name).to_string()) {
            vars_per_node += 1;
        }
    }
    assert!(
        vars_per_node > 0 && n_total % vars_per_node == 0,
        "failed to infer node-major block size from Jacobian variable names"
    );
    NodeMajorLayout::new(n_total / vars_per_node, vars_per_node)
        .expect("node-major layout should be valid")
}

fn preferred_superblock_group(jac: &Jacobian, n_total: usize) -> usize {
    let layout = infer_node_major_layout(jac, n_total);
    let n_nodes = layout.n_nodes();
    let vars_per_node = layout.vars_per_node();

    // The small 2-field damped cases factorize well with the natural
    // node-major block size. For combustion-like 6-field chains we observed
    // that `g=1` is too strict for the block solver, while `g=2` is the best
    // current structured baseline. Keep the policy deliberately simple here so
    // the benchmark answers "does the new solver help?" rather than
    // "can we tune every scenario by hand?".
    if vars_per_node >= 6 && n_nodes % 2 == 0 {
        2
    } else {
        1
    }
}

fn solve_sparse_lu(matrix: &faer::sparse::SparseColMat<usize, f64>, rhs: &[f64]) -> Vec<f64> {
    let lu = matrix
        .sp_lu()
        .expect("sparse LU factorization should succeed");
    let rhs_col = faer::Col::from_fn(rhs.len(), |index| rhs[index]);
    let solution = lu.solve(&rhs_col);
    solution.iter().copied().collect()
}

fn solve_banded_block_tridiagonal(
    jac: &Jacobian,
    assembly: &BandedAssembly,
    rhs: &[f64],
) -> BandedSolveDiagnostics {
    let layout = infer_node_major_layout(jac, assembly.n());
    let block_tridiagonal = assembly
        .to_block_tridiagonal(layout.n_blocks(), layout.block_size())
        .expect("banded assembly should match node-major block-tridiagonal layout");
    let native_attempt = build_linear_solver(
        &block_tridiagonal,
        LinearSolverConfig {
            policy: LinearSolverPolicy::ForceBlockTridiagonal,
            fallback: FallbackPolicy::Never,
            iterative_refinement_steps: 0,
        },
    );
    let (native_block_status, native_block_error) = match native_attempt {
        Ok(_) => ("ok".to_string(), "-".to_string()),
        Err(err) => ("failed".to_string(), format!("{err:?}")),
    };
    let solver = build_linear_solver(
        &block_tridiagonal,
        LinearSolverConfig {
            policy: LinearSolverPolicy::Auto,
            fallback: FallbackPolicy::ToFaerSparse,
            iterative_refinement_steps: 0,
        },
    )
    .expect("banded linear solver factorization should succeed");
    let solver_backend = solver.backend_name().to_string();
    let mut rhs_owned = rhs.to_vec();
    solver
        .solve_in_place(rhs_owned.as_mut_slice())
        .expect("banded linear solve should succeed");
    BandedSolveDiagnostics {
        solution: rhs_owned,
        solver_backend,
        layout,
        native_block_status,
        native_block_error,
    }
}

fn solve_banded_block_tridiagonal_consistent(
    jac: &Jacobian,
    assembly: &BandedAssembly,
    rhs: &[f64],
    refinement_steps: usize,
) -> BandedSolveDiagnostics {
    let node_layout = infer_node_major_layout(jac, assembly.n());
    let nodes_per_superblock = preferred_superblock_group(jac, assembly.n());
    let (layout, block_tridiagonal) = if nodes_per_superblock == 1 {
        (
            node_layout,
            assembly
                .to_block_tridiagonal(node_layout.n_blocks(), node_layout.block_size())
                .expect("banded assembly should match node-major block-tridiagonal layout"),
        )
    } else {
        let super_layout = SuperBlockLayout::new(
            node_layout.n_nodes(),
            node_layout.vars_per_node(),
            nodes_per_superblock,
        )
        .expect("superblock layout should be constructible for consistent solve");
        (
            NodeMajorLayout::new(super_layout.n_blocks(), super_layout.block_size())
                .expect("superblock metadata should map to synthetic node-major layout"),
            assembly
                .finalize_superblock_tridiagonal(&super_layout)
                .expect("banded assembly should match chosen superblock layout"),
        )
    };

    let solver = build_solver_for_system(
        LinearSystemRef::BlockTridiagonal(&block_tridiagonal),
        LinearSolverConfig {
            policy: LinearSolverPolicy::ForceBlockTridiagonalConsistent,
            fallback: FallbackPolicy::Never,
            iterative_refinement_steps: refinement_steps,
        },
    )
    .expect("consistent banded factorization should succeed");
    let mut rhs_owned = rhs.to_vec();
    solver
        .solve_in_place(rhs_owned.as_mut_slice())
        .expect("consistent banded solve should succeed");

    BandedSolveDiagnostics {
        solution: rhs_owned,
        solver_backend: solver.backend_name().to_string(),
        layout,
        native_block_status: "ok".to_string(),
        native_block_error: "-".to_string(),
    }
}

fn eval_ir_blocks_sequential(
    args: &[f64],
    blocks: &[(usize, LinearBlock)],
    output_len: usize,
) -> Vec<f64> {
    let mut out = vec![0.0; output_len];
    for (offset, block) in blocks {
        let len = block.outputs.len();
        block.eval_into(args, &mut out[*offset..(*offset + len)]);
    }
    out
}

fn eval_ir_blocks_parallel(
    args: &[f64],
    blocks: &[(usize, LinearBlock)],
    output_len: usize,
) -> Vec<f64> {
    let mut out = vec![0.0; output_len];
    let pieces: Vec<(usize, Vec<f64>)> = blocks
        .par_iter()
        .map(|(offset, block)| (*offset, block.eval(args)))
        .collect();

    for (offset, values) in pieces {
        out[offset..(offset + values.len())].copy_from_slice(values.as_slice());
    }

    out
}

fn eval_generated_bvp_residual_values(args: &[f64]) -> Vec<f64> {
    let mut out = vec![0.0; 16];
    generated_bvp_fixture::fixture_bvp_residual_chunk_0(args, &mut out[0..8]);
    generated_bvp_fixture::fixture_bvp_residual_chunk_1(args, &mut out[8..16]);
    out
}

fn eval_generated_bvp_residual_col(args: &[f64]) -> Col<f64> {
    let values = eval_generated_bvp_residual_values(args);
    ColRef::from_slice(values.as_slice()).to_owned()
}

fn eval_generated_bvp_sparse_values(args: &[f64]) -> Vec<f64> {
    let mut outputs = Vec::new();
    let chunk_fns: [fn(&[f64], &mut [f64]); 16] = [
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_0,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_1,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_2,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_3,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_4,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_5,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_6,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_7,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_8,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_9,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_10,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_11,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_12,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_13,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_14,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_15,
    ];
    let chunk_sizes = [2usize, 2, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 2];

    for (chunk_fn, chunk_size) in chunk_fns.into_iter().zip(chunk_sizes) {
        let mut chunk_out = vec![0.0; chunk_size];
        chunk_fn(args, &mut chunk_out);
        outputs.extend(chunk_out);
    }

    outputs
}

fn eval_residual_values_from_plan_and_chunks(
    args: &[f64],
    plan: &crate::symbolic::codegen::codegen_runtime_api::ResidualRuntimePlan<'_>,
    chunk_fns: &[GeneratedChunkFn],
) -> Vec<f64> {
    assert_eq!(plan.chunks.len(), chunk_fns.len());
    let mut out = vec![0.0; plan.output_len];
    for (chunk, chunk_fn) in plan.chunks.iter().zip(chunk_fns.iter().copied()) {
        let start = chunk.output_offset;
        let end = start + chunk.residuals.len();
        chunk_fn(args, &mut out[start..end]);
    }
    out
}

fn eval_sparse_values_from_plan_and_chunks(
    args: &[f64],
    plan: &crate::symbolic::codegen::codegen_runtime_api::SparseJacobianRuntimePlan<'_>,
    chunk_fns: &[GeneratedChunkFn],
) -> Vec<f64> {
    assert_eq!(plan.chunks.len(), chunk_fns.len());
    let mut out = vec![0.0; plan.nnz()];
    for (chunk, chunk_fn) in plan.chunks.iter().zip(chunk_fns.iter().copied()) {
        let start = chunk.value_offset;
        let end = start + chunk.entries.len();
        chunk_fn(args, &mut out[start..end]);
    }
    out
}

fn bind_residual_fixture_chunks<'a>(
    plan: &'a crate::symbolic::codegen::codegen_runtime_api::ResidualRuntimePlan<'a>,
    chunk_fns: &[GeneratedChunkFn],
) -> Vec<ResidualChunkBinding<'a>> {
    assert_eq!(plan.chunks.len(), chunk_fns.len());
    plan.chunks
        .iter()
        .zip(chunk_fns.iter().copied())
        .map(|(chunk, eval)| ResidualChunkBinding {
            fn_name: chunk.fn_name.as_str(),
            output_offset: chunk.output_offset,
            output_len: chunk.residuals.len(),
            eval,
        })
        .collect()
}

fn bind_sparse_fixture_chunks<'a>(
    plan: &'a crate::symbolic::codegen::codegen_runtime_api::SparseJacobianRuntimePlan<'a>,
    chunk_fns: &[GeneratedChunkFn],
) -> Vec<SparseJacobianChunkBinding<'a>> {
    assert_eq!(plan.chunks.len(), chunk_fns.len());
    plan.chunks
        .iter()
        .zip(chunk_fns.iter().copied())
        .map(|(chunk, eval)| SparseJacobianChunkBinding {
            fn_name: chunk.fn_name.as_str(),
            value_offset: chunk.value_offset,
            value_len: chunk.entries.len(),
            eval,
        })
        .collect()
}

fn sparse_strategy_label(strategy: SparseChunkingStrategy) -> &'static str {
    match strategy {
        SparseChunkingStrategy::Whole => "Whole",
        SparseChunkingStrategy::ByTargetChunkCount { .. } => "ByTargetChunkCount",
        SparseChunkingStrategy::ByNonZeroCount { .. } => "ByNonZeroCount",
        SparseChunkingStrategy::ByRowCount { .. } => "ByRowCount",
    }
}

fn residual_variant_chunk_fns(outputs_per_chunk: usize) -> &'static [GeneratedChunkFn] {
    match outputs_per_chunk {
        4 => &generated_bvp_large_chunk_variant_bindings::RESIDUAL_O4,
        8 => &generated_bvp_large_chunk_variant_bindings::RESIDUAL_O8,
        16 => &generated_bvp_large_chunk_variant_bindings::RESIDUAL_O16,
        other => panic!("unsupported residual chunk variant: {other}"),
    }
}

fn sparse_variant_chunk_fns(rows_per_chunk: usize) -> &'static [GeneratedChunkFn] {
    match rows_per_chunk {
        4 => &generated_bvp_large_chunk_variant_bindings::SPARSE_R4,
        8 => &generated_bvp_large_chunk_variant_bindings::SPARSE_R8,
        16 => &generated_bvp_large_chunk_variant_bindings::SPARSE_R16,
        other => panic!("unsupported sparse chunk variant: {other}"),
    }
}

fn residual_superblock_chunk_fns(target_chunks: usize) -> &'static [GeneratedChunkFn] {
    match target_chunks {
        2 => &generated_bvp_large_superblock_variant_bindings::RESIDUAL_T2,
        other => panic!("unsupported residual superblock variant: {other}"),
    }
}

fn sparse_superblock_chunk_fns(target_chunks: usize) -> &'static [GeneratedChunkFn] {
    match target_chunks {
        2 => &generated_bvp_large_superblock_variant_bindings::SPARSE_T2,
        other => panic!("unsupported sparse superblock variant: {other}"),
    }
}

fn stress_residual_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_bvp_fixture_bindings::RESIDUAL_CHUNKS
}

fn stress_sparse_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_bvp_fixture_bindings::SPARSE_VALUE_CHUNKS
}

fn stress_large_residual_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_large_bvp_fixture_bindings::RESIDUAL_CHUNKS
}

fn stress_large_sparse_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_large_bvp_fixture_bindings::SPARSE_VALUE_CHUNKS
}

fn stress_above_threshold_residual_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_above_threshold_bvp_fixture_bindings::RESIDUAL_CHUNKS
}

fn stress_above_threshold_sparse_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_above_threshold_bvp_fixture_bindings::SPARSE_VALUE_CHUNKS
}

#[test]
fn real_bvp_parallel2_and_generated_aot_match_exactly_on_same_inputs() {
    let mut jac = build_real_bvp_damp1_case(8);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.2 + index as f64 * 0.01)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    let residual_generated = eval_generated_bvp_residual_col(flat_args.as_slice());
    let residual_generated_dense: Vec<f64> =
        residual_generated.iter().map(|value| *value).collect();

    assert_eq!(residual_parallel2.len(), residual_generated_dense.len());
    for (lhs, rhs) in residual_parallel2
        .iter()
        .zip(residual_generated_dense.iter())
    {
        let lhs = *lhs;
        let rhs = *rhs;
        assert!(
            (lhs - rhs).abs() < 1e-11,
            "residual mismatch: {lhs} vs {rhs}"
        );
    }

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_runtime_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });

    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    let jacobian_generated_values = eval_generated_bvp_sparse_values(flat_args.as_slice());
    let jacobian_generated = sparse_runtime_plan
        .assemble_sparse_col_mat(jacobian_generated_values.as_slice())
        .to_DMatrixType();

    assert_eq!(jacobian_parallel2.shape(), jacobian_generated.shape());
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-11,
                "jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

#[test]
fn benchmark_real_bvp_pipeline_stage_timings() {
    const SAMPLES: usize = 5;
    let mut symbolic_build = Vec::with_capacity(SAMPLES);
    let mut parallel2_residual_compile = Vec::with_capacity(SAMPLES);
    let mut parallel2_jacobian_compile = Vec::with_capacity(SAMPLES);
    let mut sparse_extract = Vec::with_capacity(SAMPLES);
    let mut residual_plan = Vec::with_capacity(SAMPLES);
    let mut sparse_plan = Vec::with_capacity(SAMPLES);
    let mut module_assembly = Vec::with_capacity(SAMPLES);
    let mut emit_source = Vec::with_capacity(SAMPLES);

    for _ in 0..SAMPLES {
        let start = Instant::now();
        let mut jac = build_real_bvp_damp1_case(8);
        symbolic_build.push(start.elapsed());

        let variable_names_owned = jac.variable_string.clone();
        let variable_names: Vec<&str> = variable_names_owned
            .iter()
            .map(|name| name.as_str())
            .collect();

        let start = Instant::now();
        jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
        parallel2_residual_compile.push(start.elapsed());

        let start = Instant::now();
        jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());
        parallel2_jacobian_compile.push(start.elapsed());

        let start = Instant::now();
        let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
        sparse_extract.push(start.elapsed());

        let residual_task = ResidualTask {
            fn_name: "fixture_bvp_residual",
            residuals: &jac.vector_of_functions,
            variables: &variable_names,
            params: None,
        };

        let start = Instant::now();
        let residual_runtime =
            residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
                max_outputs_per_chunk: 8,
            });
        residual_plan.push(start.elapsed());

        let sparse_task = SparseJacobianTask {
            fn_name: "fixture_bvp_sparse_values",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };

        let start = Instant::now();
        let sparse_runtime =
            sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });
        sparse_plan.push(start.elapsed());

        let start = Instant::now();
        let mut module = CodegenModule::new("generated_bvp_fixture_perf");
        for chunk in &residual_runtime.chunks {
            module.push_residual_block_plan(&chunk.plan);
        }
        for chunk in &sparse_runtime.chunks {
            module.push_sparse_values_plan(&chunk.plan);
        }
        module_assembly.push(start.elapsed());

        let start = Instant::now();
        let source = module.emit_source();
        emit_source.push(start.elapsed());
        black_box(source.len());
    }

    println!("\n=== Real BVP AOT / parallel2 stage timings (median of {SAMPLES}) ===");
    println!(
        "symbolic build (discretize + symbolic Jacobian): {:?}",
        median_duration(symbolic_build)
    );
    println!(
        "parallel2 residual compile: {:?}",
        median_duration(parallel2_residual_compile)
    );
    println!(
        "parallel2 sparse Jacobian compile: {:?}",
        median_duration(parallel2_jacobian_compile)
    );
    println!(
        "AOT sparse entry extraction: {:?}",
        median_duration(sparse_extract)
    );
    println!(
        "AOT residual runtime plan: {:?}",
        median_duration(residual_plan)
    );
    println!(
        "AOT sparse runtime plan: {:?}",
        median_duration(sparse_plan)
    );
    println!(
        "AOT codegen module assembly: {:?}",
        median_duration(module_assembly)
    );
    println!("AOT emit_source(): {:?}", median_duration(emit_source));
}

#[test]
fn benchmark_sparse_jacobian_chunking_strategies_on_real_bvp() {
    const SAMPLES: usize = 7;

    let jac = build_real_bvp_damp1_case(8);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };

    let strategies = [
        SparseChunkingStrategy::Whole,
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 },
        SparseChunkingStrategy::ByNonZeroCount {
            max_entries_per_chunk: 3,
        },
    ];

    println!("\n=== Sparse Jacobian chunking strategies on real BVP ===");
    for strategy in strategies {
        let mut planning_samples = Vec::with_capacity(SAMPLES);
        let mut assembly_samples = Vec::with_capacity(SAMPLES);
        let mut emit_samples = Vec::with_capacity(SAMPLES);

        let runtime_plan = sparse_task.runtime_plan(strategy);
        let chunk_count = runtime_plan.chunks.len();
        let max_chunk_nnz = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.entries.len())
            .max()
            .unwrap_or(0);
        let min_chunk_nnz = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.entries.len())
            .min()
            .unwrap_or(0);

        for _ in 0..SAMPLES {
            let start = Instant::now();
            let plan = sparse_task.runtime_plan(strategy);
            black_box(plan.chunks.len());
            planning_samples.push(start.elapsed());

            let start = Instant::now();
            let mut module = CodegenModule::new("generated_strategy_probe");
            for chunk in &runtime_plan.chunks {
                module.push_sparse_values_plan(&chunk.plan);
            }
            assembly_samples.push(start.elapsed());

            let start = Instant::now();
            let source = module.emit_source();
            black_box(source.len());
            emit_samples.push(start.elapsed());
        }

        let mut module = CodegenModule::new("generated_strategy_probe");
        for chunk in &runtime_plan.chunks {
            module.push_sparse_values_plan(&chunk.plan);
        }
        let source_len = module.emit_source().len();

        println!(
            "{}: chunks={}, nnz/chunk min..max={}..{}, planning={:?}, assembly={:?}, emit={:?}, source_len={}",
            sparse_strategy_label(strategy),
            chunk_count,
            min_chunk_nnz,
            max_chunk_nnz,
            median_duration(planning_samples),
            median_duration(assembly_samples),
            median_duration(emit_samples),
            source_len,
        );
    }
}

#[derive(tabled::Tabled)]
struct BandedPerfRow {
    scenario: String,
    path: String,
    compile_ms: String,
    runtime_ns: String,
    jacobian_diff: String,
}

#[derive(tabled::Tabled)]
struct LinearSolvePerfRow {
    scenario: String,
    path: String,
    prep_ms: String,
    solver_backend: String,
    layout: String,
    native_block_status: String,
    native_block_error: String,
    total_ns: String,
    solve_diff: String,
}

struct LinearSolveCorrectnessCheck {
    scenario: String,
    path: String,
    diff: f64,
    tol: f64,
}

struct BandedSolveDiagnostics {
    solution: Vec<f64>,
    solver_backend: String,
    layout: NodeMajorLayout,
    native_block_status: String,
    native_block_error: String,
}

#[test]
fn benchmark_bvp_sparse_vs_banded_lambdify_and_ir_table() {
    let scenarios: [(&str, fn() -> Jacobian, usize, usize); 3] = [
        ("small-damp1-24", || build_real_bvp_damp1_case(24), 25, 7),
        ("combustion-100", || build_combustion_bvp_case(100), 7, 5),
        ("combustion-1000", || build_combustion_bvp_case(1000), 3, 3),
    ];
    let mut rows = Vec::new();

    for (label, build_case, iters, samples) in scenarios {
        let jac = build_case();
        let prepared = jac.prepare_sparse_aot_problem(
            "bench_residual",
            "bench_sparse_values",
            ResidualChunkingStrategy::Whole,
            SparseChunkingStrategy::Whole,
        );
        let sparse_prepared = prepared.as_prepared_problem();
        let banded_prepared = prepared.as_prepared_banded_problem();
        let args: Vec<f64> = (0..sparse_prepared.input_names().len())
            .map(|index| 0.25 + index as f64 * 1.0e-4)
            .collect();

        let sparse_entries = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_compile_begin = Instant::now();
        let sparse_fns = compile_sparse_lambdify(sparse_prepared.input_names(), &sparse_entries);
        let sparse_compile_ms = sparse_compile_begin.elapsed().as_secs_f64() * 1_000.0;

        let banded_compile_begin = Instant::now();
        let banded_fns = compile_banded_lambdify(&prepared);
        let banded_compile_ms = banded_compile_begin.elapsed().as_secs_f64() * 1_000.0;

        let sparse_ir_build_begin = Instant::now();
        let sparse_blocks = sparse_ir_blocks(&sparse_prepared.jacobian_plan);
        let sparse_ir_build_ms = sparse_ir_build_begin.elapsed().as_secs_f64() * 1_000.0;

        let banded_ir_build_begin = Instant::now();
        let banded_blocks = banded_ir_blocks(&banded_prepared.jacobian_plan);
        let banded_ir_build_ms = banded_ir_build_begin.elapsed().as_secs_f64() * 1_000.0;

        let sparse_lambdify_dense = sparse_prepared
            .jacobian_plan
            .assemble_sparse_col_mat(&eval_functions(&sparse_fns, &args))
            .to_DMatrixType();
        let banded_lambdify_dense = banded_assembly_to_dense(
            &banded_prepared
                .jacobian_plan
                .assemble_banded_assembly(&eval_functions(&banded_fns, &args)),
        );
        let sparse_ir_dense = sparse_prepared
            .jacobian_plan
            .assemble_sparse_col_mat(&eval_ir_blocks_sequential(
                args.as_slice(),
                &sparse_blocks,
                sparse_prepared.jacobian_plan.nnz(),
            ))
            .to_DMatrixType();
        let banded_ir_dense =
            banded_assembly_to_dense(&banded_prepared.jacobian_plan.assemble_banded_assembly(
                &eval_ir_blocks_sequential(
                    args.as_slice(),
                    &banded_blocks,
                    banded_prepared.jacobian_plan.nnz(),
                ),
            ));

        let banded_lambdify_diff = max_abs_diff(
            sparse_lambdify_dense.as_slice(),
            banded_lambdify_dense.as_slice(),
        );
        let sparse_ir_diff =
            max_abs_diff(sparse_lambdify_dense.as_slice(), sparse_ir_dense.as_slice());
        let banded_ir_diff =
            max_abs_diff(sparse_lambdify_dense.as_slice(), banded_ir_dense.as_slice());

        let sparse_lambdify_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_functions(&sparse_fns, &args);
                            black_box(
                                sparse_prepared
                                    .jacobian_plan
                                    .assemble_sparse_col_mat(values.as_slice()),
                            );
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_lambdify_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_functions(&banded_fns, &args);
                            black_box(
                                banded_prepared
                                    .jacobian_plan
                                    .assemble_banded_assembly(values.as_slice()),
                            );
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let sparse_ir_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_ir_blocks_sequential(
                                args.as_slice(),
                                &sparse_blocks,
                                sparse_prepared.jacobian_plan.nnz(),
                            );
                            black_box(
                                sparse_prepared
                                    .jacobian_plan
                                    .assemble_sparse_col_mat(values.as_slice()),
                            );
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_ir_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_ir_blocks_sequential(
                                args.as_slice(),
                                &banded_blocks,
                                banded_prepared.jacobian_plan.nnz(),
                            );
                            black_box(
                                banded_prepared
                                    .jacobian_plan
                                    .assemble_banded_assembly(values.as_slice()),
                            );
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );

        rows.push(BandedPerfRow {
            scenario: label.to_string(),
            path: "Lambdify-Sparse".to_string(),
            compile_ms: format!("{:.3}", sparse_compile_ms),
            runtime_ns: format!("{:.0}", sparse_lambdify_ns),
            jacobian_diff: format!("{:.3e}", 0.0),
        });
        rows.push(BandedPerfRow {
            scenario: label.to_string(),
            path: "Lambdify-Banded".to_string(),
            compile_ms: format!("{:.3}", banded_compile_ms),
            runtime_ns: format!("{:.0}", banded_lambdify_ns),
            jacobian_diff: format!("{:.3e}", banded_lambdify_diff),
        });
        rows.push(BandedPerfRow {
            scenario: label.to_string(),
            path: "AOTIR-Sparse".to_string(),
            compile_ms: format!("{:.3}", sparse_ir_build_ms),
            runtime_ns: format!("{:.0}", sparse_ir_ns),
            jacobian_diff: format!("{:.3e}", sparse_ir_diff),
        });
        rows.push(BandedPerfRow {
            scenario: label.to_string(),
            path: "AOTIR-Banded".to_string(),
            compile_ms: format!("{:.3}", banded_ir_build_ms),
            runtime_ns: format!("{:.0}", banded_ir_ns),
            jacobian_diff: format!("{:.3e}", banded_ir_diff),
        });

        let tol = if label == "combustion-1000" {
            1.0e-9
        } else {
            1.0e-11
        };
        assert!(banded_lambdify_diff <= tol);
        assert!(sparse_ir_diff <= tol);
        assert!(banded_ir_diff <= tol);
    }

    let mut builder = Builder::default();
    builder.push_record([
        "scenario",
        "path",
        "compile_ms",
        "runtime_ns",
        "jacobian_diff",
    ]);
    for row in rows {
        builder.push_record([
            row.scenario,
            row.path,
            row.compile_ms,
            row.runtime_ns,
            row.jacobian_diff,
        ]);
    }
    let mut table = builder.build();
    table.with(Style::rounded());
    println!("[BVP perf] sparse vs banded jacobian runtime table");
    println!("{table}");
}

#[test]
fn benchmark_bvp_sparse_vs_banded_linear_solve_table() {
    let scenarios: [(&str, fn() -> Jacobian, usize, usize); 3] = [
        ("small-damp1-24", || build_real_bvp_damp1_case(24), 15, 5),
        ("combustion-100", || build_combustion_bvp_case(100), 4, 4),
        ("combustion-1000", || build_combustion_bvp_case(1000), 2, 2),
    ];
    let mut rows = Vec::new();
    let mut checks = Vec::new();

    for (label, build_case, iters, samples) in scenarios {
        let jac = build_case();
        let prepared = jac.prepare_sparse_aot_problem(
            "bench_linear_residual",
            "bench_linear_sparse_values",
            ResidualChunkingStrategy::Whole,
            SparseChunkingStrategy::Whole,
        );
        let sparse_prepared = prepared.as_prepared_problem();
        let banded_prepared = prepared.as_prepared_banded_problem();
        let args: Vec<f64> = (0..sparse_prepared.input_names().len())
            .map(|index| 0.25 + index as f64 * 1.0e-4)
            .collect();
        let rhs = sample_rhs(sparse_prepared.jacobian_plan.structure.rows);

        let sparse_entries = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_compile_begin = Instant::now();
        let sparse_fns = compile_sparse_lambdify(sparse_prepared.input_names(), &sparse_entries);
        let sparse_compile_ms = sparse_compile_begin.elapsed().as_secs_f64() * 1_000.0;

        let banded_compile_begin = Instant::now();
        let banded_fns = compile_banded_lambdify(&prepared);
        let banded_compile_ms = banded_compile_begin.elapsed().as_secs_f64() * 1_000.0;

        let sparse_ir_build_begin = Instant::now();
        let sparse_blocks = sparse_ir_blocks(&sparse_prepared.jacobian_plan);
        let sparse_ir_build_ms = sparse_ir_build_begin.elapsed().as_secs_f64() * 1_000.0;

        let banded_ir_build_begin = Instant::now();
        let banded_blocks = banded_ir_blocks(&banded_prepared.jacobian_plan);
        let banded_ir_build_ms = banded_ir_build_begin.elapsed().as_secs_f64() * 1_000.0;

        let sparse_baseline_values = eval_functions(&sparse_fns, &args);
        let sparse_baseline_matrix = sparse_prepared
            .jacobian_plan
            .assemble_sparse_col_mat(sparse_baseline_values.as_slice());
        let baseline_solution = solve_sparse_lu(&sparse_baseline_matrix, rhs.as_slice());

        let banded_lambdify_diagnostics = solve_banded_block_tridiagonal(
            &jac,
            &banded_prepared
                .jacobian_plan
                .assemble_banded_assembly(eval_functions(&banded_fns, &args).as_slice()),
            rhs.as_slice(),
        );
        let banded_lambdify_consistent_diagnostics = solve_banded_block_tridiagonal_consistent(
            &jac,
            &banded_prepared
                .jacobian_plan
                .assemble_banded_assembly(eval_functions(&banded_fns, &args).as_slice()),
            rhs.as_slice(),
            1,
        );
        let sparse_ir_solution = solve_sparse_lu(
            &sparse_prepared.jacobian_plan.assemble_sparse_col_mat(
                eval_ir_blocks_sequential(
                    args.as_slice(),
                    &sparse_blocks,
                    sparse_prepared.jacobian_plan.nnz(),
                )
                .as_slice(),
            ),
            rhs.as_slice(),
        );
        let banded_ir_diagnostics = solve_banded_block_tridiagonal(
            &jac,
            &banded_prepared.jacobian_plan.assemble_banded_assembly(
                eval_ir_blocks_sequential(
                    args.as_slice(),
                    &banded_blocks,
                    banded_prepared.jacobian_plan.nnz(),
                )
                .as_slice(),
            ),
            rhs.as_slice(),
        );
        let banded_ir_consistent_diagnostics = solve_banded_block_tridiagonal_consistent(
            &jac,
            &banded_prepared.jacobian_plan.assemble_banded_assembly(
                eval_ir_blocks_sequential(
                    args.as_slice(),
                    &banded_blocks,
                    banded_prepared.jacobian_plan.nnz(),
                )
                .as_slice(),
            ),
            rhs.as_slice(),
            1,
        );

        let banded_lambdify_diff =
            max_abs_diff(&baseline_solution, &banded_lambdify_diagnostics.solution);
        let banded_lambdify_consistent_diff = max_abs_diff(
            &baseline_solution,
            &banded_lambdify_consistent_diagnostics.solution,
        );
        let sparse_ir_diff = max_abs_diff(&baseline_solution, &sparse_ir_solution);
        let banded_ir_diff = max_abs_diff(&baseline_solution, &banded_ir_diagnostics.solution);
        let banded_ir_consistent_diff = max_abs_diff(
            &baseline_solution,
            &banded_ir_consistent_diagnostics.solution,
        );

        let sparse_lambdify_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_functions(&sparse_fns, &args);
                            let matrix = sparse_prepared
                                .jacobian_plan
                                .assemble_sparse_col_mat(values.as_slice());
                            black_box(solve_sparse_lu(&matrix, rhs.as_slice()));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_lambdify_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_functions(&banded_fns, &args);
                            let assembly = banded_prepared
                                .jacobian_plan
                                .assemble_banded_assembly(values.as_slice());
                            black_box(solve_banded_block_tridiagonal(
                                &jac,
                                &assembly,
                                rhs.as_slice(),
                            ));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_lambdify_consistent_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_functions(&banded_fns, &args);
                            let assembly = banded_prepared
                                .jacobian_plan
                                .assemble_banded_assembly(values.as_slice());
                            black_box(solve_banded_block_tridiagonal_consistent(
                                &jac,
                                &assembly,
                                rhs.as_slice(),
                                1,
                            ));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let sparse_ir_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_ir_blocks_sequential(
                                args.as_slice(),
                                &sparse_blocks,
                                sparse_prepared.jacobian_plan.nnz(),
                            );
                            let matrix = sparse_prepared
                                .jacobian_plan
                                .assemble_sparse_col_mat(values.as_slice());
                            black_box(solve_sparse_lu(&matrix, rhs.as_slice()));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_ir_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_ir_blocks_sequential(
                                args.as_slice(),
                                &banded_blocks,
                                banded_prepared.jacobian_plan.nnz(),
                            );
                            let assembly = banded_prepared
                                .jacobian_plan
                                .assemble_banded_assembly(values.as_slice());
                            black_box(solve_banded_block_tridiagonal(
                                &jac,
                                &assembly,
                                rhs.as_slice(),
                            ));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );
        let banded_ir_consistent_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..iters {
                            let values = eval_ir_blocks_sequential(
                                args.as_slice(),
                                &banded_blocks,
                                banded_prepared.jacobian_plan.nnz(),
                            );
                            let assembly = banded_prepared
                                .jacobian_plan
                                .assemble_banded_assembly(values.as_slice());
                            black_box(solve_banded_block_tridiagonal_consistent(
                                &jac,
                                &assembly,
                                rhs.as_slice(),
                                1,
                            ));
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            iters,
        );

        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "Lambdify-Sparse".to_string(),
            prep_ms: format!("{:.3}", sparse_compile_ms),
            solver_backend: "faer_sparse_lu".to_string(),
            layout: "-".to_string(),
            native_block_status: "-".to_string(),
            native_block_error: "-".to_string(),
            total_ns: format!("{:.0}", sparse_lambdify_ns),
            solve_diff: format!("{:.3e}", 0.0),
        });
        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "Lambdify-Banded".to_string(),
            prep_ms: format!("{:.3}", banded_compile_ms),
            solver_backend: banded_lambdify_diagnostics.solver_backend,
            layout: format!(
                "{}x{}",
                banded_lambdify_diagnostics.layout.n_blocks(),
                banded_lambdify_diagnostics.layout.block_size()
            ),
            native_block_status: banded_lambdify_diagnostics.native_block_status,
            native_block_error: banded_lambdify_diagnostics.native_block_error,
            total_ns: format!("{:.0}", banded_lambdify_ns),
            solve_diff: format!("{:.3e}", banded_lambdify_diff),
        });
        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "Lambdify-Banded-Cons".to_string(),
            prep_ms: format!("{:.3}", banded_compile_ms),
            solver_backend: banded_lambdify_consistent_diagnostics.solver_backend,
            layout: format!(
                "{}x{}",
                banded_lambdify_consistent_diagnostics.layout.n_blocks(),
                banded_lambdify_consistent_diagnostics.layout.block_size()
            ),
            native_block_status: banded_lambdify_consistent_diagnostics.native_block_status,
            native_block_error: banded_lambdify_consistent_diagnostics.native_block_error,
            total_ns: format!("{:.0}", banded_lambdify_consistent_ns),
            solve_diff: format!("{:.3e}", banded_lambdify_consistent_diff),
        });
        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "AOTIR-Sparse".to_string(),
            prep_ms: format!("{:.3}", sparse_ir_build_ms),
            solver_backend: "faer_sparse_lu".to_string(),
            layout: "-".to_string(),
            native_block_status: "-".to_string(),
            native_block_error: "-".to_string(),
            total_ns: format!("{:.0}", sparse_ir_ns),
            solve_diff: format!("{:.3e}", sparse_ir_diff),
        });
        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "AOTIR-Banded".to_string(),
            prep_ms: format!("{:.3}", banded_ir_build_ms),
            solver_backend: banded_ir_diagnostics.solver_backend,
            layout: format!(
                "{}x{}",
                banded_ir_diagnostics.layout.n_blocks(),
                banded_ir_diagnostics.layout.block_size()
            ),
            native_block_status: banded_ir_diagnostics.native_block_status,
            native_block_error: banded_ir_diagnostics.native_block_error,
            total_ns: format!("{:.0}", banded_ir_ns),
            solve_diff: format!("{:.3e}", banded_ir_diff),
        });
        rows.push(LinearSolvePerfRow {
            scenario: label.to_string(),
            path: "AOTIR-Banded-Cons".to_string(),
            prep_ms: format!("{:.3}", banded_ir_build_ms),
            solver_backend: banded_ir_consistent_diagnostics.solver_backend,
            layout: format!(
                "{}x{}",
                banded_ir_consistent_diagnostics.layout.n_blocks(),
                banded_ir_consistent_diagnostics.layout.block_size()
            ),
            native_block_status: banded_ir_consistent_diagnostics.native_block_status,
            native_block_error: banded_ir_consistent_diagnostics.native_block_error,
            total_ns: format!("{:.0}", banded_ir_consistent_ns),
            solve_diff: format!("{:.3e}", banded_ir_consistent_diff),
        });

        let tol = if label == "combustion-1000" {
            1.0e-8
        } else {
            1.0e-10
        };
        checks.push(LinearSolveCorrectnessCheck {
            scenario: label.to_string(),
            path: "Lambdify-Banded".to_string(),
            diff: banded_lambdify_diff,
            tol,
        });
        checks.push(LinearSolveCorrectnessCheck {
            scenario: label.to_string(),
            path: "Lambdify-Banded-Cons".to_string(),
            diff: banded_lambdify_consistent_diff,
            tol,
        });
        checks.push(LinearSolveCorrectnessCheck {
            scenario: label.to_string(),
            path: "AOTIR-Sparse".to_string(),
            diff: sparse_ir_diff,
            tol,
        });
        checks.push(LinearSolveCorrectnessCheck {
            scenario: label.to_string(),
            path: "AOTIR-Banded".to_string(),
            diff: banded_ir_diff,
            tol,
        });
        checks.push(LinearSolveCorrectnessCheck {
            scenario: label.to_string(),
            path: "AOTIR-Banded-Cons".to_string(),
            diff: banded_ir_consistent_diff,
            tol,
        });
    }

    let mut builder = Builder::default();
    builder.push_record([
        "scenario",
        "path",
        "prep_ms",
        "solver_backend",
        "layout",
        "native_block_status",
        "native_block_error",
        "total_ns",
        "solve_diff",
    ]);
    for row in rows {
        builder.push_record([
            row.scenario,
            row.path,
            row.prep_ms,
            row.solver_backend,
            row.layout,
            row.native_block_status,
            row.native_block_error,
            row.total_ns,
            row.solve_diff,
        ]);
    }
    let mut table = builder.build();
    table.with(Style::rounded());
    println!("[BVP perf] sparse vs banded linear solve table");
    println!("{table}");

    let failing: Vec<_> = checks
        .into_iter()
        .filter(|check| check.diff > check.tol)
        .collect();
    if !failing.is_empty() {
        for check in &failing {
            println!(
                "[BVP perf diagnostic] correctness mismatch scenario={} path={} diff={:.6e} tol={:.6e}",
                check.scenario, check.path, check.diff, check.tol
            );
        }
    }
}

fn block_column_nonzero_counts(block: &[f64], bs: usize) -> Vec<usize> {
    (0..bs)
        .map(|col| {
            (0..bs)
                .filter(|&row| block[row * bs + col].abs() > 1.0e-14)
                .count()
        })
        .collect()
}

fn block_column_max_abs(block: &[f64], bs: usize) -> Vec<f64> {
    (0..bs)
        .map(|col| {
            (0..bs)
                .map(|row| block[row * bs + col].abs())
                .fold(0.0_f64, f64::max)
        })
        .collect()
}

fn format_block(block: &[f64], bs: usize) -> String {
    let mut out = String::new();
    for row in 0..bs {
        let line = (0..bs)
            .map(|col| format!("{:>12.4e}", block[row * bs + col]))
            .collect::<Vec<_>>()
            .join(" ");
        out.push_str(&line);
        out.push('\n');
    }
    out
}

fn format_block_column_stats(block: &[f64], bs: usize) -> String {
    let nnz = block_column_nonzero_counts(block, bs);
    let max_abs = block_column_max_abs(block, bs);
    nnz.iter()
        .zip(max_abs.iter())
        .map(|(nnz, max_abs)| format!("{nnz}:{max_abs:.2e}"))
        .collect::<Vec<_>>()
        .join(", ")
}

fn dense_linf_diff(lhs: &[Vec<f64>], rhs: &[Vec<f64>]) -> f64 {
    lhs.iter()
        .zip(rhs.iter())
        .flat_map(|(lhs_row, rhs_row)| lhs_row.iter().zip(rhs_row.iter()))
        .map(|(lhs, rhs)| (lhs - rhs).abs())
        .fold(0.0_f64, f64::max)
}

fn dense_matmul_square(lhs: &[Vec<f64>], rhs: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let n = lhs.len();
    let mut out = vec![vec![0.0; n]; n];
    for i in 0..n {
        for k in 0..n {
            let lhs_ik = lhs[i][k];
            if lhs_ik == 0.0 {
                continue;
            }
            for j in 0..n {
                out[i][j] += lhs_ik * rhs[k][j];
            }
        }
    }
    out
}

fn permutation_from_swap_pivots(pivots: &[usize]) -> Vec<usize> {
    let mut perm: Vec<usize> = (0..pivots.len()).collect();
    for (k, &p) in pivots.iter().enumerate() {
        perm.swap(k, p);
    }
    perm
}

fn apply_block_row_permutations(
    dense: &[Vec<f64>],
    block_size: usize,
    pivots_by_block: &[Vec<usize>],
) -> Vec<Vec<f64>> {
    let n = dense.len();
    let mut out = dense.to_vec();
    for (block_idx, pivots) in pivots_by_block.iter().enumerate() {
        let perm = permutation_from_swap_pivots(pivots);
        let row0 = block_idx * block_size;
        let mut permuted_rows = vec![vec![0.0; n]; block_size];
        for local_row in 0..block_size {
            permuted_rows[local_row] = dense[row0 + perm[local_row]].clone();
        }
        for local_row in 0..block_size {
            out[row0 + local_row] = permuted_rows[local_row].clone();
        }
    }
    out
}

fn split_dense_lu_block(lu: &[f64], bs: usize) -> (Vec<f64>, Vec<f64>) {
    let mut l = vec![0.0; bs * bs];
    let mut u = vec![0.0; bs * bs];
    for i in 0..bs {
        for j in 0..bs {
            let idx = i * bs + j;
            if i > j {
                l[idx] = lu[idx];
            } else if i == j {
                l[idx] = 1.0;
                u[idx] = lu[idx];
            } else {
                u[idx] = lu[idx];
            }
        }
    }
    (l, u)
}

fn reconstruct_block_lu_factors(lu: &BlockTridiagonalLu) -> (Vec<Vec<f64>>, Vec<Vec<f64>>) {
    let nb = lu.n_blocks();
    let bs = lu.block_size();
    let n = lu.n();
    let mut l = vec![vec![0.0; n]; n];
    let mut u = vec![vec![0.0; n]; n];

    for block_idx in 0..nb {
        let row0 = block_idx * bs;
        let col0 = block_idx * bs;
        let (l_diag, u_diag) = split_dense_lu_block(&lu.diag_lu_blocks()[block_idx], bs);

        for i in 0..bs {
            for j in 0..bs {
                l[row0 + i][col0 + j] = l_diag[i * bs + j];
                u[row0 + i][col0 + j] = u_diag[i * bs + j];
            }
        }

        if block_idx + 1 < nb {
            let lower = &lu.lower_mult_blocks()[block_idx];
            let upper = &lu.upper_blocks()[block_idx];
            let lower_row0 = (block_idx + 1) * bs;
            let lower_col0 = block_idx * bs;
            let upper_row0 = block_idx * bs;
            let upper_col0 = (block_idx + 1) * bs;

            for i in 0..bs {
                for j in 0..bs {
                    l[lower_row0 + i][lower_col0 + j] = lower[i * bs + j];
                    u[upper_row0 + i][upper_col0 + j] = upper[i * bs + j];
                }
            }
        }
    }

    (l, u)
}

/// Focused structural diagnostic for the combustion banded path.
///
/// This is intentionally ignored by default because it prints large block data.
/// We use it when the linear-solve table reports a native block fallback and we
/// want to inspect the first diagonal block under the current node-major layout.
pub(crate) fn run_diagnose_combustion_banded_native_zero_pivot() {
    for n in [100usize, 1000usize] {
        let jac = build_combustion_bvp_case(n);
        let prepared = jac.prepare_sparse_aot_problem(
            "diag_combustion_residual",
            "diag_combustion_sparse_values",
            ResidualChunkingStrategy::Whole,
            SparseChunkingStrategy::Whole,
        );
        let banded_prepared = prepared.as_prepared_banded_problem();
        let args: Vec<f64> = (0..banded_prepared.input_names().len())
            .map(|index| 0.25 + index as f64 * 1.0e-4)
            .collect();
        let banded_fns = compile_banded_lambdify(&prepared);
        let assembly = banded_prepared
            .jacobian_plan
            .assemble_banded_assembly(eval_functions(&banded_fns, &args).as_slice());
        let layout = infer_node_major_layout(&jac, assembly.n());
        let block = assembly
            .to_block_tridiagonal(layout.n_blocks(), layout.block_size())
            .expect("combustion assembly should convert to block-tridiagonal");

        let first_diag = block.diag_block(0).expect("first diagonal block exists");
        let first_upper = block.upper_block(0).expect("first upper block exists");

        println!(
            "[BVP banded diagnostic] combustion-{n}: layout={}x{}",
            layout.n_blocks(),
            layout.block_size()
        );
        println!(
            "[BVP banded diagnostic] combustion-{n}: first diag column nnz = {:?}",
            block_column_nonzero_counts(first_diag, layout.block_size())
        );
        println!(
            "[BVP banded diagnostic] combustion-{n}: first diag column max_abs = {:?}",
            block_column_max_abs(first_diag, layout.block_size())
        );
        println!(
            "[BVP banded diagnostic] combustion-{n}: first diag block\n{}",
            format_block(first_diag, layout.block_size())
        );
        println!(
            "[BVP banded diagnostic] combustion-{n}: first upper block\n{}",
            format_block(first_upper, layout.block_size())
        );
    }
}

/// Superblock experiment for the combustion banded path.
///
/// The goal is not to benchmark speed, but to answer one structural question:
/// does grouping several neighboring nodes into one solver block remove the
/// local singularity that makes the plain `6x6` node-major block fail?
pub(crate) fn run_diagnose_combustion_superblock_groupings() {
    #[derive(Debug)]
    struct SuperblockRow {
        scenario: String,
        group: String,
        solver: String,
        layout: String,
        first_diag_cols: String,
        status: String,
        error: String,
        factor_residual: String,
        solve_diff: String,
        block_residual: String,
    }

    let mut rows = Vec::new();

    for n in [100usize, 1000usize] {
        let jac = build_combustion_bvp_case(n);
        let prepared = jac.prepare_sparse_aot_problem(
            "diag_combustion_superblock_residual",
            "diag_combustion_superblock_sparse_values",
            ResidualChunkingStrategy::Whole,
            SparseChunkingStrategy::Whole,
        );
        let input_names = prepared.as_prepared_problem().input_names().to_vec();
        let args: Vec<f64> = (0..input_names.len())
            .map(|index| 0.25 + index as f64 * 1.0e-4)
            .collect();

        let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_fns =
            compile_sparse_lambdify(input_names.as_slice(), sparse_entries_owned.as_slice());
        let sparse_values = eval_functions(&sparse_fns, &args);
        let sparse_prepared = prepared.as_prepared_problem();
        let sparse_matrix = sparse_prepared
            .jacobian_plan
            .assemble_sparse_col_mat(sparse_values.as_slice());
        let rhs = sample_rhs(sparse_matrix.nrows());
        let sparse_solution = solve_sparse_lu(&sparse_matrix, rhs.as_slice());

        let banded_prepared = prepared.as_prepared_banded_problem();
        let banded_fns = compile_banded_lambdify(&prepared);
        let banded_values = eval_functions(&banded_fns, &args);
        let assembly = banded_prepared
            .jacobian_plan
            .assemble_banded_assembly(banded_values.as_slice());

        let node_layout = infer_node_major_layout(&jac, assembly.n());
        let n_nodes = node_layout.n_nodes();
        let vars_per_node = node_layout.vars_per_node();

        for nodes_per_superblock in [1usize, 2, 3, 4] {
            if n_nodes % nodes_per_superblock != 0 {
                rows.push(SuperblockRow {
                    scenario: format!("combustion-{n}"),
                    group: nodes_per_superblock.to_string(),
                    solver: "n/a".to_string(),
                    layout: format!("skip/{vars_per_node}"),
                    first_diag_cols: "-".to_string(),
                    status: "skipped".to_string(),
                    error: "n_nodes is not divisible by group".to_string(),
                    factor_residual: "-".to_string(),
                    solve_diff: "-".to_string(),
                    block_residual: "-".to_string(),
                });
                continue;
            }

            let layout = SuperBlockLayout::new(n_nodes, vars_per_node, nodes_per_superblock)
                .expect("superblock layout should be valid");
            let block = assembly
                .finalize_superblock_tridiagonal(&layout)
                .expect("superblock layout should convert to uniform block-tridiagonal storage");
            let first_diag = block.diag_block(0).expect("first diagonal block exists");

            let mut legacy_lu = BlockTridiagonalLu::new(block.n_blocks(), block.block_size())
                .expect("native block LU workspace should be constructible");
            let mut consistent_lu =
                BlockTridiagonalLuConsistent::new(block.n_blocks(), block.block_size())
                    .expect("consistent block LU workspace should be constructible");

            rows.push(SuperblockRow {
                scenario: format!("combustion-{n}"),
                group: nodes_per_superblock.to_string(),
                solver: "legacy".to_string(),
                layout: format!("{}x{}", layout.n_blocks(), layout.block_size()),
                first_diag_cols: format_block_column_stats(first_diag, layout.block_size()),
                status: match legacy_lu.factor_from(&block) {
                    Ok(()) => "ok".to_string(),
                    Err(_) => "failed".to_string(),
                },
                error: "-".to_string(),
                factor_residual: "-".to_string(),
                solve_diff: "-".to_string(),
                block_residual: "-".to_string(),
            });

            if let Err(err) = legacy_lu.factor_from(&block) {
                if let Some(last) = rows.last_mut() {
                    last.error = format!("{err:?}");
                }
            } else if let Some(last) = rows.last_mut() {
                let dense_a = block.to_dense();
                let dense_pa = apply_block_row_permutations(
                    dense_a.as_slice(),
                    block.block_size(),
                    legacy_lu.diag_pivots(),
                );
                let (dense_l, dense_u) = reconstruct_block_lu_factors(&legacy_lu);
                let dense_lu = dense_matmul_square(dense_l.as_slice(), dense_u.as_slice());
                let mut rhs_native = rhs.clone();
                legacy_lu.solve_in_place(rhs_native.as_mut_slice()).expect(
                    "legacy superblock solve should succeed after successful factorization",
                );
                let block_ax = crate::somelinalg::banded::ops::dense_matvec(
                    dense_a.as_slice(),
                    rhs_native.as_slice(),
                );
                last.factor_residual = format!(
                    "{:.3e}",
                    dense_linf_diff(dense_pa.as_slice(), dense_lu.as_slice())
                );
                last.solve_diff = format!("{:.3e}", max_abs_diff(&rhs_native, &sparse_solution));
                last.block_residual = format!(
                    "{:.3e}",
                    crate::somelinalg::banded::ops::vec_diff_linf(
                        block_ax.as_slice(),
                        rhs.as_slice(),
                    )
                );
            }

            rows.push(SuperblockRow {
                scenario: format!("combustion-{n}"),
                group: nodes_per_superblock.to_string(),
                solver: "consistent".to_string(),
                layout: format!("{}x{}", layout.n_blocks(), layout.block_size()),
                first_diag_cols: format_block_column_stats(first_diag, layout.block_size()),
                status: match consistent_lu.factor_from(&block) {
                    Ok(()) => "ok".to_string(),
                    Err(_) => "failed".to_string(),
                },
                error: "-".to_string(),
                factor_residual: "-".to_string(),
                solve_diff: "-".to_string(),
                block_residual: "-".to_string(),
            });
            if let Err(err) = consistent_lu.factor_from(&block) {
                if let Some(last) = rows.last_mut() {
                    last.error = format!("{err:?}");
                }
            } else if let Some(last) = rows.last_mut() {
                let dense_a = block.to_dense();
                let mut rhs_native = rhs.clone();
                consistent_lu
                    .solve_in_place(rhs_native.as_mut_slice())
                    .expect(
                        "consistent superblock solve should succeed after successful factorization",
                    );
                let block_ax = crate::somelinalg::banded::ops::dense_matvec(
                    dense_a.as_slice(),
                    rhs_native.as_slice(),
                );
                last.factor_residual = format!(
                    "{:.3e}",
                    consistent_lu
                        .factor_residual_inf(&block)
                        .expect("consistent factor residual should be available")
                );
                last.solve_diff = format!("{:.3e}", max_abs_diff(&rhs_native, &sparse_solution));
                last.block_residual = format!(
                    "{:.3e}",
                    crate::somelinalg::banded::ops::vec_diff_linf(
                        block_ax.as_slice(),
                        rhs.as_slice(),
                    )
                );
            }
        }
    }

    let mut builder = Builder::default();
    builder.push_record([
        "scenario",
        "group",
        "solver",
        "layout",
        "first_diag_cols",
        "status",
        "error",
        "factor_residual",
        "solve_diff",
        "block_residual",
    ]);
    for row in rows {
        builder.push_record([
            row.scenario,
            row.group,
            row.solver,
            row.layout,
            row.first_diag_cols,
            row.status,
            row.error,
            row.factor_residual,
            row.solve_diff,
            row.block_residual,
        ]);
    }
    let mut table = builder.build();
    table.with(Style::rounded());
    println!("[BVP superblock diagnostic] combustion superblock grouping table");
    println!("{table}");
}

#[test]
fn benchmark_real_bvp_parallel2_vs_generated_aot_hot_path() {
    const SAMPLES: usize = 7;
    const RESIDUAL_ITERS: usize = 800;
    const JACOBIAN_ITERS: usize = 250;

    let mut jac = build_real_bvp_damp1_case(8);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.2 + index as f64 * 0.01)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_expected = jac.residiual_function.call(0.0, &input_col);
    let residual_generated = eval_generated_bvp_residual_col(flat_args.as_slice());
    let residual_expected_dense = residual_expected.to_DVectorType();
    let residual_generated_dense: Vec<f64> =
        residual_generated.iter().map(|value| *value).collect();
    assert_eq!(
        residual_expected_dense.len(),
        residual_generated_dense.len()
    );
    for (lhs, rhs) in residual_expected_dense
        .iter()
        .zip(residual_generated_dense.iter())
    {
        let lhs = *lhs;
        let rhs = *rhs;
        assert!(
            (lhs - rhs).abs() < 1e-11,
            "residual mismatch: {lhs} vs {rhs}"
        );
    }

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let generated_sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });
    let mut jacobian_fun = jac
        .jac_function
        .take()
        .expect("parallel2 Jacobian should be compiled");
    let jacobian_expected = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
    let jacobian_generated_values = eval_generated_bvp_sparse_values(flat_args.as_slice());
    let jacobian_generated = generated_sparse_plan
        .assemble_sparse_col_mat(jacobian_generated_values.as_slice())
        .to_DMatrixType();
    assert_eq!(jacobian_expected.shape(), jacobian_generated.shape());
    for i in 0..jacobian_expected.nrows() {
        for j in 0..jacobian_expected.ncols() {
            let lhs = jacobian_expected[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-11,
                "jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }

    let mut parallel2_residual_samples = Vec::with_capacity(SAMPLES);
    let mut aot_residual_samples = Vec::with_capacity(SAMPLES);
    let mut parallel2_jacobian_samples = Vec::with_capacity(SAMPLES);
    let mut aot_jacobian_values_samples = Vec::with_capacity(SAMPLES);
    let mut aot_jacobian_solver_samples = Vec::with_capacity(SAMPLES);

    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = jac.residiual_function.call(0.0, &input_col);
            black_box(out.get_val(0));
            black_box(out.len());
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = eval_generated_bvp_residual_col(flat_args.as_slice());
            black_box(out[0]);
            black_box(out.len());
        }
        aot_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = jacobian_fun.call(0.0, &input_col);
            black_box(out.shape());
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let values = eval_generated_bvp_sparse_values(flat_args.as_slice());
            black_box(values[0]);
            black_box(values.len());
        }
        aot_jacobian_values_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let values = eval_generated_bvp_sparse_values(flat_args.as_slice());
            let mat = generated_sparse_plan.assemble_sparse_col_mat(values.as_slice());
            black_box((mat.nrows(), mat.ncols(), mat.compute_nnz()));
        }
        aot_jacobian_solver_samples.push(start.elapsed());
    }

    let parallel2_residual = median_duration(parallel2_residual_samples);
    let aot_residual = median_duration(aot_residual_samples);
    let parallel2_jacobian = median_duration(parallel2_jacobian_samples);
    let aot_jacobian_values = median_duration(aot_jacobian_values_samples);
    let aot_jacobian_solver = median_duration(aot_jacobian_solver_samples);

    println!("\n=== Real BVP hot-path comparison (median of {SAMPLES}) ===");
    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(parallel2_residual, RESIDUAL_ITERS)
    );
    println!(
        "AOT generated residual chunks -> Col: {:.2} ns/call",
        per_iter_ns(aot_residual, RESIDUAL_ITERS)
    );
    println!(
        "parallel2 Jacobian -> SparseColMat: {:.2} ns/call",
        per_iter_ns(parallel2_jacobian, JACOBIAN_ITERS)
    );
    println!(
        "AOT generated Jacobian values only: {:.2} ns/call",
        per_iter_ns(aot_jacobian_values, JACOBIAN_ITERS)
    );
    println!(
        "AOT generated Jacobian values + SparseColMat assembly: {:.2} ns/call",
        per_iter_ns(aot_jacobian_solver, JACOBIAN_ITERS)
    );
}

#[test]
fn benchmark_real_bvp_orchestrator_runtime_paths() {
    const SAMPLES: usize = 7;
    const RESIDUAL_ITERS: usize = 800;
    const JACOBIAN_ITERS: usize = 250;

    let mut jac = build_real_bvp_damp1_case(8);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.2 + index as f64 * 0.01)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 8,
    });
    let residual_chunk_fns: [GeneratedChunkFn; 2] = [
        generated_bvp_fixture::fixture_bvp_residual_chunk_0,
        generated_bvp_fixture::fixture_bvp_residual_chunk_1,
    ];
    let residual_bindings = bind_residual_fixture_chunks(&residual_plan, &residual_chunk_fns);
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });
    let sparse_chunk_fns: [GeneratedChunkFn; 16] = [
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_0,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_1,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_2,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_3,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_4,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_5,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_6,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_7,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_8,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_9,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_10,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_11,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_12,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_13,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_14,
        generated_bvp_fixture::fixture_bvp_sparse_values_chunk_15,
    ];
    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, &sparse_chunk_fns);
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    assert_eq!(
        residual_sequential.eval(flat_args.as_slice()),
        residual_parallel.eval(flat_args.as_slice())
    );
    assert_eq!(
        sparse_sequential.eval_values(flat_args.as_slice()),
        sparse_parallel.eval_values(flat_args.as_slice())
    );

    let mut parallel2_residual_samples = Vec::with_capacity(SAMPLES);
    let mut manual_residual_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_residual_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_residual_samples = Vec::with_capacity(SAMPLES);

    let mut parallel2_jacobian_samples = Vec::with_capacity(SAMPLES);
    let mut manual_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_sparse_samples = Vec::with_capacity(SAMPLES);

    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = jac.residiual_function.call(0.0, &input_col);
            black_box(out.get_val(0));
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = eval_generated_bvp_residual_values(flat_args.as_slice());
            black_box(out[0]);
        }
        manual_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = jac
                .jac_function
                .as_mut()
                .expect("parallel2 Jacobian should exist")
                .call(0.0, &input_col);
            black_box(out.shape());
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = eval_generated_bvp_sparse_values(flat_args.as_slice());
            black_box(out[0]);
        }
        manual_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!("\n=== Real BVP runtime path comparison (median of {SAMPLES}) ===");
    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "manual fixture residual harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "sequential orchestrator residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel orchestrator residual (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel2 Jacobian values+matrix: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), JACOBIAN_ITERS)
    );
    println!(
        "manual fixture sparse values harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "sequential orchestrator sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "parallel orchestrator sparse values (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), JACOBIAN_ITERS)
    );
}

#[test]
fn large_real_bvp_parallel2_and_generated_aot_match_exactly_on_same_inputs() {
    let n_steps = generated_bvp_large_fixture_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 8,
    });
    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        &generated_bvp_large_fixture_bindings::RESIDUAL_CHUNKS,
    );
    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    assert_eq!(residual_parallel2.len(), residual_generated.len());
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        let lhs = *lhs;
        let rhs = *rhs;
        assert!(
            (lhs - rhs).abs() < 1e-11,
            "large residual mismatch: {lhs} vs {rhs}"
        );
    }

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });
    let values_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        &generated_bvp_large_fixture_bindings::SPARSE_VALUE_CHUNKS,
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(values_generated.as_slice())
        .to_DMatrixType();
    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    assert_eq!(jacobian_parallel2.shape(), jacobian_generated.shape());
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-11,
                "large jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

#[test]
fn benchmark_large_real_bvp_orchestrator_runtime_paths() {
    const SAMPLES: usize = 5;
    const RESIDUAL_ITERS: usize = 180;
    const JACOBIAN_ITERS: usize = 40;

    let n_steps = generated_bvp_large_fixture_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 8,
    });
    let residual_bindings = bind_residual_fixture_chunks(
        &residual_plan,
        &generated_bvp_large_fixture_bindings::RESIDUAL_CHUNKS,
    );
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 });
    let sparse_bindings = bind_sparse_fixture_chunks(
        &sparse_plan,
        &generated_bvp_large_fixture_bindings::SPARSE_VALUE_CHUNKS,
    );
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    assert_eq!(
        residual_sequential.eval(flat_args.as_slice()),
        residual_parallel.eval(flat_args.as_slice())
    );
    assert_eq!(
        sparse_sequential.eval_values(flat_args.as_slice()),
        sparse_parallel.eval_values(flat_args.as_slice())
    );

    let mut parallel2_residual_samples = Vec::with_capacity(SAMPLES);
    let mut manual_residual_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_residual_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_residual_samples = Vec::with_capacity(SAMPLES);

    let mut parallel2_jacobian_samples = Vec::with_capacity(SAMPLES);
    let mut manual_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_sparse_samples = Vec::with_capacity(SAMPLES);

    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = jac.residiual_function.call(0.0, &input_col);
            black_box(out.get_val(0));
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = eval_residual_values_from_plan_and_chunks(
                flat_args.as_slice(),
                &residual_plan,
                &generated_bvp_large_fixture_bindings::RESIDUAL_CHUNKS,
            );
            black_box(out[0]);
        }
        manual_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = jac
                .jac_function
                .as_mut()
                .expect("parallel2 Jacobian should exist")
                .call(0.0, &input_col);
            black_box(out.shape());
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = eval_sparse_values_from_plan_and_chunks(
                flat_args.as_slice(),
                &sparse_plan,
                &generated_bvp_large_fixture_bindings::SPARSE_VALUE_CHUNKS,
            );
            black_box(out[0]);
        }
        manual_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!(
        "\n=== Larger real BVP runtime path comparison (n_steps = {}, median of {SAMPLES}) ===",
        n_steps
    );
    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "manual fixture residual harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "sequential orchestrator residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel orchestrator residual (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel2 Jacobian values+matrix: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), JACOBIAN_ITERS)
    );
    println!(
        "manual fixture sparse values harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "sequential orchestrator sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "parallel orchestrator sparse values (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), JACOBIAN_ITERS)
    );
}

#[test]
fn huge_real_bvp_parallel2_and_generated_aot_match_exactly_on_same_inputs() {
    let n_steps = generated_bvp_huge_fixture_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.00075)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp128_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 64,
    });
    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        &generated_bvp_huge_fixture_bindings::RESIDUAL_CHUNKS,
    );
    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    assert_eq!(residual_parallel2.len(), residual_generated.len());
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        let lhs = *lhs;
        let rhs = *rhs;
        assert!(
            (lhs - rhs).abs() < 1e-11,
            "huge residual mismatch: {lhs} vs {rhs}"
        );
    }

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp128_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 64 });
    let values_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        &generated_bvp_huge_fixture_bindings::SPARSE_VALUE_CHUNKS,
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(values_generated.as_slice())
        .to_DMatrixType();
    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    assert_eq!(jacobian_parallel2.shape(), jacobian_generated.shape());
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-11,
                "huge jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

#[test]
fn benchmark_huge_real_bvp_orchestrator_runtime_paths() {
    const SAMPLES: usize = 5;
    const RESIDUAL_ITERS: usize = 50;
    const JACOBIAN_ITERS: usize = 16;

    let n_steps = generated_bvp_huge_fixture_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.00075)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp128_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 64,
    });
    let residual_bindings = bind_residual_fixture_chunks(
        &residual_plan,
        &generated_bvp_huge_fixture_bindings::RESIDUAL_CHUNKS,
    );
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp128_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 64 });
    let sparse_bindings = bind_sparse_fixture_chunks(
        &sparse_plan,
        &generated_bvp_huge_fixture_bindings::SPARSE_VALUE_CHUNKS,
    );
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    assert_eq!(
        residual_sequential.eval(flat_args.as_slice()),
        residual_parallel.eval(flat_args.as_slice())
    );
    assert_eq!(
        sparse_sequential.eval_values(flat_args.as_slice()),
        sparse_parallel.eval_values(flat_args.as_slice())
    );

    let mut parallel2_residual_samples = Vec::with_capacity(SAMPLES);
    let mut manual_residual_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_residual_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_residual_samples = Vec::with_capacity(SAMPLES);

    let mut parallel2_jacobian_samples = Vec::with_capacity(SAMPLES);
    let mut manual_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut sequential_sparse_samples = Vec::with_capacity(SAMPLES);
    let mut parallel_sparse_samples = Vec::with_capacity(SAMPLES);

    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = jac.residiual_function.call(0.0, &input_col);
            black_box(out.get_val(0));
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = eval_residual_values_from_plan_and_chunks(
                flat_args.as_slice(),
                &residual_plan,
                &generated_bvp_huge_fixture_bindings::RESIDUAL_CHUNKS,
            );
            black_box(out[0]);
        }
        manual_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = jac
                .jac_function
                .as_mut()
                .expect("parallel2 Jacobian should exist")
                .call(0.0, &input_col);
            black_box(out.shape());
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = eval_sparse_values_from_plan_and_chunks(
                flat_args.as_slice(),
                &sparse_plan,
                &generated_bvp_huge_fixture_bindings::SPARSE_VALUE_CHUNKS,
            );
            black_box(out[0]);
        }
        manual_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!(
        "\n=== Huge real BVP runtime path comparison (n_steps = {}, median of {SAMPLES}) ===",
        n_steps
    );
    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "manual fixture residual harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "sequential orchestrator residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel orchestrator residual (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), RESIDUAL_ITERS)
    );
    println!(
        "parallel2 Jacobian values+matrix: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), JACOBIAN_ITERS)
    );
    println!(
        "manual fixture sparse values harness: {:.2} ns/call",
        per_iter_ns(median_duration(manual_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "sequential orchestrator sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), JACOBIAN_ITERS)
    );
    println!(
        "parallel orchestrator sparse values (forced parallel): {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), JACOBIAN_ITERS)
    );
}

#[test]
fn benchmark_large_real_bvp_residual_chunk_sizes() {
    const SAMPLES: usize = 5;
    let n_steps = generated_bvp_large_fixture_bindings::N_STEPS;
    let jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };

    let strategies = [
        ResidualChunkingStrategy::Whole,
        ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 2 },
        ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
        ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: 1,
        },
        ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: 2,
        },
        ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: 4,
        },
        ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: 8,
        },
        ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: 16,
        },
    ];

    let recommended =
        recommended_residual_chunking_for_parallelism(jac.vector_of_functions.len(), 4);

    println!(
        "\n=== Larger real BVP residual chunk-size sweep (n_steps = {}, outputs = {}) ===",
        n_steps,
        jac.vector_of_functions.len()
    );
    println!(
        "recommended residual strategy (4 chunks/worker): {:?}",
        recommended
    );

    for strategy in strategies {
        let mut planning_samples = Vec::with_capacity(SAMPLES);
        let mut assembly_samples = Vec::with_capacity(SAMPLES);
        let mut emit_samples = Vec::with_capacity(SAMPLES);

        let runtime_plan = residual_task.runtime_plan(strategy);
        let chunk_count = runtime_plan.chunks.len();
        let max_chunk_outputs = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.residuals.len())
            .max()
            .unwrap_or(0);
        let min_chunk_outputs = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.residuals.len())
            .min()
            .unwrap_or(0);

        for _ in 0..SAMPLES {
            let start = Instant::now();
            let plan = residual_task.runtime_plan(strategy);
            black_box(plan.chunks.len());
            planning_samples.push(start.elapsed());

            let start = Instant::now();
            let mut module = CodegenModule::new("generated_large_residual_strategy_probe");
            for chunk in &runtime_plan.chunks {
                module.push_residual_block_plan(&chunk.plan);
            }
            assembly_samples.push(start.elapsed());

            let start = Instant::now();
            let source = module.emit_source();
            black_box(source.len());
            emit_samples.push(start.elapsed());
        }

        let mut module = CodegenModule::new("generated_large_residual_strategy_probe");
        for chunk in &runtime_plan.chunks {
            module.push_residual_block_plan(&chunk.plan);
        }
        let source_len = module.emit_source().len();

        println!(
            "{strategy:?}: chunks={}, outputs/chunk min..max={}..{}, planning={:?}, assembly={:?}, emit={:?}, source_len={}",
            chunk_count,
            min_chunk_outputs,
            max_chunk_outputs,
            median_duration(planning_samples),
            median_duration(assembly_samples),
            median_duration(emit_samples),
            source_len,
        );
    }
}

#[test]
fn benchmark_large_real_bvp_sparse_chunk_sizes() {
    const SAMPLES: usize = 5;
    let n_steps = generated_bvp_large_fixture_bindings::N_STEPS;
    let jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);

    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };

    let strategies = [
        SparseChunkingStrategy::Whole,
        SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 2 },
        SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 4 },
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 1 },
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 2 },
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 4 },
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 8 },
        SparseChunkingStrategy::ByRowCount { rows_per_chunk: 16 },
    ];

    let recommended = recommended_row_chunking_for_parallelism(jac.vector_of_functions.len(), 4);

    println!(
        "\n=== Larger real BVP sparse chunk-size sweep (n_steps = {}, rows = {}, nnz = {}) ===",
        n_steps,
        jac.vector_of_functions.len(),
        sparse_entries.len()
    );
    println!(
        "recommended sparse strategy (4 chunks/worker): {:?}",
        recommended
    );

    for strategy in strategies {
        let mut planning_samples = Vec::with_capacity(SAMPLES);
        let mut assembly_samples = Vec::with_capacity(SAMPLES);
        let mut emit_samples = Vec::with_capacity(SAMPLES);

        let runtime_plan = sparse_task.runtime_plan(strategy);
        let chunk_count = runtime_plan.chunks.len();
        let max_chunk_nnz = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.entries.len())
            .max()
            .unwrap_or(0);
        let min_chunk_nnz = runtime_plan
            .chunks
            .iter()
            .map(|chunk| chunk.entries.len())
            .min()
            .unwrap_or(0);

        for _ in 0..SAMPLES {
            let start = Instant::now();
            let plan = sparse_task.runtime_plan(strategy);
            black_box(plan.chunks.len());
            planning_samples.push(start.elapsed());

            let start = Instant::now();
            let mut module = CodegenModule::new("generated_large_sparse_strategy_probe");
            for chunk in &runtime_plan.chunks {
                module.push_sparse_values_plan(&chunk.plan);
            }
            assembly_samples.push(start.elapsed());

            let start = Instant::now();
            let source = module.emit_source();
            black_box(source.len());
            emit_samples.push(start.elapsed());
        }

        let mut module = CodegenModule::new("generated_large_sparse_strategy_probe");
        for chunk in &runtime_plan.chunks {
            module.push_sparse_values_plan(&chunk.plan);
        }
        let source_len = module.emit_source().len();

        println!(
            "{strategy:?}: chunks={}, nnz/chunk min..max={}..{}, planning={:?}, assembly={:?}, emit={:?}, source_len={}",
            chunk_count,
            min_chunk_nnz,
            max_chunk_nnz,
            median_duration(planning_samples),
            median_duration(assembly_samples),
            median_duration(emit_samples),
            source_len,
        );
    }
}

#[test]
fn large_real_bvp_compiled_chunk_variants_match_parallel2() {
    let n_steps = generated_bvp_large_chunk_variant_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();

    for outputs_per_chunk in [4usize, 8, 16] {
        let residual_task = ResidualTask {
            fn_name: "fixture_variant_residual",
            residuals: &jac.vector_of_functions,
            variables: &variable_names,
            params: None,
        };
        let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: outputs_per_chunk,
        });
        let residual_generated = eval_residual_values_from_plan_and_chunks(
            flat_args.as_slice(),
            &residual_plan,
            residual_variant_chunk_fns(outputs_per_chunk),
        );
        for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
            assert!(
                (*lhs - *rhs).abs() < 1e-11,
                "residual variant O{outputs_per_chunk} mismatch: {lhs} vs {rhs}"
            );
        }
    }

    for rows_per_chunk in [4usize, 8, 16] {
        let sparse_task = SparseJacobianTask {
            fn_name: "fixture_variant_sparse",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };
        let sparse_plan =
            sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk });
        let values_generated = eval_sparse_values_from_plan_and_chunks(
            flat_args.as_slice(),
            &sparse_plan,
            sparse_variant_chunk_fns(rows_per_chunk),
        );
        let jacobian_generated = sparse_plan
            .assemble_sparse_col_mat(values_generated.as_slice())
            .to_DMatrixType();
        for i in 0..jacobian_parallel2.nrows() {
            for j in 0..jacobian_parallel2.ncols() {
                let lhs = jacobian_parallel2[(i, j)];
                let rhs = jacobian_generated[(i, j)];
                assert!(
                    (lhs - rhs).abs() < 1e-11,
                    "sparse variant R{rows_per_chunk} mismatch at ({i}, {j}): {lhs} vs {rhs}"
                );
            }
        }
    }
}

#[test]
fn benchmark_large_real_bvp_compiled_chunk_variants_runtime() {
    const SAMPLES: usize = 5;
    const RESIDUAL_ITERS: usize = 180;
    const JACOBIAN_ITERS: usize = 40;

    let n_steps = generated_bvp_large_chunk_variant_bindings::N_STEPS;
    let jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);

    println!(
        "\n=== Larger real BVP compiled chunk-variant runtime sweep (n_steps = {}) ===",
        n_steps
    );

    for outputs_per_chunk in [4usize, 8, 16] {
        let residual_task = ResidualTask {
            fn_name: "fixture_variant_residual",
            residuals: &jac.vector_of_functions,
            variables: &variable_names,
            params: None,
        };
        let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
            max_outputs_per_chunk: outputs_per_chunk,
        });
        let residual_bindings = bind_residual_fixture_chunks(
            &residual_plan,
            residual_variant_chunk_fns(outputs_per_chunk),
        );
        let residual_sequential =
            SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
        let residual_parallel = ParallelResidualExecutor::with_config(
            &residual_plan,
            residual_bindings,
            benchmark_parallel_config(),
        );

        let mut sequential_samples = Vec::with_capacity(SAMPLES);
        let mut parallel_samples = Vec::with_capacity(SAMPLES);
        for _ in 0..SAMPLES {
            let start = Instant::now();
            for _ in 0..RESIDUAL_ITERS {
                let out = residual_sequential.eval(flat_args.as_slice());
                black_box(out[0]);
            }
            sequential_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..RESIDUAL_ITERS {
                let out = residual_parallel.eval(flat_args.as_slice());
                black_box(out[0]);
            }
            parallel_samples.push(start.elapsed());
        }

        println!(
            "Residual O{}: outputs={}, chunks={}, jobs={}, fallback={}, sequential={:.2} ns/call, parallel={:.2} ns/call",
            outputs_per_chunk,
            residual_plan.output_len,
            residual_plan.chunks.len(),
            residual_parallel.job_count(),
            residual_parallel.uses_sequential_fallback(),
            per_iter_ns(median_duration(sequential_samples), RESIDUAL_ITERS),
            per_iter_ns(median_duration(parallel_samples), RESIDUAL_ITERS),
        );
    }

    for rows_per_chunk in [4usize, 8, 16] {
        let sparse_task = SparseJacobianTask {
            fn_name: "fixture_variant_sparse",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };
        let sparse_plan =
            sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk });
        let sparse_bindings =
            bind_sparse_fixture_chunks(&sparse_plan, sparse_variant_chunk_fns(rows_per_chunk));
        let sparse_sequential =
            SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
        let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
            &sparse_plan,
            sparse_bindings,
            benchmark_parallel_config(),
        );

        let mut sequential_samples = Vec::with_capacity(SAMPLES);
        let mut parallel_samples = Vec::with_capacity(SAMPLES);
        for _ in 0..SAMPLES {
            let start = Instant::now();
            for _ in 0..JACOBIAN_ITERS {
                let out = sparse_sequential.eval_values(flat_args.as_slice());
                black_box(out[0]);
            }
            sequential_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..JACOBIAN_ITERS {
                let out = sparse_parallel.eval_values(flat_args.as_slice());
                black_box(out[0]);
            }
            parallel_samples.push(start.elapsed());
        }

        println!(
            "Sparse R{}: nnz={}, chunks={}, jobs={}, fallback={}, sequential={:.2} ns/call, parallel={:.2} ns/call",
            rows_per_chunk,
            sparse_plan.nnz(),
            sparse_plan.chunks.len(),
            sparse_parallel.job_count(),
            sparse_parallel.uses_sequential_fallback(),
            per_iter_ns(median_duration(sequential_samples), JACOBIAN_ITERS),
            per_iter_ns(median_duration(parallel_samples), JACOBIAN_ITERS),
        );
    }
}

#[test]
fn large_real_bvp_superblock_variants_match_parallel2() {
    let n_steps = generated_bvp_large_superblock_variant_bindings::N_STEPS;
    let mut jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual_t2",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task
        .runtime_plan(ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        residual_superblock_chunk_fns(2),
    );
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        assert!(
            (*lhs - *rhs).abs() < 1e-11,
            "residual superblock T2 mismatch: {lhs} vs {rhs}"
        );
    }

    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_t2",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let values_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        sparse_superblock_chunk_fns(2),
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(values_generated.as_slice())
        .to_DMatrixType();
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-11,
                "sparse superblock T2 mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

#[test]
fn benchmark_large_real_bvp_superblock_variants_runtime() {
    const SAMPLES: usize = 5;
    const RESIDUAL_ITERS: usize = 180;
    const JACOBIAN_ITERS: usize = 40;

    let n_steps = generated_bvp_large_superblock_variant_bindings::N_STEPS;
    let jac = build_real_bvp_damp1_case(n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.15 + index as f64 * 0.0025)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);

    println!(
        "\n=== Larger real BVP superblock runtime sweep (n_steps = {}) ===",
        n_steps
    );

    let residual_task = ResidualTask {
        fn_name: "fixture_bvp32_residual_t2",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task
        .runtime_plan(ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let residual_bindings =
        bind_residual_fixture_chunks(&residual_plan, residual_superblock_chunk_fns(2));
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let mut residual_seq_samples = Vec::with_capacity(SAMPLES);
    let mut residual_par_samples = Vec::with_capacity(SAMPLES);
    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        residual_seq_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..RESIDUAL_ITERS {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        residual_par_samples.push(start.elapsed());
    }

    println!(
        "Residual T2: outputs={}, chunks={}, jobs={}, fallback={}, sequential={:.2} ns/call, parallel={:.2} ns/call",
        residual_plan.output_len,
        residual_plan.chunks.len(),
        residual_parallel.job_count(),
        residual_parallel.uses_sequential_fallback(),
        per_iter_ns(median_duration(residual_seq_samples), RESIDUAL_ITERS),
        per_iter_ns(median_duration(residual_par_samples), RESIDUAL_ITERS),
    );

    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_bvp32_sparse_t2",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 2 });
    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, sparse_superblock_chunk_fns(2));
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    let mut sparse_seq_samples = Vec::with_capacity(SAMPLES);
    let mut sparse_par_samples = Vec::with_capacity(SAMPLES);
    for _ in 0..SAMPLES {
        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sparse_seq_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..JACOBIAN_ITERS {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sparse_par_samples.push(start.elapsed());
    }

    println!(
        "Sparse T2: nnz={}, chunks={}, jobs={}, fallback={}, sequential={:.2} ns/call, parallel={:.2} ns/call",
        sparse_plan.nnz(),
        sparse_plan.chunks.len(),
        sparse_parallel.job_count(),
        sparse_parallel.uses_sequential_fallback(),
        per_iter_ns(median_duration(sparse_seq_samples), JACOBIAN_ITERS),
        per_iter_ns(median_duration(sparse_par_samples), JACOBIAN_ITERS),
    );
}

pub(crate) fn run_benchmark_stress_multifield_bvp_pipeline() {
    let samples = stress_samples(7);
    let residual_iters = stress_residual_iters(60);
    let jacobian_iters = stress_jacobian_iters(20);
    let field_count = env_usize("RST_STRESS_FIELDS", 8);
    let n_steps = env_usize("RST_STRESS_STEPS", 96);
    let estimated_unknowns = field_count * n_steps;
    let estimated_dense_cells = estimated_unknowns * estimated_unknowns;

    println!(
        "\n=== Stress multifield BVP pipeline ===\nfields={}, steps={}, estimated_unknowns={}, estimated_dense_jacobian_cells={}",
        field_count, n_steps, estimated_unknowns, estimated_dense_cells
    );
    println!(
        "samples={}, residual_iters={}, jacobian_iters={}",
        samples, residual_iters, jacobian_iters
    );

    let start = Instant::now();
    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let symbolic_build = start.elapsed();

    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let start = Instant::now();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    let parallel2_residual_compile = start.elapsed();

    let start = Instant::now();
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());
    let parallel2_jacobian_compile = start.elapsed();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "stress_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "stress_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };

    let residual_strategy =
        recommended_residual_chunking_for_parallelism(jac.vector_of_functions.len(), 4);
    let sparse_strategy =
        recommended_row_chunking_for_parallelism(jac.vector_of_functions.len(), 4);

    let mut extract_samples = Vec::with_capacity(samples);
    let mut residual_plan_samples = Vec::with_capacity(samples);
    let mut sparse_plan_samples = Vec::with_capacity(samples);
    let mut assembly_samples = Vec::with_capacity(samples);
    let mut emit_samples = Vec::with_capacity(samples);

    for _ in 0..samples {
        let start = Instant::now();
        let sparse_entries_repeat_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_entries_repeat = borrowed_sparse_entries(&sparse_entries_repeat_owned);
        black_box(sparse_entries_repeat.len());
        extract_samples.push(start.elapsed());

        let start = Instant::now();
        let residual_plan = residual_task.runtime_plan(residual_strategy);
        black_box(residual_plan.chunks.len());
        residual_plan_samples.push(start.elapsed());

        let start = Instant::now();
        let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
        black_box(sparse_plan.chunks.len());
        sparse_plan_samples.push(start.elapsed());

        let residual_plan = residual_task.runtime_plan(residual_strategy);
        let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
        let start = Instant::now();
        let module =
            residual_and_sparse_module("generated_stress_probe", &residual_plan, &sparse_plan);
        assembly_samples.push(start.elapsed());

        let start = Instant::now();
        let source = module.emit_source();
        black_box(source.len());
        emit_samples.push(start.elapsed());
    }

    let residual_plan = residual_task.runtime_plan(residual_strategy);
    let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
    let module = residual_and_sparse_module("generated_stress_probe", &residual_plan, &sparse_plan);
    let source_len = module.emit_source().len();

    let mut residual_runtime_samples = Vec::with_capacity(samples);
    let mut jacobian_runtime_samples = Vec::with_capacity(samples);
    for _ in 0..samples {
        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = jac
                .residiual_function
                .call(0.0, &input_col)
                .to_DVectorType();
            black_box(out[0]);
        }
        residual_runtime_samples.push(start.elapsed());

        let jacobian_fun = jac
            .jac_function
            .as_mut()
            .expect("parallel2 Jacobian should be compiled");
        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let mat = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
            black_box(mat[(0, 0)]);
        }
        jacobian_runtime_samples.push(start.elapsed());
    }

    println!("symbolic build: {:?}", symbolic_build);
    println!(
        "parallel2 residual compile: {:?}",
        parallel2_residual_compile
    );
    println!(
        "parallel2 Jacobian compile: {:?}",
        parallel2_jacobian_compile
    );
    println!(
        "AOT sparse extraction: {:?}",
        median_duration(extract_samples)
    );
    println!(
        "AOT residual planning ({residual_strategy:?}): {:?}",
        median_duration(residual_plan_samples)
    );
    println!(
        "AOT sparse planning ({sparse_strategy:?}): {:?}",
        median_duration(sparse_plan_samples)
    );
    println!(
        "AOT module assembly: {:?}",
        median_duration(assembly_samples)
    );
    println!(
        "AOT emit_source: {:?} (source_len={})",
        median_duration(emit_samples),
        source_len
    );
    println!(
        "parallel2 residual runtime: {:.2} ns/call",
        per_iter_ns(median_duration(residual_runtime_samples), residual_iters)
    );
    println!(
        "parallel2 Jacobian runtime: {:.2} ns/call",
        per_iter_ns(median_duration(jacobian_runtime_samples), jacobian_iters)
    );
    println!(
        "stress plans: residual_chunks={}, sparse_chunks={}, nnz={}",
        residual_plan.chunks.len(),
        sparse_plan.chunks.len(),
        sparse_plan.nnz()
    );
}

/// Scaling sweep for larger multifield systems.
///
/// This benchmark is intentionally different from the checked-in compiled
/// fixture benchmarks above:
/// - those benchmarks measure emitted-and-compiled AOT Rust functions for one
///   fixed generated case;
/// - this benchmark measures the AOT pipeline at the lowered IR level so we
///   can compare sequential vs parallel AOT execution on arbitrary larger
///   stress sizes without generating a new fixture module for every size.
///
/// In other words:
/// - "compiled AOT fixture benchmark" answers "how fast is one concrete
///   generated Rust module?";
/// - "IR-level AOT scaling benchmark" answers "how does sequential vs
///   parallel AOT behave as the symbolic problem gets much larger?".
pub(crate) fn run_benchmark_stress_multifield_scaling_sweep() {
    let samples = stress_samples(5);
    let residual_iters = stress_residual_iters(24);
    let jacobian_iters = stress_jacobian_iters(8);
    let field_counts = stress_env_list_usize("RST_STRESS_FIELD_LIST", &[8, 16, 24, 32]);
    let n_steps = env_usize("RST_STRESS_SWEEP_STEPS", 128);

    println!(
        "\n=== Stress multifield scaling sweep ===\nfields={:?}, steps={}, samples={}, residual_iters={}, jacobian_iters={}",
        field_counts, n_steps, samples, residual_iters, jacobian_iters
    );

    for field_count in field_counts {
        let estimated_unknowns = field_count * n_steps;
        let estimated_dense_cells = estimated_unknowns * estimated_unknowns;
        println!(
            "\n--- field_count={}, steps={}, estimated_unknowns={}, estimated_dense_jacobian_cells={} ---",
            field_count, n_steps, estimated_unknowns, estimated_dense_cells
        );

        let start = Instant::now();
        let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
        let symbolic_build = start.elapsed();

        let variable_names_owned = jac.variable_string.clone();
        let variable_names: Vec<&str> = variable_names_owned
            .iter()
            .map(|name| name.as_str())
            .collect();

        let variable_values: Vec<f64> = (0..jac.variable_string.len())
            .map(|index| 0.05 + index as f64 * 0.0005)
            .collect();
        let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
        let flat_args = runtime_args.flatten();
        let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

        let start = Instant::now();
        jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
        let parallel2_residual_compile = start.elapsed();

        let start = Instant::now();
        jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());
        let parallel2_jacobian_compile = start.elapsed();

        let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
        let residual_task = ResidualTask {
            fn_name: "stress_scaling_residual",
            residuals: &jac.vector_of_functions,
            variables: &variable_names,
            params: None,
        };
        let sparse_task = SparseJacobianTask {
            fn_name: "stress_scaling_sparse_values",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };

        let residual_strategy =
            recommended_residual_chunking_for_parallelism(jac.vector_of_functions.len(), 4);
        let sparse_strategy =
            recommended_row_chunking_for_parallelism(jac.vector_of_functions.len(), 4);
        let residual_plan = residual_task.runtime_plan(residual_strategy);
        let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
        let residual_ir_blocks = residual_ir_blocks(&residual_plan);
        let sparse_ir_blocks = sparse_ir_blocks(&sparse_plan);

        let residual_ir_seq = eval_ir_blocks_sequential(
            flat_args.as_slice(),
            &residual_ir_blocks,
            residual_plan.output_len,
        );
        let residual_ir_par = eval_ir_blocks_parallel(
            flat_args.as_slice(),
            &residual_ir_blocks,
            residual_plan.output_len,
        );
        assert_eq!(residual_ir_seq.len(), residual_ir_par.len());
        for (lhs, rhs) in residual_ir_seq.iter().zip(residual_ir_par.iter()) {
            assert!(
                (*lhs - *rhs).abs() < 1e-10,
                "stress residual IR mismatch: {lhs} vs {rhs}"
            );
        }

        let sparse_ir_seq =
            eval_ir_blocks_sequential(flat_args.as_slice(), &sparse_ir_blocks, sparse_plan.nnz());
        let sparse_ir_par =
            eval_ir_blocks_parallel(flat_args.as_slice(), &sparse_ir_blocks, sparse_plan.nnz());
        assert_eq!(sparse_ir_seq.len(), sparse_ir_par.len());
        for (lhs, rhs) in sparse_ir_seq.iter().zip(sparse_ir_par.iter()) {
            assert!(
                (*lhs - *rhs).abs() < 1e-10,
                "stress sparse IR mismatch: {lhs} vs {rhs}"
            );
        }

        let mut residual_plan_samples = Vec::with_capacity(samples);
        let mut sparse_plan_samples = Vec::with_capacity(samples);
        let mut assembly_samples = Vec::with_capacity(samples);
        let mut emit_samples = Vec::with_capacity(samples);
        let mut residual_runtime_samples = Vec::with_capacity(samples);
        let mut jacobian_runtime_samples = Vec::with_capacity(samples);
        let mut residual_aot_seq_samples = Vec::with_capacity(samples);
        let mut residual_aot_par_samples = Vec::with_capacity(samples);
        let mut sparse_aot_seq_samples = Vec::with_capacity(samples);
        let mut sparse_aot_par_samples = Vec::with_capacity(samples);

        for _ in 0..samples {
            let start = Instant::now();
            let residual_plan = residual_task.runtime_plan(residual_strategy);
            black_box(residual_plan.chunks.len());
            residual_plan_samples.push(start.elapsed());

            let start = Instant::now();
            let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
            black_box(sparse_plan.chunks.len());
            sparse_plan_samples.push(start.elapsed());

            let residual_plan = residual_task.runtime_plan(residual_strategy);
            let sparse_plan = sparse_task.runtime_plan(sparse_strategy);
            let start = Instant::now();
            let module = residual_and_sparse_module(
                "generated_stress_scaling_probe",
                &residual_plan,
                &sparse_plan,
            );
            assembly_samples.push(start.elapsed());

            let start = Instant::now();
            let source = module.emit_source();
            black_box(source.len());
            emit_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..residual_iters {
                let out = jac
                    .residiual_function
                    .call(0.0, &input_col)
                    .to_DVectorType();
                black_box(out[0]);
            }
            residual_runtime_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..residual_iters {
                let out = eval_ir_blocks_sequential(
                    flat_args.as_slice(),
                    &residual_ir_blocks,
                    residual_plan.output_len,
                );
                black_box(out[0]);
            }
            residual_aot_seq_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..residual_iters {
                let out = eval_ir_blocks_parallel(
                    flat_args.as_slice(),
                    &residual_ir_blocks,
                    residual_plan.output_len,
                );
                black_box(out[0]);
            }
            residual_aot_par_samples.push(start.elapsed());

            let jacobian_fun = jac
                .jac_function
                .as_mut()
                .expect("parallel2 Jacobian should be compiled");
            let start = Instant::now();
            for _ in 0..jacobian_iters {
                let mat = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
                black_box(mat[(0, 0)]);
            }
            jacobian_runtime_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..jacobian_iters {
                let out = eval_ir_blocks_sequential(
                    flat_args.as_slice(),
                    &sparse_ir_blocks,
                    sparse_plan.nnz(),
                );
                black_box(out[0]);
            }
            sparse_aot_seq_samples.push(start.elapsed());

            let start = Instant::now();
            for _ in 0..jacobian_iters {
                let out = eval_ir_blocks_parallel(
                    flat_args.as_slice(),
                    &sparse_ir_blocks,
                    sparse_plan.nnz(),
                );
                black_box(out[0]);
            }
            sparse_aot_par_samples.push(start.elapsed());
        }
        let module = residual_and_sparse_module(
            "generated_stress_scaling_probe",
            &residual_plan,
            &sparse_plan,
        );
        let source_len = module.emit_source().len();

        println!("symbolic build: {:?}", symbolic_build);
        println!(
            "parallel2 residual compile: {:?}",
            parallel2_residual_compile
        );
        println!(
            "parallel2 Jacobian compile: {:?}",
            parallel2_jacobian_compile
        );
        println!(
            "AOT residual planning ({residual_strategy:?}): {:?}",
            median_duration(residual_plan_samples)
        );
        println!(
            "AOT sparse planning ({sparse_strategy:?}): {:?}",
            median_duration(sparse_plan_samples)
        );
        println!(
            "AOT module assembly: {:?}",
            median_duration(assembly_samples)
        );
        println!(
            "AOT emit_source: {:?} (source_len={})",
            median_duration(emit_samples),
            source_len
        );
        println!(
            "parallel2 residual runtime: {:.2} ns/call",
            per_iter_ns(median_duration(residual_runtime_samples), residual_iters)
        );
        println!(
            "sequential AOT IR residual runtime: {:.2} ns/call",
            per_iter_ns(median_duration(residual_aot_seq_samples), residual_iters)
        );
        println!(
            "parallel AOT IR residual runtime: {:.2} ns/call",
            per_iter_ns(median_duration(residual_aot_par_samples), residual_iters)
        );
        println!(
            "parallel2 Jacobian runtime: {:.2} ns/call",
            per_iter_ns(median_duration(jacobian_runtime_samples), jacobian_iters)
        );
        println!(
            "sequential AOT IR sparse values runtime: {:.2} ns/call",
            per_iter_ns(median_duration(sparse_aot_seq_samples), jacobian_iters)
        );
        println!(
            "parallel AOT IR sparse values runtime: {:.2} ns/call",
            per_iter_ns(median_duration(sparse_aot_par_samples), jacobian_iters)
        );
        println!(
            "stress scaling plans: residual_chunks={}, sparse_chunks={}, nnz={}",
            residual_plan.chunks.len(),
            sparse_plan.chunks.len(),
            sparse_plan.nnz()
        );
    }
}

pub(crate) fn run_benchmark_stress_multifield_parallel_crossover() {
    let samples = stress_samples(7);
    let residual_iters = stress_residual_iters(24);
    let jacobian_iters = stress_jacobian_iters(8);
    let field_counts = stress_env_list_usize("RST_STRESS_FIELD_LIST", &[8, 16, 24, 32, 48, 64]);
    let n_steps = env_usize("RST_STRESS_SWEEP_STEPS", 128);

    println!(
        "\n=== Stress multifield AOT IR crossover search ===\nfields={:?}, steps={}, samples={}, residual_iters={}, jacobian_iters={}",
        field_counts, n_steps, samples, residual_iters, jacobian_iters
    );

    let mut residual_crossover = None;
    let mut sparse_crossover = None;

    for field_count in field_counts {
        let jac = build_multifield_stress_bvp_case(field_count, n_steps);
        let variable_names_owned = jac.variable_string.clone();
        let variable_names: Vec<&str> = variable_names_owned
            .iter()
            .map(|name| name.as_str())
            .collect();

        let variable_values: Vec<f64> = (0..jac.variable_string.len())
            .map(|index| 0.05 + index as f64 * 0.0005)
            .collect();
        let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
        let flat_args = runtime_args.flatten();

        let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
        let residual_task = ResidualTask {
            fn_name: "stress_crossover_residual",
            residuals: &jac.vector_of_functions,
            variables: &variable_names,
            params: None,
        };
        let sparse_task = SparseJacobianTask {
            fn_name: "stress_crossover_sparse_values",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };

        let residual_plan = residual_task.runtime_plan(
            recommended_residual_chunking_for_parallelism(jac.vector_of_functions.len(), 4),
        );
        let sparse_plan = sparse_task.runtime_plan(recommended_row_chunking_for_parallelism(
            jac.vector_of_functions.len(),
            4,
        ));

        let residual_ir_blocks = residual_ir_blocks(&residual_plan);
        let sparse_ir_blocks = sparse_ir_blocks(&sparse_plan);

        let residual_seq_probe = eval_ir_blocks_sequential(
            flat_args.as_slice(),
            &residual_ir_blocks,
            residual_plan.output_len,
        );
        let residual_par_probe = eval_ir_blocks_parallel(
            flat_args.as_slice(),
            &residual_ir_blocks,
            residual_plan.output_len,
        );
        for (lhs, rhs) in residual_seq_probe.iter().zip(residual_par_probe.iter()) {
            assert!(
                (*lhs - *rhs).abs() < 1e-10,
                "residual crossover mismatch: {lhs} vs {rhs}"
            );
        }

        let sparse_seq_probe =
            eval_ir_blocks_sequential(flat_args.as_slice(), &sparse_ir_blocks, sparse_plan.nnz());
        let sparse_par_probe =
            eval_ir_blocks_parallel(flat_args.as_slice(), &sparse_ir_blocks, sparse_plan.nnz());
        for (lhs, rhs) in sparse_seq_probe.iter().zip(sparse_par_probe.iter()) {
            assert!(
                (*lhs - *rhs).abs() < 1e-10,
                "sparse crossover mismatch: {lhs} vs {rhs}"
            );
        }

        let residual_seq_samples: Vec<Duration> = (0..samples)
            .map(|_| {
                let start = Instant::now();
                for _ in 0..residual_iters {
                    let out = eval_ir_blocks_sequential(
                        flat_args.as_slice(),
                        &residual_ir_blocks,
                        residual_plan.output_len,
                    );
                    black_box(out[0]);
                }
                start.elapsed()
            })
            .collect();
        let residual_par_samples: Vec<Duration> = (0..samples)
            .map(|_| {
                let start = Instant::now();
                for _ in 0..residual_iters {
                    let out = eval_ir_blocks_parallel(
                        flat_args.as_slice(),
                        &residual_ir_blocks,
                        residual_plan.output_len,
                    );
                    black_box(out[0]);
                }
                start.elapsed()
            })
            .collect();

        let sparse_seq_samples: Vec<Duration> = (0..samples)
            .map(|_| {
                let start = Instant::now();
                for _ in 0..jacobian_iters {
                    let out = eval_ir_blocks_sequential(
                        flat_args.as_slice(),
                        &sparse_ir_blocks,
                        sparse_plan.nnz(),
                    );
                    black_box(out[0]);
                }
                start.elapsed()
            })
            .collect();
        let sparse_par_samples: Vec<Duration> = (0..samples)
            .map(|_| {
                let start = Instant::now();
                for _ in 0..jacobian_iters {
                    let out = eval_ir_blocks_parallel(
                        flat_args.as_slice(),
                        &sparse_ir_blocks,
                        sparse_plan.nnz(),
                    );
                    black_box(out[0]);
                }
                start.elapsed()
            })
            .collect();

        let residual_seq_ns = per_iter_ns(median_duration(residual_seq_samples), residual_iters);
        let residual_par_ns = per_iter_ns(median_duration(residual_par_samples), residual_iters);
        let sparse_seq_ns = per_iter_ns(median_duration(sparse_seq_samples), jacobian_iters);
        let sparse_par_ns = per_iter_ns(median_duration(sparse_par_samples), jacobian_iters);

        if residual_crossover.is_none() && residual_par_ns < residual_seq_ns {
            residual_crossover = Some(field_count);
        }
        if sparse_crossover.is_none() && sparse_par_ns < sparse_seq_ns {
            sparse_crossover = Some(field_count);
        }

        println!(
            "field_count={}, unknowns={}, residual seq/par = {:.2}/{:.2} ns, sparse seq/par = {:.2}/{:.2} ns",
            field_count,
            field_count * n_steps,
            residual_seq_ns,
            residual_par_ns,
            sparse_seq_ns,
            sparse_par_ns
        );
    }

    println!(
        "residual crossover field_count: {}",
        residual_crossover
            .map(|value| value.to_string())
            .unwrap_or_else(|| "not reached".to_string())
    );
    println!(
        "sparse crossover field_count: {}",
        sparse_crossover
            .map(|value| value.to_string())
            .unwrap_or_else(|| "not reached".to_string())
    );
}

pub(crate) fn run_benchmark_compiled_stress_multifield_aot_runtime() {
    let samples = stress_samples(9);
    let residual_iters = stress_residual_iters(120);
    let jacobian_iters = stress_jacobian_iters(48);
    let field_count = generated_stress_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_bvp_fixture_bindings::N_STEPS;

    println!(
        "\n=== Compiled stress multifield AOT runtime ===\nfields={}, steps={}",
        field_count, n_steps
    );
    println!(
        "samples={}, residual_iters={}, jacobian_iters={}",
        samples, residual_iters, jacobian_iters
    );

    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "fixture_stress_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 48,
    });
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 48 });

    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        stress_residual_chunk_fns(),
    );
    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        assert!(
            (*lhs - *rhs).abs() < 1e-10,
            "stress residual mismatch: {lhs} vs {rhs}"
        );
    }

    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    let sparse_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        stress_sparse_chunk_fns(),
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(sparse_generated.as_slice())
        .to_DMatrixType();
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-10,
                "stress Jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
    let residual_bindings =
        bind_residual_fixture_chunks(&residual_plan, stress_residual_chunk_fns());
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, stress_sparse_chunk_fns());
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    let mut parallel2_residual_samples = Vec::with_capacity(samples);
    let mut sequential_residual_samples = Vec::with_capacity(samples);
    let mut parallel_residual_samples = Vec::with_capacity(samples);

    let mut parallel2_jacobian_samples = Vec::with_capacity(samples);
    let mut sequential_sparse_samples = Vec::with_capacity(samples);
    let mut parallel_sparse_samples = Vec::with_capacity(samples);

    for _ in 0..samples {
        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = jac
                .residiual_function
                .call(0.0, &input_col)
                .to_DVectorType();
            black_box(out[0]);
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let jacobian_fun = jac
            .jac_function
            .as_mut()
            .expect("parallel2 Jacobian should stay compiled");
        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
            black_box(out[(0, 0)]);
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), residual_iters)
    );
    println!(
        "sequential AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), residual_iters)
    );
    println!(
        "parallel AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), residual_iters)
    );
    println!(
        "parallel2 Jacobian: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), jacobian_iters)
    );
    println!(
        "sequential AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), jacobian_iters)
    );
    println!(
        "parallel AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), jacobian_iters)
    );
    println!(
        "compiled stress plans: residual_chunks={}, sparse_chunks={}, jobs_residual={}, jobs_sparse={}",
        residual_plan.chunks.len(),
        sparse_plan.chunks.len(),
        residual_parallel.job_count(),
        sparse_parallel.job_count()
    );
}

#[test]
fn compiled_large_stress_fixture_matches_parallel2() {
    let field_count = generated_stress_large_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_large_bvp_fixture_bindings::N_STEPS;

    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "fixture_stress_large_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_large_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 96,
    });
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 96 });

    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        stress_large_residual_chunk_fns(),
    );
    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        assert!(
            (*lhs - *rhs).abs() < 1e-10,
            "large stress residual mismatch: {lhs} vs {rhs}"
        );
    }

    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    let sparse_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        stress_large_sparse_chunk_fns(),
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(sparse_generated.as_slice())
        .to_DMatrixType();
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-10,
                "large stress Jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

pub(crate) fn run_benchmark_compiled_large_stress_multifield_aot_runtime() {
    let samples = stress_samples(7);
    let residual_iters = stress_residual_iters(80);
    let jacobian_iters = stress_jacobian_iters(32);
    let field_count = generated_stress_large_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_large_bvp_fixture_bindings::N_STEPS;

    println!(
        "\n=== Compiled larger stress multifield AOT runtime ===\nfields={}, steps={}",
        field_count, n_steps
    );
    println!(
        "samples={}, residual_iters={}, jacobian_iters={}",
        samples, residual_iters, jacobian_iters
    );

    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "fixture_stress_large_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_large_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 96,
    });
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 96 });

    let residual_bindings =
        bind_residual_fixture_chunks(&residual_plan, stress_large_residual_chunk_fns());
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, stress_large_sparse_chunk_fns());
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    let mut parallel2_residual_samples = Vec::with_capacity(samples);
    let mut sequential_residual_samples = Vec::with_capacity(samples);
    let mut parallel_residual_samples = Vec::with_capacity(samples);

    let mut parallel2_jacobian_samples = Vec::with_capacity(samples);
    let mut sequential_sparse_samples = Vec::with_capacity(samples);
    let mut parallel_sparse_samples = Vec::with_capacity(samples);

    for _ in 0..samples {
        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = jac
                .residiual_function
                .call(0.0, &input_col)
                .to_DVectorType();
            black_box(out[0]);
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let jacobian_fun = jac
            .jac_function
            .as_mut()
            .expect("parallel2 Jacobian should stay compiled");
        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
            black_box(out[(0, 0)]);
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), residual_iters)
    );
    println!(
        "sequential AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), residual_iters)
    );
    println!(
        "parallel AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), residual_iters)
    );
    println!(
        "parallel2 Jacobian: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), jacobian_iters)
    );
    println!(
        "sequential AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), jacobian_iters)
    );
    println!(
        "parallel AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), jacobian_iters)
    );
    println!(
        "compiled large stress plans: residual_chunks={}, sparse_chunks={}, jobs_residual={}, jobs_sparse={}",
        residual_plan.chunks.len(),
        sparse_plan.chunks.len(),
        residual_parallel.job_count(),
        sparse_parallel.job_count()
    );
}

#[test]
fn compiled_above_threshold_stress_fixture_matches_parallel2() {
    let field_count = generated_stress_above_threshold_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_above_threshold_bvp_fixture_bindings::N_STEPS;

    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "fixture_stress_above_threshold_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_above_threshold_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 96,
    });
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 96 });

    let residual_generated = eval_residual_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &residual_plan,
        stress_above_threshold_residual_chunk_fns(),
    );
    let residual_parallel2 = jac
        .residiual_function
        .call(0.0, &input_col)
        .to_DVectorType();
    for (lhs, rhs) in residual_parallel2.iter().zip(residual_generated.iter()) {
        assert!(
            (*lhs - *rhs).abs() < 1e-10,
            "above-threshold stress residual mismatch: {lhs} vs {rhs}"
        );
    }

    let jacobian_parallel2 = jac
        .jac_function
        .as_mut()
        .expect("parallel2 Jacobian should be compiled")
        .call(0.0, &input_col)
        .to_DMatrixType();
    let sparse_generated = eval_sparse_values_from_plan_and_chunks(
        flat_args.as_slice(),
        &sparse_plan,
        stress_above_threshold_sparse_chunk_fns(),
    );
    let jacobian_generated = sparse_plan
        .assemble_sparse_col_mat(sparse_generated.as_slice())
        .to_DMatrixType();
    for i in 0..jacobian_parallel2.nrows() {
        for j in 0..jacobian_parallel2.ncols() {
            let lhs = jacobian_parallel2[(i, j)];
            let rhs = jacobian_generated[(i, j)];
            assert!(
                (lhs - rhs).abs() < 1e-10,
                "above-threshold stress Jacobian mismatch at ({i}, {j}): {lhs} vs {rhs}"
            );
        }
    }
}

pub(crate) fn run_benchmark_compiled_above_threshold_stress_multifield_aot_runtime() {
    let samples = stress_samples(7);
    let residual_iters = stress_residual_iters(80);
    let jacobian_iters = stress_jacobian_iters(32);
    let field_count = generated_stress_above_threshold_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_above_threshold_bvp_fixture_bindings::N_STEPS;

    println!(
        "\n=== Compiled above-threshold stress multifield AOT runtime ===\nfields={}, steps={}",
        field_count, n_steps
    );
    println!(
        "samples={}, residual_iters={}, jacobian_iters={}",
        samples, residual_iters, jacobian_iters
    );

    let mut jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();
    jac.lambdify_residual_Col_parallel2("x", variable_names.clone());
    jac.lambdify_jacobian_SparseColMat_parallel2("x", variable_names.clone());

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();
    let input_col = ColRef::from_slice(flat_args.as_slice()).to_owned();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let residual_task = ResidualTask {
        fn_name: "fixture_stress_above_threshold_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_above_threshold_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 96,
    });
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 96 });

    let residual_bindings =
        bind_residual_fixture_chunks(&residual_plan, stress_above_threshold_residual_chunk_fns());
    let residual_sequential =
        SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let residual_parallel = ParallelResidualExecutor::with_config(
        &residual_plan,
        residual_bindings,
        benchmark_parallel_config(),
    );

    let sparse_bindings =
        bind_sparse_fixture_chunks(&sparse_plan, stress_above_threshold_sparse_chunk_fns());
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
        &sparse_plan,
        sparse_bindings,
        benchmark_parallel_config(),
    );

    let mut parallel2_residual_samples = Vec::with_capacity(samples);
    let mut sequential_residual_samples = Vec::with_capacity(samples);
    let mut parallel_residual_samples = Vec::with_capacity(samples);
    let mut parallel2_jacobian_samples = Vec::with_capacity(samples);
    let mut sequential_sparse_samples = Vec::with_capacity(samples);
    let mut parallel_sparse_samples = Vec::with_capacity(samples);

    for _ in 0..samples {
        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = jac
                .residiual_function
                .call(0.0, &input_col)
                .to_DVectorType();
            black_box(out[0]);
        }
        parallel2_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_sequential.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_residual_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..residual_iters {
            let out = residual_parallel.eval(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_residual_samples.push(start.elapsed());

        let jacobian_fun = jac
            .jac_function
            .as_mut()
            .expect("parallel2 Jacobian should stay compiled");
        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = jacobian_fun.call(0.0, &input_col).to_DMatrixType();
            black_box(out[(0, 0)]);
        }
        parallel2_jacobian_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_sequential.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        sequential_sparse_samples.push(start.elapsed());

        let start = Instant::now();
        for _ in 0..jacobian_iters {
            let out = sparse_parallel.eval_values(flat_args.as_slice());
            black_box(out[0]);
        }
        parallel_sparse_samples.push(start.elapsed());
    }

    println!(
        "parallel2 residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_residual_samples), residual_iters)
    );
    println!(
        "sequential AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_residual_samples), residual_iters)
    );
    println!(
        "parallel AOT residual: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_residual_samples), residual_iters)
    );
    println!(
        "parallel2 Jacobian: {:.2} ns/call",
        per_iter_ns(median_duration(parallel2_jacobian_samples), jacobian_iters)
    );
    println!(
        "sequential AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_sparse_samples), jacobian_iters)
    );
    println!(
        "parallel AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(parallel_sparse_samples), jacobian_iters)
    );
    println!(
        "compiled above-threshold stress plans: residual_chunks={}, sparse_chunks={}, jobs_residual={}, jobs_sparse={}",
        residual_plan.chunks.len(),
        sparse_plan.chunks.len(),
        residual_parallel.job_count(),
        sparse_parallel.job_count()
    );
}

pub(crate) fn run_benchmark_compiled_stress_sparse_parallel_job_counts() {
    let samples = stress_samples(9);
    let jacobian_iters = stress_jacobian_iters(48);
    let field_count = generated_stress_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_bvp_fixture_bindings::N_STEPS;

    println!(
        "\n=== Compiled stress sparse parallel job-count sweep ===\nfields={}, steps={}",
        field_count, n_steps
    );
    println!("samples={}, jacobian_iters={}", samples, jacobian_iters);

    let jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned
        .iter()
        .map(|name| name.as_str())
        .collect();

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|index| 0.05 + index as f64 * 0.0005)
        .collect();
    let runtime_args = RuntimeArguments::new(variable_values.as_slice(), None);
    let flat_args = runtime_args.flatten();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 48 });

    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, stress_sparse_chunk_fns());
    let sparse_sequential =
        SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());

    let sequential_samples: Vec<Duration> = (0..samples)
        .map(|_| {
            let start = Instant::now();
            for _ in 0..jacobian_iters {
                let out = sparse_sequential.eval_values(flat_args.as_slice());
                black_box(out[0]);
            }
            start.elapsed()
        })
        .collect();

    println!(
        "sequential AOT sparse values: {:.2} ns/call",
        per_iter_ns(median_duration(sequential_samples), jacobian_iters)
    );

    for sparse_jobs in [1usize, 2, 4] {
        let sparse_parallel = ParallelSparseJacobianExecutor::with_config(
            &sparse_plan,
            sparse_bindings.clone(),
            benchmark_sparse_parallel_config_with_jobs(sparse_jobs),
        );

        let samples: Vec<Duration> = (0..samples)
            .map(|_| {
                let start = Instant::now();
                for _ in 0..jacobian_iters {
                    let out = sparse_parallel.eval_values(flat_args.as_slice());
                    black_box(out[0]);
                }
                start.elapsed()
            })
            .collect();

        println!(
            "parallel AOT sparse values (requested_jobs={}, actual_jobs={}, fallback={}): {:.2} ns/call",
            sparse_jobs,
            sparse_parallel.job_count(),
            sparse_parallel.uses_sequential_fallback(),
            per_iter_ns(median_duration(samples), jacobian_iters)
        );
    }
}

// ============================================================================
// PARALLEL-NO-BENEFIT DIAGNOSTICS
//
// Three root causes are probed:
//
//  1. RAYON OVERHEAD BASELINE
//     Measures raw rayon::join cost with no-op closures.
//     If sparse eval time is within ~5x of this, parallelism cannot help
//     at this problem size regardless of chunk count.
//
//  2. CHUNK GRANULARITY + FALLBACK
//     Sweeps requested_jobs 1..=2*threads on the compiled stress fixture.
//     Prints actual_jobs, uses_sequential_fallback(), nnz/job, and ns/call.
//     Exposes whether the Auto fallback is silently killing the parallel path
//     and whether each job has enough work to amortize scheduling cost.
//
//  3. PROBLEM-SIZE CROSSOVER (IR-level, no compiled fixture needed)
//     Sweeps field_count at fixed n_steps using LinearBlock IR evaluation.
//     Prints the first field_count where parallel IR beats sequential IR.
//     Answers: "is the problem simply too small?"
// ============================================================================

/// Measures raw rayon::join overhead with no-op closures.
/// Run with: cargo test diagnose_rayon_overhead -- --nocapture --ignored
/// Measures raw rayon::join cost with no-op closures (2-way and 4-way nested).
/// If the sparse eval time per call is within ~5× of this overhead, parallelism is
/// structurally impossible to help at that problem size — the scheduler costs more than the math.
pub(crate) fn run_diagnose_rayon_overhead_baseline() {
    const ITERS: usize = 10_000;
    const SAMPLES: usize = 9;
    let workers = rayon::current_num_threads();

    let seq_timings: Vec<Duration> = (0..SAMPLES)
        .map(|_| {
            let t = Instant::now();
            for _ in 0..ITERS {
                black_box(0u64);
            }
            t.elapsed()
        })
        .collect();

    let join2_timings: Vec<Duration> = (0..SAMPLES)
        .map(|_| {
            let t = Instant::now();
            for _ in 0..ITERS {
                rayon::join(|| black_box(0u64), || black_box(0u64));
            }
            t.elapsed()
        })
        .collect();

    let join4_timings: Vec<Duration> = (0..SAMPLES)
        .map(|_| {
            let t = Instant::now();
            for _ in 0..ITERS {
                rayon::join(
                    || rayon::join(|| black_box(0u64), || black_box(0u64)),
                    || rayon::join(|| black_box(0u64), || black_box(0u64)),
                );
            }
            t.elapsed()
        })
        .collect();

    let seq_ns = per_iter_ns(median_duration(seq_timings), ITERS);
    let join2_ns = per_iter_ns(median_duration(join2_timings), ITERS);
    let join4_ns = per_iter_ns(median_duration(join4_timings), ITERS);

    println!(
        "\n=== Rayon overhead baseline (workers={workers}, iters={ITERS}) ===\
         \nno-op sequential:          {seq_ns:.1} ns/call\
         \nno-op join (2-way):        {join2_ns:.1} ns/call  (overhead +{:.1} ns)\
         \nno-op join (4-way nested): {join4_ns:.1} ns/call  (overhead +{:.1} ns)\
         \n\nIf sparse eval time is within ~5x of join overhead, parallelism\
         \ncannot help at this problem size regardless of chunk count.",
        join2_ns - seq_ns,
        join4_ns - seq_ns,
    );
}

/// Sweeps requested_jobs on the compiled stress fixture and prints Auto-policy
/// diagnostics plus forced-parallel timings for comparison.
/// Run with: cargo test diagnose_chunk_granularity -- --nocapture --ignored
/// diagnose_chunk_granularity_and_fallback
///Uses the compiled stress fixture (8 fields × 96 steps) and sweeps requested_jobs from 1 to 2 × thread_count.
/// For each configuration it prints:
/// actual_jobs — how many jobs the weighted splitter actually produced (capped by chunk count)
///uses_sequential_fallback() — whether the Auto policy silently fell back even with Never override
/// nnz/job — work per job; if this is < ~50, rayon overhead dominates
/// ns/call — actual timing with <-- faster marker
pub(crate) fn run_diagnose_chunk_granularity_and_fallback() {
    let samples = stress_samples(9);
    let jacobian_iters = stress_jacobian_iters(200);
    let max_threads = rayon::current_num_threads();

    let field_count = generated_stress_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_bvp_fixture_bindings::N_STEPS;

    let jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned.iter().map(|s| s.as_str()).collect();
    let flat_args = RuntimeArguments::new(
        &(0..jac.variable_string.len())
            .map(|i| 0.05 + i as f64 * 0.0005)
            .collect::<Vec<_>>(),
        None,
    )
    .flatten();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
    let sparse_task = SparseJacobianTask {
        fn_name: "diag_sparse",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };
    // Same chunking as the compiled stress fixture.
    let sparse_plan =
        sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount { rows_per_chunk: 48 });
    let nnz = sparse_plan.nnz();
    let chunk_count = sparse_plan.chunks.len();
    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, stress_sparse_chunk_fns());
    let auto_probe = ParallelSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let min_work = auto_probe.min_work_per_parallel_job();
    let nnz_per_chunk = auto_probe.work_per_chunk();

    println!(
        "\n=== Chunk granularity / fallback diagnostic ===\
         \nfields={field_count}, steps={n_steps}, nnz={nnz}, chunks={chunk_count}\
         \nrayon_threads={max_threads}, samples={samples}, jacobian_iters={jacobian_iters}\
         \nAuto policy derives min_work/job from measured rayon overhead.\
         \nmin_work/job={min_work}, nnz/chunk={nnz_per_chunk}\n"
    );
    println!(
        "{:<16} {:<12} {:<10} {:<10} {:<10} {:<12}",
        "mode", "jobs", "fallback", "nnz/job", "nnz/chunk", "ns/call"
    );

    // Sequential baseline.
    let seq = SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());
    let seq_ns = per_iter_ns(
        median_duration(
            (0..samples)
                .map(|_| {
                    let t = Instant::now();
                    for _ in 0..jacobian_iters {
                        black_box(seq.eval_values(&flat_args)[0]);
                    }
                    t.elapsed()
                })
                .collect(),
        ),
        jacobian_iters,
    );
    println!(
        "{:<16} {:<12} {:<10} {:<10} {:<10} {:.2}",
        "sequential", "-", "-", nnz, "-", seq_ns
    );

    let auto_ns = per_iter_ns(
        median_duration(
            (0..samples)
                .map(|_| {
                    let t = Instant::now();
                    for _ in 0..jacobian_iters {
                        black_box(auto_probe.eval_values(&flat_args)[0]);
                    }
                    t.elapsed()
                })
                .collect(),
        ),
        jacobian_iters,
    );
    println!(
        "{:<16} {:<12} {:<10} {:<10} {:<10} {:.2}  {}",
        "auto",
        auto_probe.job_count(),
        auto_probe.uses_sequential_fallback(),
        auto_probe.work_per_job(),
        auto_probe.work_per_chunk(),
        auto_ns,
        if auto_ns < seq_ns { "<-- faster" } else { "" }
    );

    println!("\nForced sweep (`ParallelFallbackPolicy::Never`) for diagnostics:");
    println!(
        "{:<16} {:<12} {:<10} {:<10} {:<10} {:<12}",
        "requested", "actual", "fallback", "nnz/job", "nnz/chunk", "ns/call"
    );

    // Parallel sweep with Never fallback so we always exercise the parallel path.
    for requested_jobs in 1..=(2 * max_threads).min(chunk_count * 2).max(4) {
        let config = ParallelExecutorConfig {
            jobs_per_worker: 1,
            max_residual_jobs: None,
            max_sparse_jobs: Some(requested_jobs),
            fallback_policy: ParallelFallbackPolicy::Never,
        };
        let par = ParallelSparseJacobianExecutor::with_config(
            &sparse_plan,
            sparse_bindings.clone(),
            config,
        );
        let ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..jacobian_iters {
                            black_box(par.eval_values(&flat_args)[0]);
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            jacobian_iters,
        );
        println!(
            "{:<16} {:<12} {:<10} {:<10} {:<10} {:.2}  {}",
            requested_jobs,
            par.job_count(),
            par.uses_sequential_fallback(),
            par.work_per_job(),
            par.work_per_chunk(),
            ns,
            if ns < seq_ns { "<-- faster" } else { "" }
        );
    }
}

#[derive(Clone, Debug)]
struct IrAmplificationRow {
    matrix: &'static str,
    stage: &'static str,
    strategy: String,
    chunks: usize,
    outputs: usize,
    instructions: usize,
    max_block_instructions: usize,
    temps: usize,
    source_kb: f64,
}

fn ir_amplification_row(
    matrix: &'static str,
    stage: &'static str,
    strategy: impl Into<String>,
    chunks: usize,
    blocks: Vec<(usize, LinearBlock)>,
    source_len: usize,
) -> IrAmplificationRow {
    let mut outputs = 0usize;
    let mut instructions = 0usize;
    let mut max_block_instructions = 0usize;
    let mut temps = 0usize;
    for (_, block) in blocks {
        outputs += block.outputs.len();
        instructions += block.instructions.len();
        max_block_instructions = max_block_instructions.max(block.instructions.len());
        temps += block.num_temps;
    }

    IrAmplificationRow {
        matrix,
        stage,
        strategy: strategy.into(),
        chunks,
        outputs,
        instructions,
        max_block_instructions,
        temps,
        source_kb: source_len as f64 / 1024.0,
    }
}

fn ir_amplification_baselines(
    rows: &[IrAmplificationRow],
) -> HashMap<(&'static str, &'static str), usize> {
    rows.iter()
        .filter(|row| row.strategy == "whole")
        .map(|row| ((row.matrix, row.stage), row.instructions.max(1)))
        .collect()
}

/// Prints whether generated chunking inflates IR/source work before runtime
/// scheduling and FFI even start.
pub(crate) fn run_diagnose_combustion_chunk_ir_amplification() {
    let n_steps = env_usize("BVP_CHUNK_IR_STEPS", 1000);
    let target_chunks = [
        ("whole", 1usize),
        ("chunk4", 4usize),
        ("chunk8", 8usize),
        ("chunk16", 16usize),
    ];

    let jac =
        build_combustion_bvp_case_with_backend(n_steps, BvpSymbolicAssemblyBackend::AtomView);
    let mut rows = Vec::new();

    for (label, chunks) in target_chunks {
        let residual_strategy = if chunks == 1 {
            ResidualChunkingStrategy::Whole
        } else {
            ResidualChunkingStrategy::ByTargetChunkCount {
                target_chunks: chunks,
            }
        };
        let sparse_strategy = if chunks == 1 {
            SparseChunkingStrategy::Whole
        } else {
            SparseChunkingStrategy::ByTargetChunkCount {
                target_chunks: chunks,
            }
        };

        let prepared = jac.prepare_sparse_aot_problem(
            "diagnose_combustion_residual",
            "diagnose_combustion_jacobian_values",
            residual_strategy,
            sparse_strategy,
        );
        let sparse_prepared = prepared.as_prepared_problem();
        let banded_prepared = prepared.as_prepared_banded_problem();

        rows.push(ir_amplification_row(
            "Shared",
            "Residual",
            label,
            sparse_prepared.residual_plan.chunks.len(),
            residual_ir_blocks(&sparse_prepared.residual_plan),
            residual_module("diag_residual", &sparse_prepared.residual_plan)
                .emit_source()
                .len(),
        ));
        rows.push(ir_amplification_row(
            "Sparse",
            "Jacobian",
            label,
            sparse_prepared.jacobian_plan.chunks.len(),
            sparse_ir_blocks(&sparse_prepared.jacobian_plan),
            sparse_module("diag_sparse", &sparse_prepared.jacobian_plan)
                .emit_source()
                .len(),
        ));
        rows.push(ir_amplification_row(
            "Banded",
            "Jacobian",
            label,
            banded_prepared.jacobian_plan.chunks.len(),
            banded_ir_blocks(&banded_prepared.jacobian_plan),
            banded_module("diag_banded", &banded_prepared.jacobian_plan)
                .emit_source()
                .len(),
        ));
    }

    let baselines = ir_amplification_baselines(&rows);
    println!(
        "\n=== Combustion BVP chunk IR/source amplification ===\
         \nn_steps={n_steps}, unknowns={}, residuals={}, sparse_nnz={}\
         \nThis does not execute callbacks. It measures generated IR/source size before scheduling/FFI.",
        jac.vector_of_variables.len(),
        jac.vector_of_functions.len(),
        jac.symbolic_jacobian_sparse_entries_owned().len(),
    );
    println!(
        "{:<8} | {:<8} | {:<7} | {:>6} | {:>8} | {:>10} | {:>8} | {:>8} | {:>9} | {:>8}",
        "matrix",
        "stage",
        "strategy",
        "chunks",
        "outputs",
        "instr",
        "amp",
        "max_blk",
        "temps",
        "src_kb",
    );
    println!("{}", "-".repeat(106));
    for row in &rows {
        let baseline = *baselines
            .get(&(row.matrix, row.stage))
            .expect("whole baseline must exist for every matrix/stage");
        let amp = row.instructions as f64 / baseline as f64;
        println!(
            "{:<8} | {:<8} | {:<7} | {:>6} | {:>8} | {:>10} | {:>8.3} | {:>8} | {:>9} | {:>8.1}",
            row.matrix,
            row.stage,
            row.strategy,
            row.chunks,
            row.outputs,
            row.instructions,
            amp,
            row.max_block_instructions,
            row.temps,
            row.source_kb,
        );
    }

    println!(
        "\nInterpretation: amp > 1 means chunking generated more total straight-line IR than the whole callback. \
         If that happens, runtime parallelism must first pay back duplicated arithmetic/source work plus FFI/scheduler overhead."
    );
}
/*
rayon::join costs ~2000 ns per call

Sequential eval costs ~16900 ns — only 8× the join overhead

With 16 chunks and 20 workers, each job gets ~191 values — not enough work

Even at field_count=2 (nnz=573), parallel IR already wins — meaning the problem is the compiled
 chunk granularity , not the problem size itself
The root cause: the compiled stress fixture uses rows_per_chunk=48, producing only 16 chunks for 3056 nnz.
 With 20 workers, jobs get ~191 values each — too little to amortize the 2000 ns join cost.
 The crossover test shows parallel wins even at small sizes when chunks are sized correctly via IR.

The fix is to increase chunk granularity so each job gets enough work. Let me check what chunk
count would break even:

*/

/// IR-level field_count sweep: finds the first problem size where parallel IR
/// beats sequential IR for sparse Jacobian values. No compiled fixture needed.
/// Run with: cargo test diagnose_problem_size_crossover -- --nocapture --ignored
/// IR-level sweep over field_count ∈ [2,4,6,8,12,16,24,32,48,64] at n_steps=96.
/// Uses LinearBlock IR evaluation (no compiled fixture needed for each size) and prints the
/// first field_count where parallel IR beats sequential IR. Answers the question:
/// is the problem simply too small for any parallelism to help?
pub(crate) fn run_diagnose_problem_size_crossover() {
    let samples = stress_samples(7);
    let jacobian_iters = stress_jacobian_iters(100);
    let n_steps = 96usize;
    let field_counts = stress_env_list_usize(
        "RST_DIAG_FIELD_COUNTS",
        &[2, 4, 6, 8, 12, 16, 24, 32, 48, 64],
    );

    println!(
        "\n=== Problem-size crossover sweep (n_steps={n_steps}) ===\
         \nsamples={samples}, jacobian_iters={jacobian_iters}\n"
    );
    println!(
        "{:<8} {:<10} {:<8} {:<14} {:<14} {:<10}",
        "fields", "unknowns", "nnz", "seq_ns/call", "par_ns/call", "par_wins"
    );

    let mut first_crossover: Option<usize> = None;

    for field_count in field_counts {
        let jac = build_multifield_stress_bvp_case(field_count, n_steps);
        let variable_names_owned = jac.variable_string.clone();
        let variable_names: Vec<&str> = variable_names_owned.iter().map(|s| s.as_str()).collect();
        let flat_args = RuntimeArguments::new(
            &(0..jac.variable_string.len())
                .map(|i| 0.05 + i as f64 * 0.0005)
                .collect::<Vec<_>>(),
            None,
        )
        .flatten();

        let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
        let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);
        let sparse_task = SparseJacobianTask {
            fn_name: "diag_crossover",
            shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
            entries: &sparse_entries,
            variables: &variable_names,
            params: None,
        };

        // 2 chunks per worker gives rayon room to load-balance.
        let strategy = recommended_row_chunking_for_parallelism(jac.vector_of_functions.len(), 2);
        let sparse_plan = sparse_task.runtime_plan(strategy);
        let nnz = sparse_plan.nnz();
        let unknowns = jac.variable_string.len();

        let ir_blocks = sparse_ir_blocks(&sparse_plan);

        let seq_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..jacobian_iters {
                            black_box(eval_ir_blocks_sequential(&flat_args, &ir_blocks, nnz)[0]);
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            jacobian_iters,
        );

        let par_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..jacobian_iters {
                            black_box(eval_ir_blocks_parallel(&flat_args, &ir_blocks, nnz)[0]);
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            jacobian_iters,
        );

        let par_wins = par_ns < seq_ns;
        if par_wins && first_crossover.is_none() {
            first_crossover = Some(field_count);
        }

        println!(
            "{:<8} {:<10} {:<8} {:<14.1} {:<14.1} {}",
            field_count,
            unknowns,
            nnz,
            seq_ns,
            par_ns,
            if par_wins { "YES <--" } else { "no" }
        );
    }

    println!(
        "\nFirst crossover: {}",
        first_crossover
            .map(|f| format!("field_count={f}"))
            .unwrap_or_else(|| "not found in sweep — increase field_counts range".to_string())
    );
}

/*
 Here's the diagnosis:

Root cause: too many jobs relative to work per job

The numbers are unambiguous:

fact	value
rayon::join overhead	~2000 ns
sequential eval	~16900 ns
ratio	~8× — barely above the "cannot help" threshold
chunks in compiled fixture	16 (from rows_per_chunk=48)
nnz per job at 16 jobs	191
best parallel result	~22700 ns (job=8, nnz/job=382) — still 34% slower than sequential
With 20 workers and jobs_per_worker=1, the executor tries to create 20 jobs from 16 chunks —
 it caps at 16. Each job then does ~191 values of math, but costs ~2000 ns just to schedule.
  The math is too cheap relative to the scheduler.

The crossover test confirms parallelism can work at this problem size — it wins at field_count=2 (nnz=573)
 with IR blocks. The difference is that IR blocks are sized by the problem, not by a fixed rows_per_chunk=48 that produces too many tiny chunks.

Three concrete fixes:

1. Raise DEFAULT_MIN_SPARSE_VALUES_FOR_PARALLEL in codegen_orchestrator.rs

The current threshold of 256 is too low. With 2000 ns join overhead and ~5 ns per value, you need at
least ~2000 values total before parallelism can break even:
1) instead of ParallelFallbackPolicy::Auto => output_len < DEFAULT_MIN_RESIDUAL_OUTPUTS_FOR_PARALLEL,
"
 ParallelFallbackPolicy::Auto => {
                // Fall back if total work is below threshold OR each job gets
                // too little work to amortize rayon scheduling overhead.
                // At ~2000 ns/join and ~5 ns/value, a job needs ~400 values.
                nnz < DEFAULT_MIN_SPARSE_VALUES_FOR_PARALLEL
                    || nnz / job_count < 400
            }""
2) instead of const DEFAULT_MIN_RESIDUAL_OUTPUTS_FOR_PARALLEL: usize = 128;
const DEFAULT_MIN_SPARSE_VALUES_FOR_PARALLEL: usize = 256;
"
// Minimum work thresholds before the Auto policy allows parallel execution.
// Derived from measured rayon::join overhead (~2000 ns on this machine) and
// typical per-value cost (~5 ns): need at least ~400 values/job to break even,
// times a minimum of 2 jobs = ~800 values total. Use 2000 as a conservative
// floor that accounts for variance and non-uniform chunk costs.
const DEFAULT_MIN_RESIDUAL_OUTPUTS_FOR_PARALLEL: usize = 512;
const DEFAULT_MIN_SPARSE_VALUES_FOR_PARALLEL: usize = 2000;"
*/

// ============================================================================
// XLARGE FIXTURE: sequential AOT vs parallel AOT only
// field_count=32, n_steps=256, ~8192 unknowns, nnz=32704, sparse_chunks=11
// ============================================================================

use crate::symbolic::codegen::testing_fixtures::test_codegen_generated_stress_xlarge_bvp_fixtures::generated_stress_xlarge_bvp_fixture_bindings;

fn xlarge_residual_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_xlarge_bvp_fixture_bindings::RESIDUAL_CHUNKS
}

fn xlarge_sparse_chunk_fns() -> &'static [GeneratedChunkFn] {
    &generated_stress_xlarge_bvp_fixture_bindings::SPARSE_VALUE_CHUNKS
}

/// Sequential AOT vs parallel AOT on the xlarge fixture.
///
/// This is the first fixture large enough that parallel should clearly win:
/// - 8192 unknowns, nnz=32704, sparse_chunks=11
/// - each sparse chunk has ~3k nnz × ~5 ns = ~15 µs of math
/// - rayon join overhead ~2 µs → 7.5× ratio, well above break-even
///
/// Run with: cargo test benchmark_compiled_xlarge_stress_seq_vs_par -- --nocapture --ignored
pub(crate) fn run_benchmark_compiled_xlarge_stress_seq_vs_par() {
    let samples = stress_samples(7);
    let residual_iters = stress_residual_iters(50);
    let jacobian_iters = stress_jacobian_iters(20);

    let field_count = generated_stress_xlarge_bvp_fixture_bindings::FIELD_COUNT;
    let n_steps = generated_stress_xlarge_bvp_fixture_bindings::N_STEPS;
    let nnz = generated_stress_xlarge_bvp_fixture_bindings::NNZ;
    let n_sparse_chunks = generated_stress_xlarge_bvp_fixture_bindings::N_SPARSE_CHUNKS;
    let n_residual_chunks = generated_stress_xlarge_bvp_fixture_bindings::N_RESIDUAL_CHUNKS;

    println!(
        "\n=== XLarge compiled sparse AOT: sequential vs parallel ===\
         \nfields={field_count}, steps={n_steps}, unknowns={}, nnz={nnz}\
         \nresidual_chunks={n_residual_chunks}, sparse_chunks={n_sparse_chunks}\
         \nsamples={samples}, residual_iters={residual_iters}, jacobian_iters={jacobian_iters}",
        field_count * n_steps
    );

    let jac = build_multifield_stress_bvp_case(field_count, n_steps);
    let variable_names_owned = jac.variable_string.clone();
    let variable_names: Vec<&str> = variable_names_owned.iter().map(|s| s.as_str()).collect();

    let variable_values: Vec<f64> = (0..jac.variable_string.len())
        .map(|i| 0.05 + i as f64 * 0.0005)
        .collect();
    let flat_args = RuntimeArguments::new(&variable_values, None).flatten();

    let sparse_entries_owned = jac.symbolic_jacobian_sparse_entries_owned();
    let sparse_entries = borrowed_sparse_entries(&sparse_entries_owned);

    let residual_task = ResidualTask {
        fn_name: "fixture_stress_xlarge_residual",
        residuals: &jac.vector_of_functions,
        variables: &variable_names,
        params: None,
    };
    let sparse_task = SparseJacobianTask {
        fn_name: "fixture_stress_xlarge_sparse_values",
        shape: (jac.vector_of_functions.len(), jac.vector_of_variables.len()),
        entries: &sparse_entries,
        variables: &variable_names,
        params: None,
    };

    // Use the same chunking the fixture was generated with.
    let residual_plan = residual_task.runtime_plan(ResidualChunkingStrategy::ByOutputCount {
        max_outputs_per_chunk: 256,
    });
    let sparse_plan = sparse_task.runtime_plan(SparseChunkingStrategy::ByRowCount {
        rows_per_chunk: 800,
    });

    let residual_bindings =
        bind_residual_fixture_chunks(&residual_plan, xlarge_residual_chunk_fns());
    let sparse_bindings = bind_sparse_fixture_chunks(&sparse_plan, xlarge_sparse_chunk_fns());

    let residual_seq = SequentialResidualExecutor::new(&residual_plan, residual_bindings.clone());
    let sparse_seq = SequentialSparseJacobianExecutor::new(&sparse_plan, sparse_bindings.clone());

    // Correctness check before timing.
    assert_eq!(
        residual_seq.eval(&flat_args).len(),
        residual_plan.output_len
    );
    assert_eq!(sparse_seq.eval_values(&flat_args).len(), sparse_plan.nnz());

    // Sweep job counts: 1 (=sequential), 2, 4, n_sparse_chunks.
    let job_counts = [1usize, 2, 4, n_sparse_chunks];

    println!(
        "\n{:<16} {:<12} {:<12} {:<14} {:<14}",
        "config", "res_jobs", "sparse_jobs", "res_ns/call", "sparse_ns/call"
    );

    // Sequential baseline.
    let seq_res_ns = per_iter_ns(
        median_duration(
            (0..samples)
                .map(|_| {
                    let t = Instant::now();
                    for _ in 0..residual_iters {
                        black_box(residual_seq.eval(&flat_args)[0]);
                    }
                    t.elapsed()
                })
                .collect(),
        ),
        residual_iters,
    );
    let seq_sparse_ns = per_iter_ns(
        median_duration(
            (0..samples)
                .map(|_| {
                    let t = Instant::now();
                    for _ in 0..jacobian_iters {
                        black_box(sparse_seq.eval_values(&flat_args)[0]);
                    }
                    t.elapsed()
                })
                .collect(),
        ),
        jacobian_iters,
    );
    println!(
        "{:<16} {:<12} {:<12} {:<14.0} {:<14.0}",
        "sequential", "-", "-", seq_res_ns, seq_sparse_ns
    );

    for &jobs in &job_counts {
        let config = ParallelExecutorConfig {
            jobs_per_worker: 1,
            max_residual_jobs: Some(jobs),
            max_sparse_jobs: Some(jobs),
            fallback_policy: ParallelFallbackPolicy::Never,
        };
        let res_par = ParallelResidualExecutor::with_config(
            &residual_plan,
            residual_bindings.clone(),
            config,
        );
        let sparse_par = ParallelSparseJacobianExecutor::with_config(
            &sparse_plan,
            sparse_bindings.clone(),
            config,
        );

        // Correctness.
        assert_eq!(
            res_par.eval(&flat_args),
            residual_seq.eval(&flat_args),
            "parallel residual mismatch at jobs={jobs}"
        );
        assert_eq!(
            sparse_par.eval_values(&flat_args),
            sparse_seq.eval_values(&flat_args),
            "parallel sparse mismatch at jobs={jobs}"
        );

        let res_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..residual_iters {
                            black_box(res_par.eval(&flat_args)[0]);
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            residual_iters,
        );
        let sparse_ns = per_iter_ns(
            median_duration(
                (0..samples)
                    .map(|_| {
                        let t = Instant::now();
                        for _ in 0..jacobian_iters {
                            black_box(sparse_par.eval_values(&flat_args)[0]);
                        }
                        t.elapsed()
                    })
                    .collect(),
            ),
            jacobian_iters,
        );

        println!(
            "{:<16} {:<12} {:<12} {:<14.0} {:<14.0}  res_speedup={:.2}x  sparse_speedup={:.2}x  (fallback: res={}, sparse={})",
            format!("par-jobs{jobs}"),
            res_par.job_count(),
            sparse_par.job_count(),
            res_ns,
            sparse_ns,
            seq_res_ns / res_ns,
            seq_sparse_ns / sparse_ns,
            res_par.uses_sequential_fallback(),
            sparse_par.uses_sequential_fallback(),
        );
    }
}