RustedSciThe 0.3.35

#![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::CodegenIR::{CodegenModule, LinearBlock};
use crate::symbolic::codegen_adapters::{
    residual_and_sparse_module, residual_ir_blocks, sparse_ir_blocks,
};
use crate::symbolic::codegen_orchestrator::{
    GeneratedChunkFn, ParallelResidualExecutor, ParallelSparseJacobianExecutor,
    ResidualChunkBinding, SequentialResidualExecutor, SequentialSparseJacobianExecutor,
    SparseJacobianChunkBinding,
};
use crate::symbolic::codegen_runtime_api::{
    ResidualChunkingStrategy, RuntimeArguments, recommended_residual_chunking_for_parallelism,
    recommended_row_chunking_for_parallelism,
};
use crate::symbolic::codegen_tasks::{
    ResidualTask, SparseChunkingStrategy, SparseExprEntry, SparseJacobianTask,
};
use crate::symbolic::codegen_test_support::{
    benchmark_parallel_config, benchmark_sparse_parallel_config_with_jobs,
    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::symbolic::symbolic_engine::Expr;
use crate::symbolic::symbolic_functions_BVP::Jacobian;
use crate::symbolic::test_codegen_generated_bvp_fixtures::generated_bvp_fixture;
use crate::symbolic::test_codegen_generated_bvp_huge_fixtures::generated_bvp_huge_fixture_bindings;
use crate::symbolic::test_codegen_generated_bvp_large_chunk_variants::generated_bvp_large_chunk_variant_bindings;
use crate::symbolic::test_codegen_generated_bvp_large_fixtures::generated_bvp_large_fixture_bindings;
use crate::symbolic::test_codegen_generated_bvp_large_superblock_variants::generated_bvp_large_superblock_variant_bindings;
use crate::symbolic::test_codegen_generated_stress_above_threshold_bvp_fixtures::generated_stress_above_threshold_bvp_fixture_bindings;
use crate::symbolic::test_codegen_generated_stress_bvp_fixtures::generated_stress_bvp_fixture_bindings;
use crate::symbolic::test_codegen_generated_stress_large_bvp_fixtures::generated_stress_large_bvp_fixture_bindings;
use faer::col::{Col, ColRef};
use rayon::prelude::*;
use std::collections::HashMap;
use std::hint::black_box;
use std::time::{Duration, Instant};

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
}

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 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_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_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_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_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,
        );
    }
}

#[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),
    );
}

#[test]
#[ignore = "stress-scale diagnostic; run explicitly with optional RST_STRESS_FIELDS / RST_STRESS_STEPS"]
fn 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()
    );
}

#[test]
#[ignore = "stress-scale field-count sweep for larger multifield systems"]
/// 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?".
fn 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()
        );
    }
}

#[test]
#[ignore = "stress-scale crossover search for sequential vs parallel AOT IR"]
fn 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())
    );
}

#[test]
#[ignore = "stress-scale compiled AOT runtime diagnostic"]
fn 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}"
            );
        }
    }
}

#[test]
#[ignore = "compiled larger stress AOT runtime diagnostic"]
fn 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}"
            );
        }
    }
}

#[test]
#[ignore = "compiled above-threshold stress AOT runtime diagnostic"]
fn 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()
    );
}

#[test]
#[ignore = "stress-scale sparse parallel job sweep"]
fn 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)
        );
    }
}